[DISCUSSION] Spurious pcolormesh wrapping is fixed

LiamBindle commented 4 years ago

Hi everyone,

I mentioned this on Slack last week, but plotting GCHP full global data with pcolormesh() can be difficult because you'll get horizontal streaks for grid-boxes that cross the antimeridian. Below are examples.

TLDR: Use cartopy version 0.19 or greater if you want to plot GCHP data for the entire globe.

Plotting CS data

Set up the figure

ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_global()
ax.coastlines()

Plot a face that doesn't cross the antimeridian (looks good)

ds = xr.open_dataset('GCHP.SpeciesConc.nc4')
plt.pcolormesh(
    ds.lons.isel(nf=4).values, 
    ds.lats.isel(nf=4).values, 
    ds.SpeciesConc_NO2.isel(nf=4, lev=0, time=0).values,  
    vmax=8e-9
)

Now, plot a face that does cross the antimeridian (results in horizontal streaking)

plt.pcolormesh(
    ds.lons.isel(nf=3).values, 
    ds.lats.isel(nf=3).values, 
    ds.SpeciesConc_NO2.isel(nf=3, lev=0, time=0).values,  
    vmax=8e-9
)

Plotting SG data

This is illustrated a bit better with stretched-grids.

Again, set up the figure

ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_global()
ax.coastlines()

Plot a face that doesn't cross the antimeridian

ds = xr.open_dataset('GCHP.SpeciesConc.nc4')
plt.pcolormesh(
    ds.lons.isel(nf=0).values, 
    ds.lats.isel(nf=0).values, 
    ds.SpeciesConc_NO2.isel(nf=0, lev=0, time=0).values, 
    vmax=8e-9
)

Next, plot a face that does cross the AM

plt.pcolormesh(
    ds.lons.isel(nf=1).values, 
    ds.lats.isel(nf=1).values, 
    ds.SpeciesConc_NO2.isel(nf=1, lev=0, time=0).values, 
    vmax=8e-9
)

The issue and workarounds

The PlateCarree projection is a 2D space unaware of wrapping at the antimeridian. Work arounds include:

Don't plot data that's close to the antimeridian (e.g., GCPy doesn't plot grid-boxes that are within 2 deg of the AM; see here)
Manually draw polygons in a gnominic projection for grid-boxes that cross the AM (this is what I've prefered doing because it properly wraps the data around the AM)

Related issues

A recent fix

Fixed in https://github.com/SciTools/cartopy/pull/1622 (thanks @htonchia and @greglucas). This merge was on August 19, 2020. IIUC this fix will be released in cartopy version 0.19.

LiamBindle commented 4 years ago

Just following up on my prefered workaround. Here's a minimal example.

I forgot to mention that with this method you have to plot each grid-box that crosses the AM as a polygon and manually specify its color. Since this is a minimal example, I omitted logic to only do this for grid-boxes that cross the AM, and instead I do it for all grid-boxes on nf=1. This can be sped up by using masked arrays with pcolormesh() and then only drawing polygons for grid-boxes that cross the AM.

Note: This method also requires grid-box corners.

import shapely.geometry
import pyproj

# Get grid corners in epsg:4326
lons = grid.xe.isel(nf=1).values   # Note: grid.xe and grid.ye are the grid's corners
lats = grid.ye.isel(nf=1).values

# Project corners to gnomonic space centered on the face nf=1
gno = ccrs.Gnomonic(20, -150)
xe, ye = pyproj.transform('epsg:4326', gno.proj4_init, lons, lats, always_xy=True)

# We have to manually assign colors to polygons, tf. create norm and cmap
cmap = plt.get_cmap('viridis')
norm = plt.Normalize(vmin=0, vmax=8e-9)

for j in range(ds.dims['Ydim']):
    for i in range(ds.dims['Xdim']):
        fill=cmap(norm(ds.SpeciesConc_NO2.isel(lev=0, nf=1, Ydim=j, Xdim=i).values[0])) # color of polygon
        polygon = shapely.geometry.Polygon([(xe[j,i], ye[j,i]), (xe[j+1,i], ye[j+1,i]), (xe[j+1,i+1], ye[j+1,i+1]), (xe[j,i+1], ye[j,i+1])]) # create a polygon
        ax.add_geometries( # add the polygon to the plot
            [polygon],
            gno, edgecolor=fill, facecolor=fill
        )

A full example of this method can be found here

This can make plots that look like this

LiamBindle commented 4 years ago

I just tried cartopy version 0.18.1.dev113+g3565f78 (which includes SciTools/cartopy#1622). Here is the result.

ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_global()
ax.coastlines()

plt.pcolormesh(
    ds.lons.isel(nf=0).values, 
    ds.lats.isel(nf=0).values, 
    ds.SpeciesConc_NO2.isel(nf=0, lev=0, time=0).values, 
    vmax=8e-9
)

plt.pcolormesh(
    ds.lons.isel(nf=1).values, 
    ds.lats.isel(nf=1).values, 
    ds.SpeciesConc_NO2.isel(nf=1, lev=0, time=0).values, 
    vmax=8e-9
)

So SciTools/cartopy#1622 does fix the horizontal streaking that we see. Thanks again @htonchia and @greglucas, this is terrific!

WilliamDowns commented 4 years ago

Good to know this will be officially fixed with Cartopy 0.19, which will make it easy to include in GCPy once we can confirm that our plotting works correctly in that new version.

geoschem / GCHP