Merge branch 'develop' into feature/interp_util_func

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/hazard/tc_tracks.py

@@ -1509,7 +1509,6 @@ def to_geodataframe(self, as_points=False, split_lines_antimeridian=True):
         return gdf
 
     @staticmethod
-    @numba.jit(forceobj=True)
     def _one_interp_data(track, time_step_h, land_geom=None):
         """Interpolate values of one track.
 

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",