fix issue #69 Warn.zeropadding for islands (#70)

* fix issue #69 Warn.zeropadding for islands

* Clearify description of Warn.zeropadding

Co-authored-by: leonie-villiger <87335708+leonie-villiger@users.noreply.github.com>

* Correct copying error in Warn.zeropadding

check_lat was calculated using the lon variables due to a copying error. This was corrected and it is now calculated correctly using the lat variables

Co-authored-by: leonie-villiger <87335708+leonie-villiger@users.noreply.github.com>

* Compare absolute values in check for grid resolution in Warn.zeropadding

* simply coord_rec creation in Warn.zeropadding

* correct lint in Warn.zeropadding

* include coord in test of Warn.zeropadding

---------

Co-authored-by: Thomas Roosli <thomas.roeoesli@meteoswiss.ch>
Co-authored-by: leonie-villiger <87335708+leonie-villiger@users.noreply.github.com>

climada_petals/engine/test/test_warn.py

@@ -333,6 +333,32 @@ def test_zeropadding(self):
         self.assertEqual(reduced_matrix[0, 2], 2)
         self.assertEqual(reduced_matrix[2, 2], 6)
         self.assertEqual(np.sum(reduced_matrix), np.sum(data))
+        self.assertEqual(reduced_matrix.size*2, coord.size)
+        np.testing.assert_array_equal(np.unique(coord[:,1]), np.arange(0,3,1))
+        np.testing.assert_array_equal(np.unique(coord[:,1]), np.arange(0,3,1))
+
+    def test_zeropadding_island(self):
+        row = np.array([1, 1, 8, 9])
+        col = np.array([3, 4, 7, 8])
+        data = np.array([20, 25, 28, 32])
+
+        reduced_matrix, coord = Warn.zeropadding(row, col, data)
+
+        self.assertEqual(reduced_matrix.shape, (9, 6))
+        self.assertEqual(reduced_matrix[0, 0], 20)
+        self.assertEqual(reduced_matrix[8, 5], 32)
+        self.assertEqual(np.sum(reduced_matrix), np.sum(data))
+        self.assertEqual(reduced_matrix.size*2, coord.size)
+        np.testing.assert_array_equal(np.unique(coord[:,0]), np.arange(1,10,1))
+        np.testing.assert_array_equal(np.unique(coord[:,1]), np.arange(3,9,1))
+
+    def test_zeropadding_resolution_error(self):
+        row = np.array([0, 0, 1.125, 2.3333, 2.3333, 2.3333])
+        col = np.array([0, 2.4444, 2.4444, 0, 1.375, 2.4444])
+        data = np.array([1, 2, 3, 4, 5, 6])
+
+        with self.assertRaises(ValueError):
+            reduced_matrix, coord = Warn.zeropadding(row, col, data)
 
     def test_plot_warning(self):
         """Plots ones with geo_scatteR_categorical"""

climada_petals/engine/warn.py

@@ -31,6 +31,7 @@
 import skimage
 
 from climada.util.plot import geo_scatter_categorical
+from climada.util.coordinates import get_resolution as u_get_resolution
 
 LOGGER = logging.getLogger(__name__)
 
@@ -461,10 +462,13 @@ def _group_cosmo_ensembles(ensembles, quant_nr):
         return single_map
 
     @staticmethod
-    def zeropadding(lat, lon, val):
-        """Produces a rectangle shaped map from a non-rectangular map (e.g., country). Therefore,
-        a rectangular around the countries' boundary is shaped and padded with zeros where no values
-        are defined.
+    def zeropadding(lat, lon, val, res_rel_error=0.01):
+        """Produces a rectangular shaped map from a non-rectangular map (e.g., country). For this,
+        a regular grid is created in the rectangle enclosing the non-rectangular map. The values are
+        mapped onto the regular grid according to their coordinates. The regular gird is filled with
+        zeros where no values are defined.
+        are defined. This only works if the lat lon values of the non-rectangular map can be
+        accurately represented on a grid with a regular resolution.
 
         Parameters
         ----------
@@ -474,6 +478,10 @@ def zeropadding(lat, lon, val):
             Longitudes of values of map.
         val : list
             Values of quantity of interest at every coordinate given.
+        res_rel_error : float
+            defines the relative error in terms of grid resolution by which the lat lon values
+            of the returned coord_rec can maximally differ from the provided lat lon values.
+            Default: 0.01
 
         Returns
         ----------
@@ -486,16 +494,43 @@ def zeropadding(lat, lon, val):
         lat = np.round(lat, decimals=12)
         lon = np.round(lon, decimals=12)
 
-        un_y = np.sort(np.unique(lat))
-        un_x = np.sort(np.unique(lon))
-
-        i = ((lat - min(lat)) / abs(un_y[1] - un_y[0])).astype(int)
-        j = ((lon - min(lon)) / abs(un_x[1] - un_x[0])).astype(int)
-        map_rec = np.zeros((len(np.unique(lat)), len(np.unique(lon))))
+        grid_res_lon, grid_res_lat = u_get_resolution(lon, lat)
+        # check if lat and lon can be represented accurately enough
+        # with a grid with grid resolution grid_res_lat and grid_res_lon
+        check_lat = (
+            np.abs(np.mod(np.diff(np.sort(np.unique(lat))-lat.min()), grid_res_lat).max())
+            <
+            np.abs(grid_res_lat * res_rel_error)
+            )
+        check_lon = (
+            np.abs(np.mod(np.diff(np.sort(np.unique(lon))-lon.min()), grid_res_lon).max())
+            <
+            np.abs(grid_res_lon * res_rel_error)
+            )
+        if not check_lon or not check_lat:
+            raise ValueError('The provided lat and lon values cannot be ' +
+                             'represented accurately enough by a regular ' +
+                             'grid. Provide different lat and lon or ' +
+                             'change res_rel_error.')
+        un_y = np.arange(
+            lat.min(),
+            lat.max()+abs(grid_res_lat),
+            abs(grid_res_lat)
+            )
+        un_x = np.arange(
+            lon.min(),
+            lon.max()+abs(grid_res_lon),
+            abs(grid_res_lon)
+            )
+
+        i = ((lat - min(lat)) / abs(grid_res_lat)).astype(int)
+        j = ((lon - min(lon)) / abs(grid_res_lon)).astype(int)
+        map_rec = np.zeros((len(un_y), len(un_x)))
         map_rec[i, j] = val
 
-        xx, yy = np.meshgrid(un_x, un_y)
-        coord_rec = np.vstack((yy.flatten(), xx.flatten())).transpose()
+        coord_rec = np.vstack(
+            [ar.flatten() for ar in np.meshgrid(un_y, un_x)]
+            ).transpose()
 
         return map_rec, coord_rec