diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 0005dcb..65e604f 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -13,6 +13,8 @@ Fixed `#29 `__). - fixed plotting of data with 3D bounds (see `#30 `__) +- fixed plotting of curvilinear data with 3D bounds (see + `#31 `__) v1.3.0 ====== diff --git a/psy_maps/plotters.py b/psy_maps/plotters.py index efbcfbb..aeec1ca 100755 --- a/psy_maps/plotters.py +++ b/psy_maps/plotters.py @@ -1679,6 +1679,10 @@ def _polycolor(self): # However, in order to wrap the boundaries correctly, we have to # identify the corresponding grid cells and then use the standard # matplotlib transform. + if xb.ndim > 2: + xb = xb.reshape((-1, xb.shape[-1])) + yb = yb.reshape((-1, yb.shape[-1])) + arr = arr.reshape(-1) if isinstance(t, wrap_proj_types) and 'lon_0' in proj.proj4_params: # We adopt and copy some code from the methodology of cartopy # _pcolormesh_patched method of the geoaxes. As such, we @@ -1717,9 +1721,6 @@ def _polycolor(self): yb_wrap = yb[mask] xb = xb[~mask] yb = yb[~mask] - if xb.ndim > 2: - xb = xb.reshape((-1, xb.shape[-1])) - yb = yb.reshape((-1, yb.shape[-1])) self.logger.debug('Making plot with %i cells', arr.size) transformed = proj.transform_points( t, xb.ravel(), yb.ravel())[..., :2].reshape(xb.shape + (2,)) @@ -1791,6 +1792,11 @@ def update(self, value): proj = self.ax.projection # HACK: We remove the cells at the boundary of the map projection xb_wrap = None + + if xb.ndim > 2: + xb = xb.reshape((-1, xb.shape[-1])) + yb = yb.reshape((-1, yb.shape[-1])) + if isinstance(t, wrap_proj_types) and 'lon_0' in proj.proj4_params: # See the :meth:`MapPlot2D._polycolor` method for a documentation # of the steps @@ -1822,20 +1828,48 @@ def update(self, value): transformed = proj.transform_points(t, xb.ravel(), yb.ravel()) xb = transformed[..., 0].reshape(orig_shape) yb = transformed[..., 1].reshape(orig_shape) - if xb.ndim > 2: - xb = xb.reshape((-1, xb.shape[-1])) - yb = yb.reshape((-1, yb.shape[-1])) # We insert nan values in the flattened edges arrays rather than # plotting the grid cells separately as it considerably speeds-up code # execution. n = len(xb) xb = np.c_[xb, xb[:, :1], [[np.nan]] * n].ravel() yb = np.c_[yb, yb[:, :1], [[np.nan]] * n].ravel() + if isinstance(value, dict): self._artists = self.ax.plot(xb, yb, **value) else: self._artists = self.ax.plot(xb, yb, value) + if xb_wrap is not None: + # first we drop all grid cells that have any NaN in their bounds + # as we do not know, how to handle these + valid = (~np.isnan(xb_wrap).any(-1)) & (~np.isnan(yb_wrap).any(-1)) + xb_wrap = xb_wrap[valid] + yb_wrap = yb_wrap[valid] + + if isinstance(t, ccrs.PlateCarree): + + # identify the grid cells at the boundary + xdiff2min = xb_wrap - xb_wrap.min(axis=-1, keepdims=True) + cross_world_mask = np.any(np.abs(xdiff2min) > 180, -1) + if cross_world_mask.any(): + cross_world_x = xb_wrap[cross_world_mask] + cross_world_y = yb_wrap[cross_world_mask] + cross_world_diff = xdiff2min[cross_world_mask] + xdiff2max = cross_world_x - cross_world_x.max( + axis=-1, keepdims=True) + + leftx = cross_world_x.copy() + leftx[cross_world_diff > 180] -= 360 + + rightx = cross_world_x.copy() + rightx[xdiff2max < -180] += 360 + + xb_wrap = np.r_[xb_wrap[~cross_world_mask], leftx, rightx] + yb_wrap = np.r_[ + yb_wrap[~cross_world_mask], cross_world_y, cross_world_y + ] + n = len(xb_wrap) xb_wrap = np.c_[xb_wrap, xb_wrap[:, :1], [[np.nan]] * n].ravel() yb_wrap = np.c_[yb_wrap, yb_wrap[:, :1], [[np.nan]] * n].ravel() diff --git a/tests/curvilinear-with-bounds.nc b/tests/curvilinear-with-bounds.nc new file mode 100644 index 0000000..9aa68ae Binary files /dev/null and b/tests/curvilinear-with-bounds.nc differ diff --git a/tests/test_plotters.py b/tests/test_plotters.py index 117deb2..f34e07d 100755 --- a/tests/test_plotters.py +++ b/tests/test_plotters.py @@ -1774,5 +1774,30 @@ def test_datagrid_3D_bounds(): assert abs(ymax - ymin - 52) < 2 +def test_plot_curvilinear_datagrid(tmpdir): + """Test if the there is no datagrid plotted over land + + This implicitly checks, if grid cells at the boundary are warped correctly. + The file ``'curvilinear-with-bounds.nc'`` contains a variable on a + curvilinear grid that is only defined over the ocean (derived from MPI-OM). + Within this test, we focus on a region over land far away from + the ocean (Czech Republic) where there are no grid cells. If the datagrid + is plotted correctly, it should be all white. + """ + from matplotlib.testing.compare import compare_images + fname = os.path.join(bt.test_dir, 'curvilinear-with-bounds.nc') + # make a white plot without datagrid + kws = dict(plot=None, xgrid=False, ygrid=False, map_extent='Czech Republic') + with psy.plot.mapplot(fname, **kws) as sp: + sp.export(str(tmpdir / "ref.png")) # this is all white + # now draw the datagrid, it should still be empty (as the input file only + # defines the data over the ocean) + with psy.plot.mapplot(fname, datagrid='k-', **kws) as sp: + sp.export(str(tmpdir / "test.png")) # this should be all white, too + results = compare_images( + str(tmpdir / "ref.png"), str(tmpdir / "test.png"), tol=1) + assert results is None, results + + if __name__ == '__main__': bt.RefTestProgram()