Modifications to enable accurate plotting of cubed-sphere data which …

…crosses 180E

An upcoming update to Cartopy (expected in 0.19.0) enables pcolormesh to correctly
display polygons crossing the antimeridian. This removes the need for thresholding
at that location, and means that ALL cubed sphere data is now shown on plots. Users
can recover the old functionality by passing the argument `cs_threshold=5` when
calling `plot_layer` (or the lower-level `plot_cs` function). A side effect of this
fix is that a wrapped set of polgyons is sometimes generated when plotting one of
the cubed-sphere faces. `plottools.set_clim` has been modified to account for this.

gcgridobj/plottools.py

@@ -204,7 +204,7 @@ def tick_label_gen(ticks):
 
     return im, cb
 
-def plot_layer(layer_data,hrz_grid=None,ax=None,crs_data=None,crs_plot=None,show_colorbar=True,coastlines=True):
+def plot_layer(layer_data,hrz_grid=None,ax=None,crs_data=None,crs_plot=None,show_colorbar=True,coastlines=True,**kwargs):
 
     if crs_data is None:
        crs_data = crs_data_standard
@@ -224,14 +224,14 @@ def plot_layer(layer_data,hrz_grid=None,ax=None,crs_data=None,crs_plot=None,show
     elif len(ld_shape) == 3:
        # Could this be a valid CS dataset?
        assert ld_shape[0] == 6 and (ld_shape[1] == ld_shape[2]), 'Layer data shape invalid (3D and not CS)'
-       im_obj = plot_cs(layer_data,hrz_grid=hrz_grid,ax=ax,crs_data=crs_data,crs_plot=crs_plot)
+       im_obj = plot_cs(layer_data,hrz_grid=hrz_grid,ax=ax,crs_data=crs_data,crs_plot=crs_plot,**kwargs)
     elif ld_shape[0] == 6*ld_shape[1]:
        # Assume cubed sphere
        is_cs = True
-       im_obj = plot_cs(layer_data,hrz_grid=hrz_grid,ax=ax,crs_data=crs_data,crs_plot=crs_plot)
+       im_obj = plot_cs(layer_data,hrz_grid=hrz_grid,ax=ax,crs_data=crs_data,crs_plot=crs_plot,**kwargs)
     else:
        # Assume lat-lon
-       im_obj = plot_latlon(layer_data,hrz_grid=hrz_grid,ax=ax,crs_data=crs_data,crs_plot=crs_plot)
+       im_obj = plot_latlon(layer_data,hrz_grid=hrz_grid,ax=ax,crs_data=crs_data,crs_plot=crs_plot,**kwargs)
 
     # If cubed-sphere, use the first image
     is_cs = isinstance(im_obj, list)
@@ -273,18 +273,19 @@ def plot_latlon(layer_data,hrz_grid=None,ax=None,crs_data=None,crs_plot=None,sho
 
     return im
 
-def update_cs(layer_data,im_vec,hrz_grid=None,cs_threshold=5):
+def update_cs(layer_data,im_vec,hrz_grid=None,cs_threshold=None):
     # WARNING: layer_data must be [6 x N x N]
-    if hrz_grid is None:
-       # Try to figure out the grid from the layer data
-       #n_cs, is_gmao = guess_cs_grid(layer_data.shape)    
-       #hrz_grid = cubedsphere.csgrid_GMAO(n_cs)
-       hrz_grid = regrid.guess_cs_grid(layer_data.shape)    
-    masked_data = np.ma.masked_where(np.abs(hrz_grid['lon'] - 180.0) < cs_threshold, layer_data)
+    if cs_threshold is not None:
+        if hrz_grid is None:
+            # Try to figure out the grid from the layer data
+            hrz_grid = regrid.guess_cs_grid(layer_data.shape)    
+        masked_data = np.ma.masked_where(np.abs(hrz_grid['lon'] - 180.0) < cs_threshold, layer_data)
+    else:
+        masked_data = layer_data
     for i_face in range(6):
-       im_vec[i_face].set_array(masked_data[i_face,:,:].ravel())
+        im_vec[i_face].set_array(masked_data[i_face,:,:].ravel())
 
-def plot_cs(layer_data,hrz_grid=None,ax=None,crs_data=None,crs_plot=None,show_colorbar=True,cs_threshold=5.0):
+def plot_cs(layer_data,hrz_grid=None,ax=None,crs_data=None,crs_plot=None,show_colorbar=True,cs_threshold=None):
 
     # 2019-12-17: dropped support for non-GMAO grids
     #n_cs, is_gmao = regrid.guess_n_cs(layer_data.shape)    
@@ -300,20 +301,27 @@ def plot_cs(layer_data,hrz_grid=None,ax=None,crs_data=None,crs_plot=None,show_co
        # Try to figure out the grid from the layer data
        #hrz_grid = cubedsphere.csgrid_GMAO(n_cs)
        hrz_grid = regrid.guess_cs_grid(layer_data.shape)
-
-    masked_data = np.ma.masked_where(np.abs(hrz_grid['lon'] - 180.0) < cs_threshold, layer_data)
+
+    # A pending PR for Cartopy is expected to fix an issue where cells crossing the antimeridian
+    # are incorrectly plotted. Until then (and for users with older versions of cartopy), setting
+    # cs_threshold to a non-zero value will at least mask out these cells. A value of 5 seems to
+    # work (somewhat) well, at least for C24 data.
+    if cs_threshold is not None:
+        masked_data = np.ma.masked_where(np.abs(hrz_grid['lon'] - 180.0) < cs_threshold, layer_data)
+    else:
+        masked_data = layer_data
 
     im_vec = []
     for i_face in range(6):
-       im = ax.pcolormesh(hrz_grid['lon_b'][i_face,:,:],hrz_grid['lat_b'][i_face,:,:],masked_data[i_face,:,:],transform=crs_data)
+       lon_b = np.mod(hrz_grid['lon_b'][i_face,:,:],360.0)
+       im = ax.pcolormesh(lon_b,hrz_grid['lat_b'][i_face,:,:],masked_data[i_face,:,:],transform=crs_data)
        im_vec.append(im)
 
     c_lim = [np.min(layer_data),np.max(layer_data)]
     if (c_lim[0] == c_lim[1]):
        c_lim = [c_lim[0] - 0.5,c_lim[1] + 0.5]
 
-    for im in im_vec:
-       im.set_clim(c_lim)
+    set_clim(im_vec,c_lim)
 
     return im_vec
 
@@ -389,6 +397,11 @@ def set_clim(im_obj,c_lim=None,cmap=None):
           # Assume max
           c_lim = [-c_lim,c_lim]
        im_obj.set_clim(c_lim)
+       # Since Cartopy 0.19.0, there is a hidden set of polygons for GeoQuadMesh objects
+       # if they cross the antimeridian. Their color limits need to be fixed separately
+       wrapped_obj = getattr(im_obj, "_wrapped_collection_fix", None)
+       if wrapped_obj is not None:
+          wrapped_obj.set_clim(c_lim)
        if cmap is not None:
           im_obj.set_cmap(cmap)
     return None