dankelley
Meteorology Display
Discussed in https://github.com/dankelley/oce/discussions/2190
Originally posted by mdupilka February 17, 2024 Hi, I use the OCE package regularly to map meteorology fields for my work. Do you have any thoughts of possibly extending OCE to display meteorology data such as wind fields using wind barbs? Regards, Max
Hi. First, I made this into an issue, not a discussion. That's because I regularly look at issues, as a sort of to-do list. Also, I'm tagging @richardsc and @clayton33, who might have thoughts on whether this would be useful.
As for your question, that's something that I've thought of, but wasn't sure there would be enough interest. I am always thinking of new ways to plot data, and usually try them outside oce to start (see e.g. https://dankelley.github.io/dek_blog/2024/02/03/scotian-shelf-shaded.html) because then testing is a lot faster. (Building oce takes about 5 minutes on my machine, for example, and lots of users don't even have the C compilers, etc. required to build oce from source.)
Question: is there a source that you consider definitive, for how to draw barbs? For example, I just searched and found a top hit as https://www.weather.gov/hfo/windbarbinfo but that is in Empirical units (mps and knots) and I wonder if that is what is used outside the US, where km/h is likely used. I would need a source that is stable, non-commercial, and of course that does not need registration or the answering of questions about cookies. (CRAN does not permit using any such web citations. Packages that cite websources that cannot be reached by their algorithm get rejected. That's quite a lot of sites that you can visit easily ... it has something to do with whether sites can answer certain queries, not whether you can see them in a browser.)
If you can provide some sources for me to read, that would help. I like having something I can cite in the help for a function. For example, I might say "see page 145 of https://tc.canada.ca/sites/default/files/2020-10/aim-2020-2_met-e.pdf", which is something from the Transport Canada federal agency. Also, if I can get a list of sources I can see variants. Some will be much faster to plot than others. If it can be done with straight lines, for example, that will be a lot faster. I think it will be a slow function anyway, but waiting 20 seconds or so for a plot that is desired is better than waiting 0.5 seconds for something that is unwanted.
PS. you're aware of https://dankelley.github.io/oce/reference/drawDirectionField.html I assume? The last diagram you see there might be a start. I realize it's not the same, of course. This is how ocean people usually like to plot flow fields. But there is the problem of overlap, and that is alleviated by barbs, which are of equal length. That is where I think this could be useful.
PPS. it's too bad you didn't ask before Christmas. Last term I was teaching an oceanography class in which half the students were meteorologists. They could tell me in a flash things like (a) whether definitions vary across nations, (b) which is preferred in today's practice and (c) are there variants that would make sense to them, but not be expensive to plot. (This last is because things with line segments will be faster than things with curved segments. Faster by a lot.)
PS. I found the following but it's from a site that we cannot use as a source (it's not open enough) and it says that it comes from another site that doesn't even exist.
https://www.weather.gov/hfo/windbarbinfo seems like a source that could be cited and is likely definitive (for an audience in the USA, anyway ... I still do wonder whether knots are used in all nations.) Below is a snapshot.
So far, I see three notations for the starting location: show open circle, show filled circle or show no circle.
Interestingly, I see in the below that the angles are quantized. I never realized that before. This is fun.
Um, the image in the just-previous comment lacks barb examples for 40 and 45 knots. So it is not really rounded to nearest 5 knots, as they state. Surely something is supposed to be plotted in those categories. Can you shed light on this, @mdupilka? I imagine there is a better source out there somewhere.
Another source http://ww2010.atmos.uiuc.edu/(Gh)/guides/maps/sfcobs/wnd.rxml is similar (apart from the ugly starting point) but it lacks some categories.
I see that matplotlib has something (https://matplotlib.org/stable/gallery/images_contours_and_fields/barb_demo.html) and my guess is that some users will be familiar with that notation. I'll look into this, whilst waiting to see if others in this thread have any suggestions.
My main concerns right now are
- what to do in the speed gaps that I see in the examples?
- which of many conventions to follow
- what is a (CRAN-acceptable) source that can be cited in docs?
Below is a reprex of an initial test. If you rerun it, you'll get different results because it uses a random number. This just tests
- setting to uniform length
- discretization of angles on a "36 point compass", which I see in one of the sources noted above
- categorization into barb category (which is an unsettled issue given the speed-gaps alluded to previously)
This draws just 1 barb, and the feathers are not yet drawn. I would like to have clarity on the speed gaps before doing that. I think the scheme is a bit like "making change" on a cash register. Put another way, it's a bit like Roman numerals. (I realize that making change and Roman numerals are historical artifacts, but so are these weather maps, I think.)
Please click the word below, to see this test. (@richardsc and @clayton33 will know that I am just doing my usual thing, of starting with something crude and simple. I do lots of little modifications of tests like this. Only when things are pretty clear do I consider putting code into oce.)
met <- TRUE # set TRUE for barbs pointing to wind source
scale <- 0.5
debug <- TRUE
par(mar = c(3, 3, 4, 2)) # using top margin for notes
plot(c(-1, 1), c(-1, 1), xlab = "", ylab = "", asp = 1, type = "n")
x0 <- 0
y0 <- 0
u <- rnorm(1, sd = 10)
v <- rnorm(1, sd = 10)
uOrig <- u
vOrig <- v
if (met) {
u <- -u
v <- -v
}
angle <- 180 / pi * atan2(v, u)
# Barbs are drawn at angles on a 36-point compass
angle36 <- 10 * round(angle / 10)
# Barbs are defined in terms of knots, rounded to nearest 5
speed <- sqrt(u^2 + v^2)
knot <- speed * 1.94384
knot5 <- 5L * as.integer(round(knot / 5))
barbType <- 1 + knot5 / 5
points(x0, y0)
S <- sinpi(angle36 / 180)
C <- cospi(angle36 / 180)
speed
#> [1] 7.438785
x1 <- x0 + scale * S
y1 <- y0 + scale * C
# guiding circles to test angle and length
if (debug) {
circlex <- scale * cos(pi / 180 * seq(0, 360, 10))
circley <- scale * sin(pi / 180 * seq(0, 360, 10))
lines(circlex, circley, type = "o", pch = 20, col = 4)
}
notStill <- knot5 != 0
segments(x0[notStill], y0[notStill], x1[notStill], y1[notStill])
mtext(sprintf(
"x0=%.2g y0=%.2g u=%.2gm/s v=%.2gm/s speed=%.2gm/s knot=%.2g knot5=%g barbType=%d",
x0, y0, uOrig, vOrig, speed, knot, knot5, barbType
), line = 0)
mtext("FIXME: draw barbs", line = 1)
Created on 2024-02-18 with reprex v2.1.0
Oh, I think the diagram above at https://github.com/dankelley/oce/issues/2191#issuecomment-1951198920 answers my questions about wind-speed gaps. I'll take a shot at coding this. It really is like making change, drawing in the reverse order as you do the sums in your head.
Here's a sketch of the 'making change' type algorithm.
feathers <- function(u) {
pos <- 0
# Handle triangle feathers
while (u > 50) {
cat("draw 50 at", pos, "\n")
u <- u - 50
pos <- pos + 1
}
# Handle Long feathers
while (u > 10) {
cat("draw 10 at", pos, "\n")
u <- u - 10
pos <- pos + 1
}
# Handle short feathers
while (u > 5) {
cat("draw 5 at", pos, "\n")
u <- u - 5
pos <- pos + 1
}
}
Hi all!
First, because this was opened as a "discussion", I think it's worth discussing whether it makes sense to add more meteorology functionality to oce (this is a discussion that @dankelley and I have more and more with new feature requests, due to the burdens of maintaining what is already a very large package, and that the build requirements for users are somewhat challenging). I can see a number of different points, for and against:
-
oceis a package designed for oceanographic data analysis, so functionality that drifts too far from that might be better in a dedicated package - However, oceanography and meteorology are intrinsically linked so it could be argued that having the functionality together makes sense
- A change like the one being discussed (and that @dankelley has already started working on, lol) might be as simple as just adding a new "arrow type" to the existing
drawDirectionField()andmapDirectionField()(so, fairly minimal change) - As mentioned already,
oceis very big (bloated?) and already difficult to maintain. Adding new features might be better added as separate packages that depend onoce(e.g. as theargoFloatspackage does).
There are probably more points that could be made, but that's a start.
I think I'm mostly in the camp of number 3, and kind of like the idea of being able to have wind-barb maps. But I would caution about adding too many other meteorologically-specific functions beyond that (if that is what @mdupilka is proposing).
Oh, and I don't know if this meets any of your needs @mdupilka, but do you know about the metR package?
I'll add another question to Clark's good list:
- How many oceanographers (oce users and others) understand the barbs on a wind diagram? I didn't, until this morning. And it took me a while to figure it out, from a fair bit of web searching.
Re clark's item 3, it is not actually that minimal. There is that "making change" computation, for one thing. And the there is the matter of doing the trig to add the feathers. The existing drawDirectionField() relies on built-in arrow types, which are done in C, so the trig is (a) much faster and (b) well-tested against problem cases.
I've enjoyed playing with this, but thinking of it during my run, it doesn't feel like a good fit for oce. I'd vote -0.8 on "add this feature". Perhaps @richardsc and @clayton33 can vote on the -1 to 1 scale we use for making decisions. No rush -- I think tomorrow is a holiday here anyway.
Thanks for the discussion. I can only imagine the work required for maintain the OCE package. Thanks a lot for that. I do use the mapDirectionField functions to draw wind vectors. It works very well. I do not have any reference for the wind barbs, The https://www.weather.gov/hfo/windbarbinfo is a standard format for wind barbs. Meteorology pretty much always uses knots for speed, which is a standard for aviation. But some software does give options for km/h display. I can understand the extra work to incorporate meteorology displays would make the package even more "bloated". Tomorrow is a holiday in parts of Canada as well.
I don't know if this may be of any interest, back many years ago I wrote a C++ program to display wind winds. Here is the code for the wind barbs
const double DegToRad = 3.141593/180.;
const int BarbLength = 30;
const double BarbSpeedAngle = 60.;
const int BarbSpeedLength = 10;
const int BarbInt = 3;
double WindDirection = 90. - WndDir;
TPoint TriangleVert[3];
dc->SelectObject(*BlackPen);
dc->SelectObject(*BlackBrush);
//Draw wind barb from center of station
int Num;
if(DrawAll) //Draw all stations and ignore the passed station number
Num = NumStations;
else //Just draw for the passed station number
Num = StnNum + 1;
for(int i = StnNum; i < Num; i++)
{
if(Stations[i].y > yFrontPosition && Stations[i].spd > 0)
{
//Set the origin
int x0 = Stations[i].x;
int y0 = Stations[i].y;
//Set the end point
int x1 = x0 + BarbLength * cos(WindDirection*DegToRad);
int y1 = y0 + BarbLength * sin(WindDirection*DegToRad);
dc->MoveTo(x0, maxY - y0);
dc->LineTo(x1, maxY - y1);
//Draw the wind speed 50 knot triangles
//Number of 50 knots
int Num50Tri = int((Stations[i].spd + 2.5)/50.);
int Space = 0;
//Draw the 50 knot triangles
for(int j = 0; j < Num50Tri; j++)
{
TriangleVert[0].x = x0 + (BarbLength - j*BarbSpeedLength) * cos(WindDirection*DegToRad);
TriangleVert[0].y = y0 + (BarbLength - j*BarbSpeedLength) * sin(WindDirection*DegToRad);
TriangleVert[1].x = x0 + (BarbLength - BarbSpeedLength - j*BarbSpeedLength) * cos(WindDirection*DegToRad);
TriangleVert[1].y = y0 + (BarbLength - BarbSpeedLength - j*BarbSpeedLength) * sin(WindDirection*DegToRad);
TriangleVert[2].x = TriangleVert[1].x + BarbSpeedLength * cos((WindDirection - BarbSpeedAngle)*DegToRad);
TriangleVert[2].y = TriangleVert[1].y + BarbSpeedLength * sin((WindDirection - BarbSpeedAngle)*DegToRad);
TriangleVert[0].y = maxY - TriangleVert[0].y;
TriangleVert[1].y = maxY - TriangleVert[1].y;
TriangleVert[2].y = maxY - TriangleVert[2].y;
dc->Polygon(TriangleVert, 3);
//Set position for next 10 knot barbs
Space = (j+1)*BarbSpeedLength + BarbInt;
}
//Draw the 10 knots wind speed barbs
//Number of 10 knot barbs
int Num10Barbs = int(fmod((Stations[i].spd + 2.5), 50.)/10.);
//Draw the 10 knot barbs
for(j = 0; j < Num10Barbs; j++)
{
int x2 = x0 + (BarbLength - (j*BarbInt + Space)) * cos(WindDirection*DegToRad);
int y2 = y0 + (BarbLength - (j*BarbInt + Space)) * sin(WindDirection*DegToRad);
int x3 = x2 + BarbSpeedLength * cos((WindDirection - BarbSpeedAngle)*DegToRad);
int y3 = y2 + BarbSpeedLength * sin((WindDirection - BarbSpeedAngle)*DegToRad);
dc->MoveTo(x2, maxY - y2);
dc->LineTo(x3, maxY - y3);
//Set space for next 5 knot barbs
if(j == Num10Barbs - 1) Space += (j+1)*BarbInt;
}
//Draw the 5 knot barb if there is one
if((fmod(Stations[i].spd, 10.) + 2.5) >= 5. && (fmod(Stations[i].spd, 10.) + 2.5) < 10.)
{
int x2 = x0 + (BarbLength - Space) * cos(WindDirection*DegToRad);
int y2 = y0 + (BarbLength - Space) * sin(WindDirection*DegToRad);
int x3 = x2 + (BarbSpeedLength+2)/2 * cos((WindDirection - BarbSpeedAngle)*DegToRad);
int y3 = y2 + (BarbSpeedLength+2)/2 * sin((WindDirection - BarbSpeedAngle)*DegToRad);
dc->MoveTo(x2, maxY - y2);
dc->LineTo(x3, maxY - y3);
}
}
}
I made a trial code that does something crude. Rather than get into drawing triangles and lines, I am showing a solid square for 50kt, a big circle for 10kt and a small circle for 5kt. This is just to test my "making change" thought process.
The image shows faked data, with flow going around a circle and increasing with distance from the centre. Frankly, this kind of diagram does not really evoke much to my eye. (The barbs do not, either.). I prefer the convention of longer arrows for stronger fields, as we already have in oce.
A trick with R is that loops etc are slow. The sine and cosine values can be precomputed because the angles are on a 10 degree grid, so that will speed things up. Redrawing the view below takes 0.13s user time but 0.37s elapsed time on my machine. I don't think it's entirely feasible to draw fields with a whole lot more points than this (because over overlap) so it seems that maybe a code in R would not be too slow even if not well vectorized. Say 2X the grid density for 4X the processing time, maybe throw in 2X to draw triangles and such, but that only gets us to about 8X, or on my machine still under 4 seconds.
Still, the arguments against adding this (even cleaned up for barbs) to oce are quite strong, I think. The audience would be small, to begin with. Add to that the fact that any code that is written has to be maintained. And pretty soon we would get requests that it work on projected maps. Then options would be needed. I would want m/s instead of knots, for example.
@mdupilka I don't think you'd find it too hard to translate your C code into R, if you really wanted to. The key functions you'd want are plot(), points(), segments() and polygon(). Oh, and the aspect ratio of the plot comes into it as well, of course. Lots of fiddly details. But it could be done. Some vectorization work might be needed to get speed, and it would make sense to precompute the sine and cosine values, which is possible because I think the official schemes quantize to 10degree increments.
Perhaps @richardsc and @clayton33 will have comments on Tuesday, after the long weekend.
Thanks for working on that. Meteorologists would be the ones using the barbs as that is the standard display for wind fields. Maybe that audience is not too large. I will look into adapting my code. I need it to run as fast as possible as the charts are used operationally in our forecasting.
I decided to look at the barbs again. There are some adjustable parameters to the aesthetics (angle of barbs, length of barbs, thickness of lines) that all ought to be parameters of a function to do this work. Below is a snapshot of something I just tried. I think it looks a lot like what I see at one of the examples shown above (reproduced here too).
Trial version
An official diagram
It looks very good. Just like what is on weather charts. Just a note for speeds of 3-7 units, there is just the short feather positioned one space from the end. Like the second barb you have, only without the 10 unit feather. Just as a side note about me, I have been a professional and research meteorologist/atmospheric scientist for longer than I care to remember.
Thanks @mdupilka. I'll change that. What I'm doing is to round to nearest 5 knots, so that category would be from 2.5 to 7.5 knots.
Also, I've seen that Canadians (or maybe Canadian modellers) have a similar-looking set of symbols that represent km/hour. I think that scheme is useful for the public since all cars here show speeds in that unit, so you can imagine rolling down your window and feeling the wind, to get an idea what numbers mean. I guess Americans just reckon a knot is like a mile per hour ... and 10% is beyond our ability to intuit, so no problems.
@mdupilka how does the top row look now? The image is not very clear, because I was trying to reproduce the "official" image. (I wouldn't use such a thick line width, I don't think.)
Trial version
An official diagram
I have to say, despite the work involved in getting this implemented, I'm coming around that it makes sense to have wind barbs included in oce. Most large-scale oceanography papers, which need to be concerned with atmospheric forcing and buoyancy fluxes etc, will include maps of barometric pressure and wind, and this makes it easy to satisfy both the oceanographic (arrows) and meteorologic (barbs) people.
As an aside, a quick google suggests that you've been involved in the CMOS community in the past. I've mostly only been involved since about 2016, but I wonder if you ever crossed paths with my father, Bill Richards, who was a meteorologist with the MSC in NS and NB?
@richardsc I agree it will be good to have in oce. I'll likely add it on the weekend, I just need to get the 50-knot triangles working, then figure a scheme to do things in km/hour for Canadians (maybe others too) and then, if that all seems ok aesthetically, I can do the next thing, which is save all the plot adornments into a vector that I can use to paint all the barbs at once, for speed. So that's a few little things, but during the week when people tend to be working, I do hope to get the aesthetics right.
The top row looks very good. I would not worry about any units. Just make the feather values unitless and let the user decide what units they wish to use. The wind barbs history is from the aviation world where all speeds are in knots throughout the world, even in Canada. But, for my purposes, I actually use km/h as the folks we deal with in our forecasting do not really understand knots.
I was a chair in CMOS for several years. I did not cross paths with your father. I will ask some other met guys I know if they have. A bit more on my use of OCE, the company I contract with does fire weather forecasting for the Northwest Territories during the summer months. I produce a number of charts related to fire weather prediction and use OCE to map the fields. You may have heard of the fire season the past year in Canada and the NWT, it was insane. I worked nearly every day and I am supposed to be mostly retired. Your OCE charts play an integral part in our forecasts and briefing.
@mdupilka how does this look. It's speeds 0, 5kn, 10kn, etc. in reading order from top left.
I see some variability in how the triangles get draw, even in the images preceding this in this thread. I kind of like this look. The ones that have asymmetric triangles look kind of weird to me. Thanks for your help and patience. Dan.
(Oh, and thanks a million for those kind words about oce. It's always nice to make something that is useful in fields beyond those initially envisioned!)
That looks really good. The symmetric triangles are the standard. A light wind, less than the 5 barb could just be indicated as a line with no barb.
Thanks. If wind is zero, which I imagine never happens, is just a dot shown? Or nothing at all?
Thanks, @mdupilka. I think I have the aesthetics done, as in the snapshot below. The numbers above the locations are just for demonstation purposes -- so you can see that a speed of 2 knots becomes a line, not a circle.). I elected to draw a little circle because the large circle I see at https://github.com/dankelley/oce/issues/2191#issuecomment-1951206825 just seems, well, too large.
Q: Are you okay with the small circle? Obviously, it's trivial for me to change it.
@mdupilka if I put this into oce, my plan would be to do drawDirectionField() first, because for that we don't get into issues with map projections that crop up with mapDirectionField(). For example, consider views that show the whole world, in which case an arrow pointing say east at the east side of the rendered map should go to the edge of the world on that side, and then re-enter at the other edge of the world, in a different place. There are other issues, too.
Would drawDirectionField() be useful to you, or are you using map projections and thus need mapDirectionField()?