[DISCUSSION] Spurious pcolormesh wrapping is fixed

[geoschem]

[DISCUSSION] Spurious pcolormesh wrapping is fixed

[LiamBindle]

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!

[LiamBindle]

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

[WilliamDowns]

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.