Remove implicit use of CRS in code base

climada/hazard/trop_cyclone/trop_cyclone_windfields.py

@@ -902,6 +902,11 @@ def _coriolis_parameter(lat: np.ndarray) -> np.ndarray:
     cp : np.ndarray of same shape as input
         Coriolis parameter.
     """
+    if not u_coord.check_if_geo_coords(lat, 0):
+        raise ValueError(
+            "Input lat and lon coordinates are not geographic "
+            "or have total extents > 180° for lat or > 360° for lon."
+        )
     return 2 * V_ANG_EARTH * np.sin(np.radians(np.abs(lat)))
 
 

climada/util/coordinates.py

@@ -85,6 +85,24 @@
 are not considered."""
 
 
+def check_if_geo_coords(lat, lon):
+    """
+    Check if latitude and longitude arrays are likely in geographic coordinates,
+    testing if min/max values are within -90 to 90 for latitude and -180 to 540
+    for longitude.
+    Returns True if lat/lon ranges are in the geographic coordinates range, otherwise False.
+    """
+    lat = np.array(lat)
+    lon = np.array(lon)
+
+    # Check if latitude is within -90 to 90 and longitude is within -180 to 540
+    # and extent are smaller than 180 and 360 respectively
+    test = (
+        lat.min() >= -91 and lat.max() <= 91 and lon.min() >= -181 and lon.max() <= 541
+    ) and ((lat.max() - lat.min()) <= 181 and (lon.max() - lon.min()) <= 361)
+    return bool(test)
+
+
 def latlon_to_geosph_vector(lat, lon, rad=False, basis=False):
     """Convert lat/lon coodinates to radial vectors (on geosphere)
 
@@ -448,7 +466,12 @@ def get_gridcellarea(lat, resolution=0.5, unit="ha"):
     unit: string, optional
         unit of the output area (default: ha, alternatives: m2, km2)
     """
-
+    # first check that lat is in geographic coordinates
+    if not check_if_geo_coords(lat, 0):
+        raise ValueError(
+            "Input lat and lon coordinates are not geographic "
+            "or have total extents > 180° for lat or > 360° for lon."
+        )
     if unit == "m2":
         area = (ONE_LAT_KM * resolution) ** 2 * np.cos(np.deg2rad(lat)) * 1000000
     elif unit == "km2":
