Plot NaNs in gray color and use log color bar when appropriate (#929)

* In the plot function geo_im_from_array, NaN values in the data will be plotted in gray.
  Before, NaN value were not plotted (i.e. transparent), making them indistinguishable
  from plot regions for which there is no data (no centroids).
* In the plot function plot_from_gdf, the colorbar with will be shown on a logarithmic scale
   if a) the gdf is about return periods or impacts, b) there are no zeros in the data,
   c) the span of the data's values are at least two orders of magnitude.
* Rename 'subplots_from_gdf' to 'plots_from_gdf'.
---------

Co-authored-by: Lukas Riedel <[email protected]>

Author: ValentinGebhart

CHANGELOG.md

@@ -14,6 +14,9 @@ Code freeze date: YYYY-MM-DD
 
 ### Changed
 
+- In `climada.util.plot.geo_im_from_array`, NaNs are plotted in gray while cells with no centroid are not plotted [#929](https://github.com/CLIMADA-project/climada_python/pull/929)
+- Renamed `climada.util.plot.subplots_from_gdf` to `climada.util.plot.plot_from_gdf` [#929](https://github.com/CLIMADA-project/climada_python/pull/929)
+
 ### Fixed
 
 ### Deprecated
@@ -89,6 +92,8 @@ CLIMADA tutorials. [#872](https://github.com/CLIMADA-project/climada_python/pull
 - climada.hazard.centroids.centr.Centroids.to_default_crs
 - climada.hazard.centroids.centr.Centroids.write_csv
 - climada.hazard.centroids.centr.Centroids.write_excel
+- climada.hazard.local_return_period [#898](https://github.com/CLIMADA-project/climada_python/pull/898)
+- climada.util.plot.subplots_from_gdf [#898](https://github.com/CLIMADA-project/climada_python/pull/898)
 
 ### Deprecated
 

climada/util/plot.py

@@ -309,10 +309,15 @@ def geo_im_from_array(array_sub, coord, var_name, title,
     if not proj:
         mid_lon = 0.5 * sum(extent[:2])
         proj = ccrs.PlateCarree(central_longitude=mid_lon)
-    if 'vmin' not in kwargs:
-        kwargs['vmin'] = np.nanmin(array_sub)
-    if 'vmax' not in kwargs:
-        kwargs['vmax'] = np.nanmax(array_sub)
+
+    if "norm" in kwargs:
+        min_value = kwargs["norm"].vmin
+    else:
+        kwargs['vmin'] = kwargs.get("vmin", np.nanmin(array_sub))
+        min_value = kwargs['vmin']
+        kwargs['vmax'] = kwargs.get("vmax", np.nanmax(array_sub))
+    min_value = min_value/2 if min_value > 0 else min_value-1
+
     if axes is None:
         proj_plot = proj
         if isinstance(proj, ccrs.PlateCarree):
@@ -328,8 +333,10 @@ def geo_im_from_array(array_sub, coord, var_name, title,
     if not isinstance(axes, np.ndarray):
         axes_iter = np.array([[axes]])
 
-    if 'cmap' not in kwargs:
-        kwargs['cmap'] = CMAP_RASTER
+    # prepare colormap
+    cmap = plt.get_cmap(kwargs.pop("cmap", CMAP_RASTER))
+    cmap.set_bad("gainsboro")  # For NaNs and infs
+    cmap.set_under("white", alpha=0)  # For values below vmin
 
     # Generate each subplot
     for array_im, axis, tit, name in zip(list_arr, axes_iter.flatten(), list_tit, list_name):
@@ -340,8 +347,11 @@ def geo_im_from_array(array_sub, coord, var_name, title,
             grid_x, grid_y = np.mgrid[
                 extent[0]: extent[1]: complex(0, RESOLUTION),
                 extent[2]: extent[3]: complex(0, RESOLUTION)]
-            grid_im = griddata((coord[:, 1], coord[:, 0]), array_im,
-                               (grid_x, grid_y))
+            grid_im = griddata(
+                (coord[:, 1], coord[:, 0]),
+                array_im,
+                (grid_x, grid_y),
+                fill_value=min_value)
         else:
             grid_x = coord[:, 1].reshape((width, height)).transpose()
             grid_y = coord[:, 0].reshape((width, height)).transpose()
@@ -359,8 +369,14 @@ def geo_im_from_array(array_sub, coord, var_name, title,
         # Create colormesh, colorbar and labels in axis
         cbax = make_axes_locatable(axis).append_axes('right', size="6.5%",
                                                      pad=0.1, axes_class=plt.Axes)
-        img = axis.pcolormesh(grid_x - mid_lon, grid_y, np.squeeze(grid_im),
-                              transform=proj, **kwargs)
+        img = axis.pcolormesh(
+            grid_x - mid_lon,
+            grid_y,
+            np.squeeze(grid_im),
+            transform=proj,
+            cmap=cmap,
+            **kwargs
+        )
         cbar = plt.colorbar(img, cax=cbax, orientation='vertical')
         cbar.set_label(name)
         axis.set_title("\n".join(wrap(tit)))
@@ -876,7 +892,7 @@ def multibar_plot(ax, data, colors=None, total_width=0.8, single_width=1,
     if legend:
         ax.legend(bars, data.keys())
 
-def subplots_from_gdf(
+def plot_from_gdf(
         gdf: gpd.GeoDataFrame,
         colorbar_name: str = None,
         title_subplots: callable = None,
@@ -917,7 +933,7 @@ def subplots_from_gdf(
     # check if inputs are correct types
     if not isinstance(gdf, gpd.GeoDataFrame):
         raise ValueError("gdf is not a GeoDataFrame")
-    gdf = gdf[['geometry', *[col for col in gdf.columns if col != 'geometry']]]
+    gdf_values = gdf.drop(columns='geometry').values.T
 
     # read meta data for fig and axis labels
     if not isinstance(colorbar_name, str):
@@ -927,23 +943,32 @@ def subplots_from_gdf(
         print("Unknown subplot-title-generation function. Subplot titles will be column names.")
         title_subplots = lambda cols: [f"{col}" for col in cols]
 
-    # change default plot kwargs if plotting return periods
-    if colorbar_name.strip().startswith('Return Period'):
-        if 'cmap' not in kwargs.keys():
-            kwargs.update({'cmap': 'viridis_r'})
-        if 'norm' not in kwargs.keys():
-            kwargs.update(
-                {'norm': mpl.colors.LogNorm(
-                    vmin=gdf.values[:,1:].min(), vmax=gdf.values[:,1:].max()
-                    ),
-                'vmin': None, 'vmax': None}
-            )
+    # use log colorbar for return periods and impact
+    if (
+        colorbar_name.strip().startswith(('Return Period', 'Impact')) and
+        'norm' not in kwargs.keys() and
+        # check if there are no zeros values in gdf
+        not np.any(gdf_values == 0) and
+        # check if value range too small for logarithmic colorscale
+        (np.log10(np.nanmax(gdf_values)) - np.log10(np.nanmin(gdf_values))) > 2
+    ):
+        kwargs.update(
+            {'norm': mpl.colors.LogNorm(
+                vmin=np.nanmin(gdf_values), vmax=np.nanmax(gdf_values)
+                ),
+            'vmin': None, 'vmax': None}
+        )
+
+    # use inverted color bar for return periods
+    if (colorbar_name.strip().startswith('Return Period') and
+        'cmap' not in kwargs.keys()):
+        kwargs.update({'cmap': 'viridis_r'})
 
     axis = geo_im_from_array(
-        gdf.values[:,1:].T,
+        gdf_values,
         gdf.geometry.get_coordinates().values[:,::-1],
         colorbar_name,
-        title_subplots(gdf.columns[1:]),
+        title_subplots(np.delete(gdf.columns, np.where(gdf.columns == 'geometry'))),
         smooth=smooth,
         axes=axis,
         figsize=figsize,

climada/util/test/test_plot.py

@@ -127,7 +127,7 @@ def test_geo_bin_from_array(self):
         plt.close()
 
     def test_geo_im_from_array(self):
-        values = np.array([1, 2.0, 5, 1])
+        values = np.array([1, 2.0, 5, np.nan])
         coord = np.array([[-17, 178], [-10, 180], [-27, 175], [-16, 186]])
         var_name = 'test'
         title = 'test'
@@ -137,8 +137,8 @@ def test_geo_im_from_array(self):
                       proj=projection, smooth=True, axes=None, figsize=(9, 13), cmap=cmap)
         self.assertEqual(var_name, ax.get_title())
         colorbar = next(x.colorbar for x in ax.collections if x.colorbar)
-        self.assertAlmostEqual(np.max(values), colorbar.vmax)
-        self.assertAlmostEqual(np.min(values), colorbar.vmin)
+        self.assertAlmostEqual(np.nanmax(values), colorbar.vmax)
+        self.assertAlmostEqual(np.nanmin(values), colorbar.vmin)
         self.assertEqual(cmap, ax.collections[0].cmap.name)
         plt.close()
 
@@ -147,20 +147,38 @@ def test_geo_im_from_array(self):
                       proj=projection, smooth=True, axes=None, figsize=(9, 13), cmap=cmap)
         self.assertEqual(var_name, ax.get_title())
         colorbar = next(x.colorbar for x in ax.collections if x.colorbar)
-        self.assertAlmostEqual(np.max(values), colorbar.vmax)
-        self.assertAlmostEqual(np.min(values), colorbar.vmin)
+        self.assertAlmostEqual(np.nanmax(values), colorbar.vmax)
+        self.assertAlmostEqual(np.nanmin(values), colorbar.vmin)
         self.assertEqual(cmap, ax.collections[0].cmap.name)
         plt.close()
 
-    def test_subplots_from_gdf(self):
+    def test_plot_from_gdf_no_log(self):
+        """test plot_from_gdf() with linear color bar (because there is a 0 in data)"""
+        return_periods = gpd.GeoDataFrame(
+            data = ((2., 5.), (0., 6.), (None, 2.), (1., 1000.)),
+            columns = ('10.0', '20.0')
+        )
+        return_periods['geometry'] = (Point(45., 26.), Point(46., 26.), Point(45., 27.), Point(46., 27.))
+        colorbar_name = 'Return Periods (Years)'
+        title_subplots = lambda cols: [f'Threshold Intensity: {col} m/s' for col in cols]
+        (axis1, axis2) = u_plot.plot_from_gdf(
+            return_periods, 
+            colorbar_name=colorbar_name, 
+            title_subplots=title_subplots)
+        self.assertEqual('Threshold Intensity: 10.0 m/s', axis1.get_title())
+        self.assertEqual('Threshold Intensity: 20.0 m/s', axis2.get_title())
+        plt.close()
+
+    def test_plot_from_gdf_log(self):
+        """test plot_from_gdf() with log color bar)"""
         return_periods = gpd.GeoDataFrame(
-            data = ((2., 5.), (3., 6.), (None, 2.), (1., 7.)),
+            data = ((2., 5.), (3., 6.), (None, 2.), (1., 1000.)),
             columns = ('10.0', '20.0')
         )
         return_periods['geometry'] = (Point(45., 26.), Point(46., 26.), Point(45., 27.), Point(46., 27.))
         colorbar_name = 'Return Periods (Years)'
         title_subplots = lambda cols: [f'Threshold Intensity: {col} m/s' for col in cols]
-        (axis1, axis2) = u_plot.subplots_from_gdf(
+        (axis1, axis2) = u_plot.plot_from_gdf(
             return_periods, 
             colorbar_name=colorbar_name, 
             title_subplots=title_subplots)

doc/tutorial/climada_hazard_Hazard.ipynb

@@ -577,7 +577,7 @@
     "<a id='Part5'></a> \n",
     "## Part 5: Visualize Hazards\n",
     "\n",
-    "There are three different plot functions: `plot_intensity()`, `plot_fraction()`, and `plot_rp_intensity()`. Depending on the inputs, different properties can be visualized. Check the documentation of the functions. Using the function `local_return_period()` and the util function `subplots_from_gdf()`, one can plot local return periods for specific hazard intensities."
+    "There are three different plot functions: `plot_intensity()`, `plot_fraction()`, and `plot_rp_intensity()`. Depending on the inputs, different properties can be visualized. Check the documentation of the functions. Using the function `local_return_period()` and the util function `plot_from_gdf()`, one can plot local return periods for specific hazard intensities."
    ]
   },
   {
@@ -759,8 +759,8 @@
     "\n",
     "# 5. tropical cyclone return period maps for the threshold intensities [30, 40]\n",
     "return_periods, label, column_label = haz_tc_fl.local_return_period([30, 40])\n",
-    "from climada.util.plot import subplots_from_gdf\n",
-    "subplots_from_gdf(return_periods, colorbar_name=label, title_subplots=column_label)\n",
+    "from climada.util.plot import plot_from_gdf\n",
+    "plot_from_gdf(return_periods, colorbar_name=label, title_subplots=column_label)\n",
     "\n",
     "# 6. intensities of all the events in centroid with id 50\n",
     "haz_tc_fl.plot_intensity(centr=50)\n",