MetPy icon indicating copy to clipboard operation
MetPy copied to clipboard

ImagePlot fails on some [0, 360] data

Open dopplershift opened this issue 3 years ago • 4 comments

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: image

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: image

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].

dopplershift avatar Jun 16 '21 05:06 dopplershift

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 image

so that's lovely.

dopplershift avatar Jun 16 '21 05:06 dopplershift

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].

dopplershift avatar Jun 16 '21 05:06 dopplershift

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.

eliteuser26 avatar Aug 10 '21 16:08 eliteuser26

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

dopplershift avatar Apr 27 '22 20:04 dopplershift