@@ -507,6 +530,11 @@ def convert_wgs_to_utm(lon, lat):
     epsg_code : int
         EPSG code of UTM projection.
     """
+    if not check_if_geo_coords(lat, lon):
+        raise ValueError(
+            "Input lat and lon coordinates are not geographic "
+            "or have total extents > 180° for lat or  > 360° for lon."
+        )
     epsg_utm_base = 32601 + (0 if lat >= 0 else 100)
     return epsg_utm_base + (math.floor((lon + 180) / 6) % 60)
 
@@ -564,6 +592,11 @@ def dist_to_coast(coord_lat, lon=None, highres=False, signed=False):
             raise ValueError(
                 f"Mismatching input coordinates size: {lat.size} != {lon.size}"
             )
+    if not check_if_geo_coords(lat, lon):
+        raise ValueError(
+            "Input lat and lon coordinates are not geographic "
+            "or have total extents > 180° for lat or  > 360° for lon."
+        )
     return dist_to_coast_nasa(lat, lon, highres=highres, signed=signed)
 
 
@@ -679,6 +712,11 @@ def coord_on_land(lat, lon, land_geom=None):
         )
     if lat.size == 0:
         return np.empty((0,), dtype=bool)
+    if not check_if_geo_coords(lat, lon):
+        raise ValueError(
+            "Input lat and lon coordinates are not geographic "
+            "or have total extents > 180° for lat or > 360° for lon."
+        )
     delta_deg = 1
     lons = lon.copy()
     if land_geom is None:
@@ -796,9 +834,10 @@ def get_country_geometries(
         out = out[country_mask]
 
     if extent:
-        if extent[1] - extent[0] > 360:
+        if not check_if_geo_coords(extent[2:], extent[:2]):
             raise ValueError(
-                f"longitude extent range is greater than 360: {extent[0]} to {extent[1]}"
+                "Input lat and lon coordinates are not geographic "
+                "or have total extents > 180° for lat or > 360° for lon."
             )
 
         if extent[1] < extent[0]:
@@ -1089,6 +1128,15 @@ def match_coordinates(
 
         # assign remaining coordinates to their geographically nearest neighbor
         if threshold > 0 and exact_assign_idx.size != coords_view.size:
+            # check that coords are geographic before proceeding to nearest neighbor search
+            if not (
+                check_if_geo_coords(coords[:, 0], coords[:, 1])
+                and check_if_geo_coords(coords_to_assign[:, 0], coords_to_assign[:, 1])
+            ):
+                raise ValueError(
+                    "Input lat and lon coordinates are not geographic "
+                    "or have total extents > 180° for lat or  > 360° for lon."
+                )
             not_assigned_idx_mask = assigned_idx == -1
             assigned_idx[not_assigned_idx_mask] = nearest_neighbor_funcs[distance](
                 coords_to_assign, coords[not_assigned_idx_mask], threshold, **kwargs
@@ -1596,6 +1644,12 @@ def get_country_code(lat, lon, gridded=False):
         return np.empty((0,), dtype=int)
     LOGGER.info("Setting region_id %s points.", str(lat.size))
     if gridded:
+        # first check that input lat lon are geographic
+        if not check_if_geo_coords(lat, lon):
+            raise ValueError(
+                "Input lat and lon coordinates are not geographic "
+                "or have total extents > 180° for lat or  > 360° for lon."
+            )
         base_file = u_hdf5.read(NATEARTH_CENTROIDS[150])
         meta, region_id = base_file["meta"], base_file["region_id"]
         transform = rasterio.Affine(*meta["transform"])

climada/util/test/test_coordinates.py

@@ -342,6 +342,23 @@ def data_arrays_resampling_demo():
 class TestFunc(unittest.TestCase):
     """Test auxiliary functions"""
 
+    def test_check_geo_coords(self):
+        """Test the check_if_geo_coords function"""
+        # test correct geographical coords
+        lat, lon = np.array([0, -2, 5]), np.array([-179, 175, 178])
+        self.assertEqual(u_coord.check_if_geo_coords(lat, lon), True)
+        # test wrong lat
+        lat, lon = np.array([0, 200, 5]), np.array([-179, 175, 178])
+        self.assertEqual(u_coord.check_if_geo_coords(lat, lon), False)
+        # test wrong lon
+        lat, lon = np.array([0, -2, 5]), np.array([-400, 175, 176])
+        self.assertEqual(u_coord.check_if_geo_coords(lat, lon), False)
+        lat, lon = np.array([0, -2, 5]), np.array([700, 1000, 800])
+        self.assertEqual(u_coord.check_if_geo_coords(lat, lon), False)
+        # test wrong extent
+        lat, lon = np.array([0, -2, 5]), np.array([-179, 175, 370])
+        self.assertEqual(u_coord.check_if_geo_coords(lat, lon), False)
+
     def test_lon_normalize(self):
         """Test the longitude normalization function"""
         data = np.array([-180, 20.1, -30, 190, -350])
@@ -1460,12 +1477,14 @@ def test_get_country_geometries_fail(self):
         with self.assertRaises(ValueError) as cm:
             u_coord.get_country_geometries(extent=(-20, 350, 0, 0))
         self.assertIn(
-            "longitude extent range is greater than 360: -20 to 350", str(cm.exception)
+            "Input lat and lon coordinates are not geographic "
+            "or have total extents > 180° for lat or > 360° for lon.",
+            str(cm.exception),
         )
         with self.assertRaises(ValueError) as cm:
-            u_coord.get_country_geometries(extent=(350, -20, 0, 0))
+            u_coord.get_country_geometries(extent=(340, -20, 0, 0))
         self.assertIn(
-            "longitude extent at the left (350) is larger "
+            "longitude extent at the left (340) is larger "
             "than longitude extent at the right (-20)",
             str(cm.exception),
         )