Unidata
Fronts Plotting Support
Description Of Changes
This adds path effects, mostly used for lines, to plot fronts. Still todo:
- [ ] Fix dpi-related plotting issue
- [ ] Integrate with #1944 or rebase after merging that one
- [ ] Clean up example more
- [ ] Add more docs for methods
- [ ] Add docstring basic plotting example?
Checklist
- [x] Closes #1153
- [x] Tests added
- [ ] Fully documented
A couple of quick observations (could be wrong as well): There could be instances when decode() function does not like CODSUS formatting. All the examples below are from this surface bulletin (https://mesonet.agron.iastate.edu/wx/afos/p.php?pil=CODSUS&e=201604301800).
For example,
LOWS 1002 2220972 1014 4421147 1009 3761203 1002 3411136 1006 3921062 1005
3270946 1004 3351042 1004 4010960 1016 3930803 1020 4310722 1009 6461225 1007
6721632 997 6050376 996 5230464 1010 7710630 1010 7020544 1005 6030464 982
5521520
If I try reading the surface bulletin as it is, I'll get the following error:
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Input In [200], in <cell line: 1>()
----> 1 decode(ff)
Input In [184], in decode(fobj)
33 ret['VALID'] = info[0]
34 elif field in ('HIGHS', 'LOWS'):
---> 35 ret.setdefault(field, []).extend(parse_points(info))
36 elif field in ('WARM', 'COLD', 'OCFNT', 'STNRY', 'TROF'):
37 dest = ret.setdefault(field, [])
Input In [184], in parse_points(strs)
11 for mag_str, loc_str in zip(strs[:-1:2], strs[1::2]):
12 mag = int(mag_str)
---> 13 yield *str_to_lonlat(loc_str), mag
Input In [184], in str_to_lonlat(s)
7 def str_to_lonlat(s):
----> 8 return -int(s[3:]) / 10, int(s[:3]) / 10
ValueError: invalid literal for int() with base 10: ''
Notice that all three lines end with the pressure value. After reformatting it to include the location along with the pressure value in the same line, we get the correct output e.g.,
LOWS 1002 2220972 1014 4421147 1009 3761203 1002 3411136 1006 3921062 1005 3270946
1004 3351042 1004 4010960 1016 3930803 1020 4310722 1009 6461225 1007 6721632
997 6050376 996 5230464 1010 7710630 1010 7020544 1005 6030464 982 5521520
decode output of formatted data
'LOWS': [(-97.2, 22.2, 1002),
(-114.7, 44.2, 1014),
(-120.3, 37.6, 1009),
(-113.6, 34.1, 1002),
(-106.2, 39.2, 1006),
(-94.6, 32.7, 1005),
(-104.2, 33.5, 1004),
(-96.0, 40.1, 1004),
(-80.3, 39.3, 1016),
(-72.2, 43.1, 1020),
(-122.5, 64.6, 1009),
(-163.2, 67.2, 1007),
(-37.6, 60.5, 997),
(-46.4, 52.3, 996),
(-63.0, 77.1, 1010),
(-54.4, 70.2, 1010),
(-46.4, 60.3, 1005),
(-152.0, 55.2, 982)],
Another issue could be missing or abbreviated pressure values in surface bulletin e.g.,
HIGHS 1012 3601014 1024 3970743 4330679 1031 5260983 1024 4690627 1034 6130890
1028 5241162 1030 5281237 1029 7040890 1034 7911217 1023 4420601 1032 4511291
1031 4091323
Notice the missing pressure value after 3970743 (or maybe they meant that pressure was 1024 at 3970743 and 43306679 but decode had issues in parsing that info. Here's what it returns:
'HIGHS': [(-101.4, 36.0, 1012),
(-74.3, 39.7, 1024),
(-0.1, 10.3, 4330679),
(-0.4, 10.2, 5260983),
(-0.4, 10.3, 4690627),
(-116.2, 52.4, 1028),
(-123.7, 52.8, 1030),
(-89.0, 70.4, 1029),
(-121.7, 79.1, 1034),
(-60.1, 44.2, 1023),
(-129.1, 45.1, 1032),
(-132.3, 40.9, 1031)],
Therefore, CODSUS string may need further scrutiny/checks before plotting fronts.
Also, I think decode has issues parsing weak fronts, whenever present. For example, in this surface bulletin (https://mesonet.agron.iastate.edu/wx/afos/p.php?pil=CODSUS&e=202205060437):
COLD WK 3595 3494 3395 3295 3296
OCFNT WK 5554 5555 5657 5757 5856 5954
OCFNT WK 5954 6051 6049 6046
WARM WK 6046 5944 5841
STNRY WK 5251 5252 5153
will throw an error:
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Input In [981], in <cell line: 1>()
----> 1 decode(data)
Input In [973], in decode(fobj)
36 elif field in ('WARM', 'COLD', 'OCFNT', 'STNRY', 'TROF'):
37 dest = ret.setdefault(field, [])
---> 38 vals = parse_path(info)
39 if continuing:
40 dest[-1].extend(vals)
Input In [973], in parse_path(strs)
15 def parse_path(strs):
---> 16 return [str_to_lonlat(s) for s in strs]
Input In [973], in <listcomp>(.0)
15 def parse_path(strs):
---> 16 return [str_to_lonlat(s) for s in strs]
Input In [973], in str_to_lonlat(s)
7 def str_to_lonlat(s):
----> 8 return -int(s[3:]) / 10, int(s[:3]) / 10
ValueError: invalid literal for int() with base 10: ''
Oooh, these will be good samples to include in tests in #1944. Thanks!
I've been plotting fronts and think that it actually works pretty well in most cases. There are some things I noticed, though. I'll add comments for each that I find.
For reference, I've been using Python 3.10.4 and matplotlib 3.5.1.
The first issue I noticed was with stationary fronts. See the attached image below where there appears to be some filling going on in the curved parts of the line. This does not happen with the other types of fronts.
This snippet should be able to reproduce this:
import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate
from metpy.plots.patheffects import StationaryFront
n = 100
front_lat = [31.08, 31.76, 32.18, 32.01, 31.7, 32.44]
front_lon = [-88.08, -89.23, -90.83, -92.08, -93.27, -95.55]
index = np.arange(len(front_lon))
ii = np.linspace(0, index.max(), n)
xcurve = interpolate.interp1d(index, front_lon, 'cubic')(ii)
ycurve = interpolate.interp1d(index, front_lat, 'cubic')(ii)
plt.plot(xcurve, ycurve, path_effects=[StationaryFront(spacing=6, size=6)],
linewidth=2)
plt.show()
From what I can tell, the problem probably occurs when creating new graphic contexts to plot each segment of the front. I can take the same coordinates and spacing, etc., and not have the issue when plotting a warm or cold front. Somehow matplotlib decides to fill areas when it should not.
While fronts are not too often highly curved, there are cases, such as with highly occluded cyclones, where it can occur. This can lead to situations where the fronts barbs/scallops can appear off of the line. See the example below.
Here is the snippet that should reproduce it:
import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate
from metpy.plots.patheffects import OccludedFront
n = 100
front_lat = [44.57, 44.15, 44.18, 44.57, 45.13, 45.43, 45.14,
44.37, 43.24, 42.19, 41.43, 41.04]
front_lon = [-99.12, -99.13, -99.78, -99.94, -99.26, -98.17,
-97.06, -95.76, -94.23, -93.18, -92.68, -92.4]
index = np.arange(len(front_lon))
ii = np.linspace(0, index.max(), n)
xcurve = interpolate.interp1d(index, front_lon, 'cubic')(ii)
ycurve = interpolate.interp1d(index, front_lat, 'cubic')(ii)
plt.plot(xcurve, ycurve, path_effects=[OccludedFront(spacing=6, size=6)],
linewidth=0.5)
plt.show()
I toned down the linewidth to make it a little more obvious. That being said, most fronts, when smoothed, look pretty good. Setting a slightly thicker line width does help in a lot of these cases, too. I think what NMAP ends up doing is taking the curved path into account and just filling between the line and the barb/scallop.
To my eye, I think the symbols on the stationary front are not exactly centered within the segments. It is not a significant amount, but it looks like there is more path on the left than the right. See the attached image for an example.
This snippet should reproduce it:
import matplotlib.pyplot as plt
from metpy.plots.patheffects import StationaryFront
x = [0, 1]
y = [1, 1]
plt.plot(x, y, path_effects=[StationaryFront(spacing=5, size=15)],
linewidth=2)
plt.show()
Choosing larger linewidth values begins to reduce the amount of the scallops/pips that are visible and rounds the edges of pips. The attached image shows an example of this.
It seems unlikely that you would ever want to increase the width to something this large, but it may be something that has at least some impact on more realistic examples. In my testing, you could mitigate the issue by creating a new graphics context with a different linewidth when drawing the symbols. Of course, if you still had a relatively large linewidth, the line portion could still cover up more of the symbol than you'd like. Perhaps translating the symbol a bit based on the linewidth would be the best solution? Lastly, a larger linewidth accentuates how off-center the symbols appears.
Figure reproduced with the following snippet:
import matplotlib.pyplot as plt
from metpy.plots.patheffects import StationaryFront
x = [0, 1]
y = [1, 1]
plt.plot(x, y, path_effects=[StationaryFront(spacing=5, size=15)],
linewidth=15)
plt.show()
Thanks for examining all these corner cases, these are definitely going to help straighten out this implementation.
Ok, thanks again to @nawendt for finding some corner cases. I've been able to fix some things:
- Need to scale our size from points to pixels everywhere, this fixed the dpi issue as well as some of the spacing on stationary front
- Fixed the facecolor issue @nawendt identified
- Fixed the stationary front issues with linewidth by using that as an offset normal to the path.
- Adjusted drawing style to use 'miter' join style and 'butt' capstyle, which helped with some of the other issues like rounding and some visual "offsets"
I've incorporated some of those as tests now and you can see in the baseline images what things look like. The only outstanding issue is the aspect where the markers pull away from curved paths. Currently, this is implemented by walking along the path and stamping fixed markers that are rotated to be tangent to the path. To fix that issue, we'd have to change to do more direct generation of the markers or event a full Path of all the markers. Not a bad idea, but it would fundamentally change the approach and delay this further--therefore I'm electing to delay that so we can get something useful out there.
I'm also electing to just wait on the example on getting #1944 in, which is next on my list to wrap this up. This is useful to have as-is.
In working with the initial implementation of this, I also made some variations of the front types that represented strengthening and weakening boundaries. These front symbols are used in NMAP and described here. I know we will use them from time to time at SPC. WPC probably uses them more often on their analysis charts. Would you like that part of this PR or part of separate PR down the line? I'll have to update them to use the fixes that you added recently.
@nawendt How about as a separate PR since those are all your additions? I think what we have does everything it needs to meet what's documented for the WPC bulletins, though I'm happy to be shown additional products or samples. The GEMPAK docs are helpful, but I'm not sure what a practical sample looks like. I also want to get this across the line so long as it's considered sound technically.
OK. I'll plan on opening a PR once you finalize and merge this implementation.