From 83d0cdb8352b2cc2e933b1aa5d331af074230c0c Mon Sep 17 00:00:00 2001 From: Jordi Figueras i Ventura Date: Thu, 3 Nov 2022 16:57:15 +0100 Subject: [PATCH] ADD: Add add_grid_lines as an argument for plotting (#1311) * Added keyword add_lines to plot_grid in gridmapdisplay * added keyword add_lines into function plot_ppi_map in radarmapdisplay. General linting of gridmapdisplay and radarmapdisplay * changed keyword to add_grid_lines * Changed keyword add_lines to add_grid_lines --- pyart/graph/gridmapdisplay.py | 106 ++++++++++++++++++--------------- pyart/graph/radarmapdisplay.py | 38 +++++++----- 2 files changed, 80 insertions(+), 64 deletions(-) diff --git a/pyart/graph/gridmapdisplay.py b/pyart/graph/gridmapdisplay.py index f9f8a3701d..1ff64708f3 100644 --- a/pyart/graph/gridmapdisplay.py +++ b/pyart/graph/gridmapdisplay.py @@ -44,7 +44,8 @@ except ImportError: _LAMBERT_GRIDLINES = False -class GridMapDisplay(object): + +class GridMapDisplay(): """ A class for creating plots from a grid object using xarray with a cartopy projection. @@ -71,15 +72,15 @@ def __init__(self, grid, debug=False): if not _CARTOPY_AVAILABLE: raise MissingOptionalDependency( 'Cartopy is required to use GridMapDisplay but is not ' - + 'installed!') + 'installed!') if not _XARRAY_AVAILABLE: raise MissingOptionalDependency( 'Xarray is required to use GridMapDisplay but is not ' - + 'installed!') + 'installed!') if not _NETCDF4_AVAILABLE: raise MissingOptionalDependency( 'netCDF4 is required to use GridMapDisplay but is not ' - + 'installed!') + 'installed!') # set attributes self.grid = grid @@ -87,7 +88,7 @@ def __init__(self, grid, debug=False): self.mappables = [] self.fields = [] self.origin = 'origin' - + def plot_grid(self, field, level=0, vmin=None, vmax=None, norm=None, cmap=None, mask_outside=False, title=None, title_flag=True, axislabels=(None, None), @@ -95,7 +96,7 @@ def plot_grid(self, field, level=0, vmin=None, vmax=None, colorbar_label=None, colorbar_orient='vertical', ax=None, fig=None, lat_lines=None, lon_lines=None, projection=None, - embellish=True, ticks=None, ticklabs=None, + embellish=True, add_grid_lines=True, ticks=None, ticklabs=None, imshow=False, **kwargs): """ Plot the grid using xarray and cartopy. @@ -158,8 +159,10 @@ def plot_grid(self, field, level=0, vmin=None, vmax=None, Map projection supported by cartopy. Used for all subsequent calls to the GeoAxes object generated. Defaults to PlateCarree. embellish : bool - True by default. Set to False to supress drawinf of coastlines + True by default. Set to False to supress drawing of coastlines etc... Use for speedup when specifying shapefiles. + add_grid_lines : bool + True by default. Set to False to supress drawing of lat/lon lines Note that lat lon labels only work with certain projections. ticks : array Colorbar custom tick label locations. @@ -176,13 +179,6 @@ def plot_grid(self, field, level=0, vmin=None, vmax=None, vmin, vmax = common.parse_vmin_vmax(self.grid, field, vmin, vmax) cmap = common.parse_cmap(cmap, field) - if lon_lines is None: - lon_lines = np.linspace(np.around(ds.lon.min()-.1, decimals=2), - np.around(ds.lon.max()+.1, decimals=2), 5) - if lat_lines is None: - lat_lines = np.linspace(np.around(ds.lat.min()-.1, decimals=2), - np.around(ds.lat.max()+.1, decimals=2), 5) - # mask the data where outside the limits if mask_outside: data = ds[field].data @@ -228,7 +224,7 @@ def plot_grid(self, field, level=0, vmin=None, vmax=None, ax = plt.axes(projection=projection) # plot the grid using xarray - if norm is not None: # if norm is set do not override with vmin/vmax + if norm is not None: # if norm is set do not override with vmin/vmax vmin = vmax = None if imshow: @@ -252,6 +248,16 @@ def plot_grid(self, field, level=0, vmin=None, vmax=None, ax.add_feature(coastlines, linestyle='-', edgecolor='k', linewidth=2) + if add_grid_lines: + if lon_lines is None: + lon_lines = np.linspace( + np.around(ds.lon.min()-.1, decimals=2), + np.around(ds.lon.max()+.1, decimals=2), 5) + if lat_lines is None: + lat_lines = np.linspace( + np.around(ds.lat.min()-.1, decimals=2), + np.around(ds.lat.max()+.1, decimals=2), 5) + # labeling gridlines poses some difficulties depending on the # projection, so we need some projection-specific methods if ax.projection in [cartopy.crs.PlateCarree(), @@ -261,7 +267,7 @@ def plot_grid(self, field, level=0, vmin=None, vmax=None, xlocs=lon_lines, ylocs=lat_lines) ax.set_extent([lon_lines.min(), lon_lines.max(), lat_lines.min(), lat_lines.max()], - crs=projection) + crs=projection) ax.set_xticks(lon_lines, crs=projection) ax.set_yticks(lat_lines, crs=projection) @@ -296,8 +302,6 @@ def plot_grid(self, field, level=0, vmin=None, vmax=None, field=field, ax=ax, fig=fig, ticks=ticks, ticklabs=ticklabs) - return - def plot_crosshairs(self, lon=None, lat=None, linestyle='--', color='r', linewidth=2, ax=None): """ @@ -439,7 +443,7 @@ def plot_latitudinal_level(self, field, y_index, vmin=None, vmax=None, if len(z_1d) > 1: z_1d = _interpolate_axes_edges(z_1d) xd, yd = np.meshgrid(x_1d, z_1d) - if norm is not None: # if norm is set do not override with vmin, vmax + if norm is not None: # if norm is set do not override with vmin, vmax vmin = vmax = None pm = ax.pcolormesh( @@ -464,7 +468,6 @@ def plot_latitudinal_level(self, field, y_index, vmin=None, vmax=None, self.plot_colorbar(mappable=pm, label=colorbar_label, orientation=colorbar_orient, field=field, ax=ax, fig=fig, ticks=ticks, ticklabs=ticklabs) - return def plot_longitude_slice(self, field, lon=None, lat=None, **kwargs): """ @@ -580,7 +583,7 @@ def plot_longitudinal_level(self, field, x_index, vmin=None, vmax=None, z_1d = _interpolate_axes_edges(z_1d) xd, yd = np.meshgrid(y_1d, z_1d) - if norm is not None: # if norm is set do not override with vmin, vmax + if norm is not None: # if norm is set do not override with vmin, vmax vmin = vmax = None pm = ax.pcolormesh( @@ -605,8 +608,7 @@ def plot_longitudinal_level(self, field, x_index, vmin=None, vmax=None, self.plot_colorbar(mappable=pm, label=colorbar_label, orientation=colorbar_orient, field=field, ax=ax, fig=fig, ticks=ticks, ticklabs=ticklabs) - return - + def plot_cross_section(self, field, start, end, steps=100, interp_type='linear', x_axis=None, vmin=None, vmax=None, @@ -618,7 +620,8 @@ def plot_cross_section(self, field, start, end, ax=None, fig=None, ticks=None, ticklabs=None, **kwargs): """ - Plot a cross section through a set of given points (latitude, longitude). + Plot a cross section through a set of given points (latitude, + longitude). This uses the MetPy cross section interpolation function. @@ -629,14 +632,15 @@ def plot_cross_section(self, field, start, end, field : str Field to be plotted. start : tuple - A latitude-longitude pair designating the start point of the cross section - (units are degrees north and degrees east). + A latitude-longitude pair designating the start point of the cross + section (units are degrees north and degrees east). end : tuple - A latitude-longitude pair designating the end point of the cross section - (units are degrees north and degrees east). + A latitude-longitude pair designating the end point of the cross + section (units are degrees north and degrees east). steps: int - The number of points along the geodesic between the start and the end point - (including the end points) to use in the cross section. Defaults to 100. + The number of points along the geodesic between the start and the + end point (including the end points) to use in the cross section. + Defaults to 100. interp_type: str The interpolation method, either ‘linear’ or ‘nearest’ (see xarray.DataArray.interp() for details). Defaults to ‘linear’. @@ -701,21 +705,24 @@ def plot_cross_section(self, field, start, end, cmap = common.parse_cmap(cmap, field) # Convert the grid into an xarray object - ds= self.grid.to_xarray() + ds = self.grid.to_xarray() # Extract the proj parameters proj_params = self.grid.get_projparams() # Convert the projection information into cartopy - radar_crs = cartopy.crs.AzimuthalEquidistant(central_longitude=proj_params["lon_0"], - central_latitude=proj_params["lat_0"]) - - # Now, convert that to cf-compliant coordinate information and assign it to data + radar_crs = cartopy.crs.AzimuthalEquidistant( + central_longitude=proj_params["lon_0"], + central_latitude=proj_params["lat_0"]) + + # Now, convert that to cf-compliant coordinate information and assign + # it to data projection_info = radar_crs.to_cf() ds = ds.metpy.assign_crs(projection_info) # Calculate the cross section, which returns a dataset - ds = cross_section(ds, start, end, steps, interp_type).set_coords(('lat', 'lon')) + ds = cross_section(ds, start, end, steps, interp_type).set_coords( + ('lat', 'lon')) # Convert from meters to km for the different variables ds["z"] = ds["z"] / 1000 @@ -724,7 +731,7 @@ def plot_cross_section(self, field, start, end, if x_axis == 'y': ds["y"] = ds["y"] / 1000 ds.y.attrs["units"] = 'North South distance from radar (km)' - + if x_axis == 'x': ds["x"] = ds["x"] / 1000 ds.y.attrs["units"] = 'East West distance from radar (km)' @@ -733,8 +740,6 @@ def plot_cross_section(self, field, start, end, plot = ds[field].plot(y='z', x=x_axis, vmin=vmin, vmax=vmax, norm=norm, add_colorbar=False, ax=ax, cmap=cmap, **kwargs) - - self.mappables.append(plot) self.fields.append(field) @@ -743,7 +748,9 @@ def plot_cross_section(self, field, start, end, if title_flag: if title is None: - ax.set_title(common.generate_cross_section_title(self.grid, field, start, end)) + ax.set_title( + common.generate_cross_section_title( + self.grid, field, start, end)) else: ax.set_title(title) @@ -751,11 +758,10 @@ def plot_cross_section(self, field, start, end, self.plot_colorbar(mappable=plot, label=colorbar_label, orientation=colorbar_orient, field=field, ax=ax, fig=fig, ticks=ticks, ticklabs=ticklabs) - return - def plot_colorbar(self, mappable=None, orientation='horizontal', label=None, - cax=None, ax=None, fig=None, field=None, ticks=None, - ticklabs=None): + def plot_colorbar(self, mappable=None, orientation='horizontal', + label=None, cax=None, ax=None, fig=None, field=None, + ticks=None, ticklabs=None): """ Plot a colorbar. @@ -789,8 +795,7 @@ def plot_colorbar(self, mappable=None, orientation='horizontal', label=None, if mappable is None: if len(self.mappables) == 0: raise ValueError('mappable must be specified.') - else: - mappable = self.mappables[-1] + mappable = self.mappables[-1] if label is None: if len(self.fields) == 0: @@ -809,7 +814,6 @@ def plot_colorbar(self, mappable=None, orientation='horizontal', label=None, if ticklabs is not None: cb.set_ticklabels(ticklabs) cb.set_label(label) - return def _find_nearest_grid_indices(self, lon, lat): """ Find the nearest x, y grid indices for a given latitude and @@ -977,6 +981,7 @@ def cartopy_coastlines(self): # These methods are a hack to allow gridlines when the projection is lambert # https://nbviewer.jupyter.org/gist/ajdawson/dd536f786741e987ae4e + def find_side(ls, side): """ Given a shapely LineString which is assumed to be rectangular, return the @@ -989,6 +994,7 @@ def find_side(ls, side): 'top': [(minx, maxy), (maxx, maxy)]} return sgeom.LineString(points[side]) + def lambert_xticks(ax, ticks): """ Draw ticks on the bottom x-axis of a Lambert Conformal projection. """ def te(xy): @@ -1003,6 +1009,7 @@ def lc(t, n, b): ax.set_xticklabels([ax.xaxis.get_major_formatter()(xtick) for xtick in xticklabels]) + def lambert_yticks(ax, ticks): """ Draw ticks on the left y-axis of a Lambert Conformal projection. """ def te(xy): @@ -1017,8 +1024,11 @@ def lc(t, n, b): ax.set_yticklabels([ax.yaxis.get_major_formatter()(ytick) for ytick in yticklabels]) + def _lambert_ticks(ax, ticks, tick_location, line_constructor, tick_extractor): - """ Get the tick locations and labels for a Lambert Conformal projection. """ + """ + Get the tick locations and labels for a Lambert Conformal projection. + """ outline_patch = sgeom.LineString( ax.spines['geo'].get_path().vertices.tolist()) axis = find_side(outline_patch, tick_location) diff --git a/pyart/graph/radarmapdisplay.py b/pyart/graph/radarmapdisplay.py index 6f19fc201c..4ade26f11b 100644 --- a/pyart/graph/radarmapdisplay.py +++ b/pyart/graph/radarmapdisplay.py @@ -97,7 +97,6 @@ def __init__(self, radar, shift=(0.0, 0.0), grid_projection=None): self.ax = None self._x0 = None # x axis radar location in map coords (meters) self._y0 = None # y axis radar location in map coords (meters) - return def _check_ax(self): """ Check that a GeoAxes object exists, raise ValueError if not """ @@ -114,8 +113,9 @@ def plot_ppi_map( width=None, height=None, lon_0=None, lat_0=None, resolution='110m', shapefile=None, shapefile_kwargs=None, edges=True, gatefilter=None, - filter_transitions=True, embellish=True, raster=False, - ticks=None, ticklabs=None, alpha=None, edgecolors='face', **kwargs): + filter_transitions=True, embellish=True, add_grid_lines=True, + raster=False, ticks=None, ticklabs=None, alpha=None, + edgecolors='face', **kwargs): """ Plot a PPI volume sweep onto a geographic map. @@ -213,6 +213,8 @@ def plot_ppi_map( embellish: bool True by default. Set to False to supress drawing of coastlines etc.. Use for speedup when specifying shapefiles. + add_grid_lines : bool + True by default. Set to False to supress drawing of lat/lon lines Note that lat lon labels only work with certain projections. raster : bool False by default. Set to true to render the display as a raster @@ -232,10 +234,6 @@ def plot_ppi_map( # parse parameters vmin, vmax = parse_vmin_vmax(self._radar, field, vmin, vmax) cmap = parse_cmap(cmap, field) - if lat_lines is None: - lat_lines = np.arange(30, 46, 1) - if lon_lines is None: - lon_lines = np.arange(-110, -75, 1) lat_0 = self.loc[0] lon_0 = self.loc[1] @@ -284,13 +282,13 @@ def plot_ppi_map( + "axes with projection Lambert Conformal.", UserWarning) with warnings.catch_warnings(): - warnings.filterwarnings("ignore") + warnings.filterwarnings("ignore") ax = plt.axes(projection=projection) - if min_lon: + if min_lon is not None: ax.set_extent([min_lon, max_lon, min_lat, max_lat], crs=cartopy.crs.PlateCarree()) - elif width: + elif width is not None and height is not None: ax.set_extent([-width/2., width/2., -height/2., height/2.], crs=self.grid_projection) @@ -318,6 +316,12 @@ def plot_ppi_map( ax.coastlines(resolution=resolution) ax.add_feature(states_provinces, edgecolor='gray') + if add_grid_lines: + if lat_lines is None: + lat_lines = np.arange(30, 46, 1) + if lon_lines is None: + lon_lines = np.arange(-110, -75, 1) + # labeling gridlines poses some difficulties depending on the # projection, so we need some projection-spectific methods if ax.projection in [cartopy.crs.PlateCarree(), @@ -370,7 +374,6 @@ def plot_ppi_map( ax=ax, ticks=ticks, ticklabs=ticklabs) # keep track of this GeoAxes object for later self.ax = ax - return def plot_point(self, lon, lat, symbol='ro', label_text=None, label_offset=(None, None), **kwargs): @@ -401,7 +404,7 @@ def plot_point(self, lon, lat, symbol='ro', label_text=None, lon_offset = 0.01 if lat_offset is None: lat_offset = 0.01 - if 'transform' not in kwargs.keys(): + if 'transform' not in kwargs: kwargs['transform'] = cartopy.crs.PlateCarree() self.ax.plot(lon, lat, symbol, **kwargs) if label_text is not None: @@ -429,11 +432,12 @@ def plot_line_geo(self, line_lons, line_lats, line_style='r-', **kwargs): """ self._check_ax() - if 'transform' not in kwargs.keys(): + if 'transform' not in kwargs: kwargs['transform'] = cartopy.crs.PlateCarree() self.ax.plot(line_lons, line_lats, line_style, **kwargs) - def plot_line_xy(self, line_x, line_y, color='r', line_style='-', **kwargs): + def plot_line_xy(self, line_x, line_y, color='r', line_style='-', + **kwargs): """ Plot a line segments on the current map given radar x, y values. @@ -453,7 +457,7 @@ def plot_line_xy(self, line_x, line_y, color='r', line_style='-', **kwargs): """ self._check_ax() - if 'transform' not in kwargs.keys(): + if 'transform' not in kwargs: kwargs['transform'] = self.grid_projection self.ax.plot(line_x, line_y, color, line_style, **kwargs) @@ -539,7 +543,9 @@ def lc(t, n, b): def _lambert_ticks(ax, ticks, tick_location, line_constructor, tick_extractor): - """ Get the tick locations and labels for a Lambert Conformal projection. """ + """ + Get the tick locations and labels for a Lambert Conformal projection. + """ outline_patch = sgeom.LineString( ax.spines['geo'].get_path().vertices.tolist()) axis = find_side(outline_patch, tick_location)