MetPy
MetPy copied to clipboard
ImagePlot fails on some [0, 360] data
This example:
from metpy.plots import ImagePlot, MapPanel, PanelContainer
from metpy.units import units
import xarray as xr
ds_z500 = xr.open_dataset('/Users/rmay/Downloads/gfs.t00z.pgrb2.0p25.f024', engine='cfgrib', filter_by_keys={'typeOfLevel': 'isobaricInhPa', 'level': 500})
img = ImagePlot()
img.data = ds_z500
img.field = 't'
img.colormap = 'plasma'
panel = MapPanel()
panel.area = 'us'
panel.title = 'GFS 500mb Temp Forecast Example'
panel.plots = [img]
pc = PanelContainer()
pc.size = (10,8)
pc.panels = [panel]
pc.show()
produces the following with 1.0 and main:
Taken from a Stack Overflow question. This pretty clearly relates to this code that wraps longitude:
https://github.com/Unidata/MetPy/blob/6944e23d13aef4d1c2226b69bd623cb1ab61f09c/src/metpy/plots/declarative.py#L1087-L1102
because without it the same code gives:
The code is pretty clearly incorrect because it results in an x array where x[0]
is 0 and x[-1]
is -0.25, since this full dataset has a longitude range of [0, 360].
Of course it's not as simple as removing that bogus wrapping code, as that breaks a the test_latlon
test:
from datetime import datetime
from metpy.cbook import get_test_data
from metpy.plots import ContourPlot, ImagePlot, MapPanel, PanelContainer
from metpy.units import units
import xarray as xr
data = xr.open_dataset(get_test_data('irma_gfs_example.nc', as_file_obj=False))
img = ImagePlot()
img.data = data
img.field = 'Temperature_isobaric'
img.level = 500 * units.hPa
img.time = datetime(2017, 9, 5, 15, 0, 0)
#img.colorbar = None
contour = ContourPlot()
contour.data = data
contour.field = 'Geopotential_height_isobaric'
contour.level = img.level
contour.time = img.time
panel = MapPanel()
#panel.projection = 'lcc'
panel.area = 'global'
panel.plots = [img, contour]
pc = PanelContainer(size=(8,5))
pc.panel = panel
pc.draw()
produces
so that's lovely.
Everything works fine if I change the extent of the image data (either with extent
manually or in the data) from [250, 315] to [-110, -45 ]. It also behaves more as expected if the range is [179, 315].
I've so far traced this into CartoPy's image warping code. It seems to be an issue where we're using the fact that PlateCarree
can (sort of) handle longitudes in [0, 360], but when you transform points back into PlateCarree
, they're in [-180, 180]. I'll continue digging into CartoPy to try to figure out how to make it not completely broken (seeing as NCEP produces data that are natively [0, 360]).
In the meanwhile, though, we need to unbreak the original example with MetPy, but also keep the test case working. This probably means writing some code to properly "roll" an image from [0, 360] to [-180, 180].
I had to face this issue in the past where I had to wrestle with plotting with longitudes from -180 to 180 degrees. When plotting on a map with this set of longitudes it would warp the lines and wrap around the globe when taken from platecarree. The only solution that I found that prevented these problems was to convert the longitudes to 0 to 360 degrees and it works over the poles too. I tried both ways and looked at the results. It was quite convincing.
Jotting down a thought here:
- Add a method to our xarray accessor(s) to use CartoPy's
add_cyclic_point
(or implement our own for dependency avoidance) - Use that to avoid any need for hacks to deal with this