diff --git a/gpm_api/utils/area.py b/gpm_api/utils/area.py index 3c768651..3fe5ac9c 100644 --- a/gpm_api/utils/area.py +++ b/gpm_api/utils/area.py @@ -9,6 +9,10 @@ ### SwathDefinition --> geographic coordinates --> computation in ECEF (xyz) ### AreaDefinition --> proj coordinates --> computation in projection space ! +# quadrilateral_corners / quadmesh +# get_lonlat_corners +# get_projection_corners + def _infer_interval_breaks(coord, axis=0): """Infer the outer and inner midpoints of 2D coordinate arrays. diff --git a/gpm_api/visualization/orbit.py b/gpm_api/visualization/orbit.py index cebe95a4..1384f76e 100644 --- a/gpm_api/visualization/orbit.py +++ b/gpm_api/visualization/orbit.py @@ -191,7 +191,7 @@ def wrapper(*args, **kwargs): # - Get slices with contiguous scans and valid geolocation list_slices = get_slices_regular(da) if len(list_slices) == 0: - return ValueError("No regular scans available. Impossible to plot.") + raise ValueError("No regular scans available. Impossible to plot.") else: list_slices = [slice(0, None)] @@ -342,7 +342,7 @@ def plot_orbit_mesh( # - Define plot_kwargs to display only the mesh plot_kwargs["facecolor"] = "none" plot_kwargs["alpha"] = 1 - plot_kwargs["edgecolors"] = (edgecolors,) + plot_kwargs["edgecolors"] = (edgecolors,) # Introduce bugs in Cartopy ! plot_kwargs["linewidth"] = (linewidth,) plot_kwargs["antialiased"] = True @@ -376,12 +376,6 @@ def plot_orbit_image( check_contiguous_scans(da) _preprocess_figure_args(ax=ax, fig_kwargs=fig_kwargs) - # - Define default x and y - # if x is None: - # x = "along_track" - # if y is None: - # y = "cross_track" - # - Initialize figure if ax is None: fig, ax = plt.subplots(**fig_kwargs) diff --git a/gpm_api/visualization/plot.py b/gpm_api/visualization/plot.py index 0403147a..3a114042 100644 --- a/gpm_api/visualization/plot.py +++ b/gpm_api/visualization/plot.py @@ -59,6 +59,27 @@ def get_antimeridian_mask(lons, buffer=True): """Get mask of longitude coordinates neighbors crossing the antimeridian.""" from scipy.ndimage import binary_dilation + # Initialize mask + n_y, n_x = lons.shape + mask = np.zeros((n_y - 1, n_x - 1)) + # Check vertical edges + row_idx, col_idx = np.where(np.abs(np.diff(lons, axis=0)) > 180) + col_idx = np.clip(col_idx - 1, 0, n_x - 1) + mask[row_idx, col_idx] = 1 + # Check horizontal edges + row_idx, col_idx = np.where(np.abs(np.diff(lons, axis=1)) > 180) + row_idx = np.clip(row_idx - 1, 0, n_y - 1) + mask[row_idx, col_idx] = 1 + # Buffer by 1 in all directions to avoid plotting cells neighbour to those crossing the antimeridian + # --> This should not be needed, but it's needed to avoid cartopy bugs ! + mask = binary_dilation(mask) + return mask + + +def get_antimeridian_mask_old(lons, buffer=True): + """Get mask of longitude coordinates neighbors crossing the antimeridian.""" + from scipy.ndimage import binary_dilation + # Check vertical edges row_idx, col_idx = np.where(np.abs(np.diff(lons, axis=0)) > 180) row_idx_rev, col_idx_rev = np.where(np.abs(np.diff(lons[::-1, :], axis=0)) > 180) @@ -264,20 +285,26 @@ def _plot_cartopy_pcolormesh( If rgb=True, expect rgb dimension to be at last position. x and y must represents longitude and latitudes. """ - # - Get x, y, and array to plot + # Get x, y, and array to plot da = da.compute() x = da[x].data y = da[y].data arr = da.data - # - Infill invalid value and add mask if necessary + # Infill invalid value and add mask if necessary x, y, arr = get_valid_pcolormesh_inputs(x, y, arr, rgb=rgb) - # - Ensure arguments + # Ensure arguments if rgb: add_colorbar = False - # - Mask cells crossing the antimeridian + # Compute coordinates of cell corners for pcolormesh quadrilateral mesh + # - This enable correct masking of cells crossing the antimeridian + from gpm_api.utils.area import _get_lonlat_corners + + x, y = _get_lonlat_corners(x, y) + + # Mask cells crossing the antimeridian # --> Here we assume not invalid coordinates anymore # --> Cartopy still bugs with several projections when data cross the antimeridian # --> This flag can be unset with gpm_api.config.set({"viz_hide_antimeridian_data": False}) @@ -306,7 +333,7 @@ def _plot_cartopy_pcolormesh( cmap.set_bad(bad) plot_kwargs["cmap"] = cmap - # - Add variable field with cartopy + # Add variable field with cartopy if not rgb: p = ax.pcolormesh( x, @@ -316,11 +343,11 @@ def _plot_cartopy_pcolormesh( **plot_kwargs, ) - # - Add RGB + # Add RGB else: p = _plot_rgb_pcolormesh(x, y, arr, ax=ax, **plot_kwargs) - # - Add colorbar + # Add colorbar # --> TODO: set axis proportion in a meaningful way ... if add_colorbar: _ = plot_colorbar(p=p, ax=ax, cbar_kwargs=cbar_kwargs)