diff --git a/mapgenerator/plotting/definitions.py b/mapgenerator/plotting/definitions.py index cde97006389ffd14276d4cc647044b1b41306e49..97783d0b7834b152f96fc373863efa58732763f1 100644 --- a/mapgenerator/plotting/definitions.py +++ b/mapgenerator/plotting/definitions.py @@ -16,6 +16,8 @@ from datetime import timedelta import os import os.path import copy +from cartopy.img_transform import mesh_projection, regrid +from cartopy.crs import PlateCarree import logging @@ -58,37 +60,30 @@ def do_interpolation(data, lon, lat): """ Takes a not regular grid (curvilinear, whatever kind) and interpolate to regular one derived from the first. Very demanding operation, try to avoid """ - # lon as second dimension - # lat as first dimension LOG.info("Irregular grid, performing interpolation ...") - lat_s = lat.shape[0] - lon_s = lon.shape[1] - reglon = np.linspace(lon.min(), lon.max(), lon_s) - reglat = np.linspace(lat.min(), lat.max(), lat_s) - LOG.info("LAT: %s, LON: %s, DATA: %s", reglon.shape, reglat.shape, - data.shape) - regx, regy = np.meshgrid(reglon, reglat) - result = np.empty(data.shape)*np.nan - if ma.is_masked(data): - data = data.filled(np.nan) - + projection = PlateCarree() + new_lon, new_lat, _ = mesh_projection( + projection, data.shape[-1], data.shape[-2], + # x_extents=(lon.min(), lon.max()), + # y_extents=(lat.min(), lat.max()) + ) + + def _interpolate(to_regrid): + return regrid(to_regrid, lon, lat, projection, projection, new_lon, new_lat) + + result = np.empty_like(data) + result.fill(np.nan) # assuming that lat and lon are the latest dimensions if len(data.shape) == 2: - result = griddata((lon.ravel(), lat.ravel()), data.ravel(), (regx, - regy)) + result = _interpolate(data) elif len(data.shape) == 3: for i in range(data.shape[0]): - tmp = griddata((lon.ravel(), lat.ravel()), data[i].ravel(), (regx, - regy)) - result[i] = tmp.reshape(result[i].shape) + result[i, :, :] = _interpolate(data[i]) elif len(data.shape) == 4: for i in range(data.shape[0]): for j in range(data.shape[1]): - tmp = griddata((lon.ravel(), lat.ravel()), data[i, j].ravel(), - (regx, regy)) - result[i, j, :, :] = tmp.reshape(result[i, j, :, :].shape) - # elif len(data.shape) == 5: - return np.ma.masked_where(np.isnan(result), result), regx, regy + result[i, j, :, :] = _interpolate(data[i, j]) + return result, new_lon[0], new_lat[:, 0] def extract_coords(self, var, coord='lons'): @@ -741,9 +736,9 @@ class MapGenerator(object): if cube.attributes.get('plot_name', None): show_name = cube.attributes['plot_name'] elif cube.long_name and len(cube.long_name) < 45: - show_name = cube.long_name + show_name = f"{cube.long_name} ({cube.var_name})" elif cube.standard_name and len(cube.standard_name) < 45: - show_name = cube.standard_name + show_name = f"{cube.standard_name} ({cube.var_name})" else: show_name = cube.var_name return f"{show_name} ({cube.units})" diff --git a/mapgenerator/plotting/plotmap.py b/mapgenerator/plotting/plotmap.py index 98c5aaf252ad371b64cf7072cdbc27f2eba50885..54d59940649a78f6ec26d829e9c3e659a48f37a9 100644 --- a/mapgenerator/plotting/plotmap.py +++ b/mapgenerator/plotting/plotmap.py @@ -6,9 +6,12 @@ import matplotlib as mpl import matplotlib.pyplot as plt +import matplotlib.ticker as mticker import cartopy.crs as ccrs import cartopy.feature as cfeature import cartopy.io.shapereader as shpreader + +from cartopy.util import add_cyclic_point from matplotlib.cbook import delete_masked_points from matplotlib.path import Path import numpy as np @@ -70,6 +73,7 @@ class PlotMap(MapCross, MapDrawOptions): self.zoom_level = None self.first_image = False # Is is the first image of the batch? self.map_name = None # name of the figure to map animation colors + self.draw_labels = None def __setattr__(self, key, value): LOG.debug("SETATTR: %s - %s", key, value) @@ -208,7 +212,7 @@ class PlotMap(MapCross, MapDrawOptions): cax, kwargs = mpl.colorbar.make_axes(self.mgaxis, location=location, pad=pad, - shrink=0.8) + shrink=0.9) cbar = self.mgplot.colorbar(mco, cax=cax, ax=self.mgaxis, @@ -218,7 +222,6 @@ class PlotMap(MapCross, MapDrawOptions): extend=self.extend, drawedges=drawedges, **kwargs) - cbar.ax.tick_params(labelsize=float(self.coordsopts[1])) for lin in cbar.ax.yaxis.get_ticklines(): lin.set_visible(False) @@ -327,9 +330,16 @@ class PlotMap(MapCross, MapDrawOptions): LOG.debug("1. GLON: %s, GLAT: %s", str(glon.shape), str(glat.shape)) # if curvilinear grid do interpolation, return already gridded coords - if not is_grid_regular(glon, glat): + if is_grid_regular(glon, glat): + if len(glon.shape)==2: + glon = glon[0] + glat = glat[:, 0] + else: map_data[0], glon, glat = do_interpolation(map_data[0], glon, glat) + + map_data[0], glon = add_cyclic_point(map_data[0], coord=glon) + LOG.info("2. GLON: %s, GLAT: %s", str(glon.shape), str(glat.shape)) xloc, yloc = glon, glat @@ -553,31 +563,28 @@ class PlotMap(MapCross, MapDrawOptions): edgecolor=str(self.countropts[1]), zorder=15) - draw_labels = bool(self.projection == 'PlateCarree') - grl = self.mgaxis.gridlines(xlocs=self.lon, - ylocs=self.lat, - crs=self._crs, - draw_labels=draw_labels, - linestyle='--', - linewidth=float(self.coordsopts[0]), - color=str(self.coordsopts[2]), - zorder=20) + if self.draw_labels is None: + self.draw_labels = bool(self.projection == 'PlateCarree') - if draw_labels: + + grl = self.mgaxis.gridlines( + xlocs=self.lon, + ylocs=self.lat, + draw_labels=self.draw_labels, + linestyle='--', + linewidth=float(self.coordsopts[0]), + color=str(self.coordsopts[2]), + zorder=20, + auto_inline=True, + ) + + if self.draw_labels and bool(self.projection == 'PlateCarree'): grl.xlabels_top = False grl.ylabels_right = False + if self.draw_labels: grl.xlabel_style = {'size': float(self.coordsopts[1])} grl.ylabel_style = grl.xlabel_style - else: - grl = self.mgaxis.gridlines(xlocs=self.lon, - ylocs=self.lat, - draw_labels=True, - linestyle='--', - linewidth=float( - self.coordsopts[0]), - color=str(self.coordsopts[2]), - zorder=20) # Change axes for colorbar if self.colorbar and not self.kml and not self.kmz and \ @@ -603,7 +610,6 @@ class PlotMap(MapCross, MapDrawOptions): fullname = "{}.{}".format(fig_name, self.filefmt) LOG.info("printing %s", fullname) if self.kml or self.kmz: - # self.mgplot.axes(frameon=0) if self.save: self.mgplot.savefig(fullname, bbox_inches='tight', @@ -616,6 +622,7 @@ class PlotMap(MapCross, MapDrawOptions): if self.save: self.mgplot.savefig(fullname, + frame_on = True, bbox_inches='tight', pad_inches=.2, dpi=self.dpi) @@ -724,15 +731,17 @@ class PlotMap(MapCross, MapDrawOptions): fig_name = "%s/%s" % (self.outdir, f_name) if self.subplot is None: plt.clf() - self.mgaxis = plt.axes(projection=self._crs) + self.mgaxis = plt.axes(projection=self._crs, frame_on=True) self.mgplot = plt.gcf() else: self.mgaxis = plt.subplot(self.subplot[0], self.subplot[1], self.subplot[2], - projection=self._crs) + projection=self._crs, + frame_on=True) self.mgplot = plt.gcf() - + self.mgaxis.patch.set_edgecolor('black') + self.mgaxis.patch.set_linewidth('1') self.gen_image_map( fig_name, grid, @@ -913,7 +922,7 @@ class PlotMap(MapCross, MapDrawOptions): cax = fig.add_axes([0.05, 0.80, 0.9, 0.15]) print("...", self.bounds, "...") cbar = mpl.colorbar.ColorbarBase( - cax, # =self.mgaxis, + cax, cmap=self.cmap, norm=self.norm, values=self.bounds,