Helper function: toggle_extent_bounds (#543)

* coord_on_land: fix buffer for global grids

* u_coord.toggle_extent_bounds

* No u_coord prefix in u_coord module

Co-authored-by: Emanuel Schmid <[email protected]>

climada/hazard/base.py

@@ -848,8 +848,9 @@ def select_tight(self, buffer=NEAREST_NEIGHBOR_THRESHOLD/ONE_LAT_KM,
             cent_nz = (self.fraction != 0).sum(axis=0).nonzero()[1]
         lon_nz = self.centroids.lon[cent_nz]
         lat_nz = self.centroids.lat[cent_nz]
-        ext = u_coord.latlon_bounds(lat=lat_nz, lon=lon_nz, buffer=buffer)
-        return self.select(extent=(ext[0], ext[2], ext[1], ext[3]))
+        return self.select(extent=u_coord.toggle_extent_bounds(
+            u_coord.latlon_bounds(lat=lat_nz, lon=lon_nz, buffer=buffer)
+        ))
 
     def local_exceedance_inten(self, return_periods=(25, 50, 100, 250)):
         """Compute exceedance intensity map for given return periods.

climada/hazard/tc_tracks.py

@@ -1142,8 +1142,7 @@ def get_extent(self, deg_buffer=0.1):
         -------
         extent : tuple (lon_min, lon_max, lat_min, lat_max)
         """
-        bounds = self.get_bounds(deg_buffer=deg_buffer)
-        return (bounds[0], bounds[2], bounds[1], bounds[3])
+        return u_coord.toggle_extent_bounds(self.get_bounds(deg_buffer=deg_buffer))
 
     @property
     def extent(self):

climada/hazard/test/test_tc_tracks.py

@@ -551,11 +551,10 @@ def test_get_extent(self):
         storms = ['1988169N14259', '2002073S16161', '2002143S07157']
         tc_track = tc.TCTracks.from_ibtracs_netcdf(storm_id=storms, provider=["usa", "bom"])
         bounds = (153.585022, -23.200001, 258.714996, 17.514986)
-        extent = (bounds[0], bounds[2], bounds[1], bounds[3])
         bounds_buf = (153.485022, -23.300001, 258.814996, 17.614986)
         np.testing.assert_array_almost_equal(tc_track.bounds, bounds)
         np.testing.assert_array_almost_equal(tc_track.get_bounds(deg_buffer=0.1), bounds_buf)
-        np.testing.assert_array_almost_equal(tc_track.extent, extent)
+        np.testing.assert_array_almost_equal(tc_track.extent, u_coord.toggle_extent_bounds(bounds))
 
     def test_generate_centroids(self):
         """Test centroids generation feature."""

climada/util/coordinates.py

@@ -243,6 +243,27 @@ def latlon_bounds(lat, lon, buffer=0.0):
     lon_min, lon_max = lon_bounds(lon, buffer=buffer)
     return (lon_min, max(lat.min() - buffer, -90), lon_max, min(lat.max() + buffer, 90))
 
+
+def toggle_extent_bounds(bounds_or_extent):
+    """Convert between the "bounds" and the "extent" description of a bounding box
+
+    The difference between the two conventions is in the order in which the bounds for each
+    coordinate direction are given. To convert from one description to the other, the two central
+    entries of the 4-tuple are swapped. Hence, the conversion is symmetric.
+
+    Parameters
+    ----------
+    bounds_or_extent : tuple (a, b, c, d)
+        Bounding box of the given points in "bounds" (or "extent") convention.
+
+    Returns
+    -------
+    extent_or_bounds : tuple (a, c, b, d)
+        Bounding box of the given points in "extent" (or "bounds") convention.
+    """
+    return (bounds_or_extent[0], bounds_or_extent[2], bounds_or_extent[1], bounds_or_extent[3])
+
+
 def dist_approx(lat1, lon1, lat2, lon2, log=False, normalize=True,
                 method="equirect", units='km'):
     """Compute approximation of geodistance in specified units
@@ -687,7 +708,7 @@ def coord_on_land(lat, lon, land_geom=None):
         bounds = latlon_bounds(lat, lons, buffer=delta_deg)
         # load land geometry with appropriate same extent
         land_geom = get_land_geometry(
-            extent=(bounds[0], bounds[2], bounds[1], bounds[3]),
+            extent=toggle_extent_bounds(bounds),
             resolution=10)
     elif not land_geom.is_empty:
         # ensure lon values are within extent of provided land_geom
@@ -792,14 +813,14 @@ def get_country_geometries(country_names=None, extent=None, resolution=10):
         # rewrap longitudes unless longitude extent is already normalized (within [-180, +180])
         lon_normalized = extent[0] >= -180 and extent[1] <= 180
         if lon_normalized:
-            bbox = box(extent[0], extent[2], extent[1], extent[3])
+            bbox = box(*toggle_extent_bounds(extent))
         else:
             # split the extent box into two boxes both within [-180, +180] in longitude
             lon_left, lon_right = lon_normalize(np.array(extent[:2]))
             extent_left = (lon_left, 180, extent[2], extent[3])
             extent_right = (-180, lon_right, extent[2], extent[3])
             bbox = shapely.ops.unary_union(
-                [box(e[0], e[2], e[1], e[3]) for e in [extent_left, extent_right]]
+                [box(*toggle_extent_bounds(e)) for e in [extent_left, extent_right]]
             )
         bbox = gpd.GeoSeries(bbox, crs=DEF_CRS)
         bbox = gpd.GeoDataFrame({'geometry': bbox}, crs=DEF_CRS)

climada/util/plot.py

@@ -599,7 +599,7 @@ def add_populated_places(axis, extent, proj=ccrs.PlateCarree(), fontsize=None):
                                          name='populated_places_simple')
 
     shp = shapereader.Reader(shp_file)
-    ext_pts = list(box(extent[0], extent[2], extent[1], extent[3]).exterior.coords)
+    ext_pts = list(box(*u_coord.toggle_extent_bounds(extent)).exterior.coords)
     ext_trans = [ccrs.PlateCarree().transform_point(pts[0], pts[1], proj)
                  for pts in ext_pts]
     for rec, point in zip(shp.records(), shp.geometries()):
@@ -632,7 +632,7 @@ def add_cntry_names(axis, extent, proj=ccrs.PlateCarree(), fontsize=None):
                                          name='admin_0_countries')
 
     shp = shapereader.Reader(shp_file)
-    ext_pts = list(box(extent[0], extent[2], extent[1], extent[3]).exterior.coords)
+    ext_pts = list(box(*u_coord.toggle_extent_bounds(extent)).exterior.coords)
     ext_trans = [ccrs.PlateCarree().transform_point(pts[0], pts[1], proj)
                  for pts in ext_pts]
     for rec, point in zip(shp.records(), shp.geometries()):

climada/util/test/test_coordinates.py

@@ -258,6 +258,10 @@ def test_latlon_bounds(self):
         bounds = u_coord.latlon_bounds(lat, lon, buffer=1)
         self.assertEqual(bounds, (-180, -90, 180, 90))
 
+    def test_toggle_extent_bounds(self):
+        """Test the conversion between 'extent' and 'bounds'"""
+        self.assertEqual(u_coord.toggle_extent_bounds((0, -1, 1, 3)), (0, 1, -1, 3))
+
     def test_geosph_vector(self):
         """Test conversion from lat/lon to unit vector on geosphere"""
         data = np.array([[0, 0], [-13, 179]], dtype=np.float64)