util.coordinates: remove broken dist_to_coast function (#840)

* util.coordinates: remove dist_to_coast function

* Update climada/hazard/centroids/centr.py

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

* u_coord: dist_to_coast as wrapper for dist_to_coast_nasa

* Update climada/hazard/centroids/centr.py

Fix trailing whitespace

* Centroids.get_dist_coast: fix tests

---------

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

Author: tovogt

climada/hazard/centroids/centr.py

@@ -659,27 +659,37 @@ def get_elevation(self, topo_path):
         """
         return u_coord.read_raster_sample(topo_path, self.lat, self.lon)
 
-    def get_dist_coast(self, signed=False, precomputed=False):
+    def get_dist_coast(self, signed=False, precomputed=True):
         """Get dist_coast attribute for every pixel or point in meters.
 
+        The distances are read from a raster file containing precomputed distances (from NASA) at
+        0.01 degree (approximately 1 km) resolution.
+
         Parameters
         ----------
         signed : bool, optional
             If True, use signed distances (positive off shore and negative on land). Default: False
         precomputed : bool, optional
-            If True, use precomputed distances (from NASA). Default: False
+            Whether distances should be read from a pre-computed raster (True) or computed
+            on-the-fly (False). Default: True.
+
+            .. deprecated:: 5.0
+               Argument is ignored, because distances are not computed on-the-fly anymore.
 
         Returns
         -------
         dist : np.array
             (Signed) distance to coast in meters.
         """
+        if not precomputed:
+            LOGGER.warning(
+                "The `precomputed` argument is deprecated and will be removed in the future"
+                " because `get_dist_coast` always uses precomputed distances."
+            )
         ne_geom = self._ne_crs_geom()
-        if precomputed:
-            return u_coord.dist_to_coast_nasa(
-                ne_geom.y.values, ne_geom.x.values, highres=True, signed=signed)
-        LOGGER.debug('Computing distance to coast for %s centroids.', str(self.size))
-        return u_coord.dist_to_coast(ne_geom, signed=signed)
+        return u_coord.dist_to_coast_nasa(
+            ne_geom.y.values, ne_geom.x.values, highres=True, signed=signed,
+        )
 
     def get_meta(self, resolution=None):
         """Returns a meta raster based on the centroids bounds.

climada/hazard/centroids/test/test_centr.py

@@ -104,7 +104,7 @@ def test_init_pass(self):
         self.assertTrue(u_coord.equal_crs(self.centr.crs, DEF_CRS))
 
         # Creating Centroids with additional attributes
-        centroids = Centroids(lat=VEC_LAT, lon=VEC_LON, 
+        centroids = Centroids(lat=VEC_LAT, lon=VEC_LON,
                               region_id=REGION_ID, on_land=ON_LAND)
 
         # Checking additional attributes
@@ -126,7 +126,7 @@ def test_init_properties(self):
                       'on_land','region_id','crs',
                       'shape','size','total_bounds','coord']
         centroids = Centroids(lat=[],lon=[])
-        [self.assertTrue(hasattr(centroids,prop)) for prop in properties] 
+        [self.assertTrue(hasattr(centroids,prop)) for prop in properties]
 
     def test_init_kwargs(self):
         """ Test default crs and kwargs forwarding """
@@ -151,7 +151,7 @@ def test_init_kwargs(self):
 
     def test_from_meta_pass(self):
         expected_lon = np.array([-30.0, -20.0, -10.0]*3)
-        expected_lat = np.repeat([30.0, 20.0, 10.0],3)        
+        expected_lat = np.repeat([30.0, 20.0, 10.0],3)
         # Check metadata
         meta = dict(
             crs=DEF_CRS,
@@ -233,7 +233,7 @@ def test_from_pnt_bounds(self):
         np.testing.assert_allclose([10.0, 10.0, 9.8], centr.lat[[0, 1, width]], atol=0.1)
         # generally we assume that from_meta does not set region_ids and on_land flags
         self.assertFalse(centr.region_id)
-        self.assertFalse(centr.on_land)             
+        self.assertFalse(centr.on_land)
 
 class TestCentroidsTransformation(unittest.TestCase):
     """ Test class for coordinate transformations of Centroid objects
@@ -247,7 +247,7 @@ def setUp(self):
         self.centr = Centroids(lat=VEC_LAT,lon=VEC_LON,crs=TEST_CRS)
 
     def test_to_default_crs(self):
-        # Creating Centroids with non-default CRS and 
+        # Creating Centroids with non-default CRS and
         # inplace transformation afterwards
         centroids = Centroids(lat=VEC_LAT, lon=VEC_LON, crs=ALT_CRS)
         self.assertTrue(u_coord.equal_crs(centroids.crs, ALT_CRS))
@@ -329,7 +329,7 @@ def test_set_on_land_implementationerror(self):
 
         with self.assertRaises(NotImplementedError):
             centroids.set_on_land(source='satellite',overwrite=True)
-     
+
     def test_set_on_land_raster(self):
         """Test set_on_land"""
         centr_ras = Centroids.from_raster_file(HAZ_DEMO_FL, window=Window(0, 0, 50, 60))
@@ -736,7 +736,7 @@ def test_append_pass(self):
     def test_append(self):
         lat2,lon2 = np.array([6,7,8,9,10]),np.array([6,7,8,9,10])
         newcentr = Centroids(lat=lat2,lon=lon2)
-        newcentr.append(self.centr) 
+        newcentr.append(self.centr)
         self.assertTrue(newcentr.size == len(self.centr.lon)+len(lon2))
         np.testing.assert_array_equal(newcentr.lon,np.concatenate([lon2,self.centr.lon]))
         np.testing.assert_array_equal(newcentr.lat,np.concatenate([lat2,self.centr.lat]))
@@ -764,7 +764,7 @@ def test_remove_duplicate_pass(self):
 
     def test_remove_duplicates_dif_on_land(self):
         ### We currently expect that only the geometry of the gdf defines duplicates.
-        ### If one geometry is duplicated with differences in other attributes e.g. on_land 
+        ### If one geometry is duplicated with differences in other attributes e.g. on_land
         ### they get removed nevertheless. Only the first occurrence will be part of the new object
         ### this test is only here to guarantee this behaviour
         lat, lon = np.array([0,0,1,2,3,4,5]),np.array([0,0,1,2,3,4,5])
@@ -926,19 +926,18 @@ def test_get_closest_point(self):
         self.assertEqual(idx, 1)
 
     def test_dist_coast_pass(self):
-        """Test set_dist_coast"""
+        """Test get_dist_coast"""
         dist_coast = self.centr.get_dist_coast()
         # Just checking that the output doesnt change over time.
         REF_VALUES = np.array([
-            1605.243, 603.261, 26112.239, 2228.629, 7207.817,
-            156271.372, 661.114, 158184.4,
+            860.0, 200.0, 25610.0, 1000.0, 4685.0,
+            507500.0, 500.0, 150500.0,
         ])
-        # 
         self.assertIsInstance(dist_coast,np.ndarray)
-        np.testing.assert_array_almost_equal(dist_coast, REF_VALUES, decimal=3)
+        np.testing.assert_allclose(dist_coast, REF_VALUES, atol=1.0)
 
     def test_dist_coast_pass_raster(self):
-        """Test set_region_id"""
+        """Test get_dist_coast for centroids derived from a raster file"""
         centr_ras = Centroids.from_raster_file(HAZ_DEMO_FL, window=Window(0, 0, 50, 60))
         dist_coast = centr_ras.get_dist_coast()
         self.assertLess(abs(dist_coast[0] - 117000), 1000)
@@ -978,7 +977,7 @@ def test_area_pass_raster(self):
             981037.92674582, 981065.50487659, 981065.50487385,
         ])
         np.testing.assert_allclose(area_pixel, test_area)
-            
+
     def test_equal_pass(self):
         """Test equal"""
         centr_list = [

climada/hazard/test/test_tc_tracks.py

@@ -788,7 +788,8 @@ def test_dist_since_lf_pass(self):
 
         self.assertGreater(
             track.dist_since_lf.values[-1],
-            u_coord.dist_to_coast(track.lat.values[-1], track.lon.values[-1]) / 1000)
+            u_coord.dist_to_coast_nasa(track.lat.values[-1], track.lon.values[-1]) / 1000,
+        )
         self.assertEqual(1020.5431562223974, track['dist_since_lf'].values[-1])
 
         # check distances on land always increase, in second landfall

climada/util/coordinates.py

@@ -457,36 +457,6 @@ def grid_is_regular(coord):
         regular = True
     return regular, count_lat[0], count_lon[0]
 
-def get_coastlines(bounds=None, resolution=110):
-    """Get Polygones of coast intersecting given bounds
-
-    Parameters
-    ----------
-    bounds : tuple
-        min_lon, min_lat, max_lon, max_lat in EPSG:4326
-    resolution : float, optional
-        10, 50 or 110. Resolution in m. Default: 110m, i.e. 1:110.000.000
-
-    Returns
-    -------
-    coastlines : GeoDataFrame
-        Polygons of coast intersecting given bounds.
-    """
-    resolution = nat_earth_resolution(resolution)
-    shp_file = shapereader.natural_earth(resolution=resolution,
-                                         category='physical',
-                                         name='coastline')
-    coast_df = gpd.read_file(shp_file)
-    coast_df.crs = NE_CRS
-    if bounds is None:
-        return coast_df[['geometry']]
-    tot_coast = np.zeros(1)
-    while not np.any(tot_coast):
-        tot_coast = coast_df.envelope.intersects(box(*bounds))
-        bounds = (bounds[0] - 20, bounds[1] - 20,
-                  bounds[2] + 20, bounds[3] + 20)
-    return coast_df[tot_coast][['geometry']]
-
 def convert_wgs_to_utm(lon, lat):
     """Get EPSG code of UTM projection for input point in EPSG 4326
 
@@ -505,49 +475,30 @@ def convert_wgs_to_utm(lon, lat):
     epsg_utm_base = 32601 + (0 if lat >= 0 else 100)
     return epsg_utm_base + (math.floor((lon + 180) / 6) % 60)
 
-def utm_zones(wgs_bounds):
-    """Get EPSG code and bounds of UTM zones covering specified region
-
-    Parameters
-    ----------
-    wgs_bounds : tuple
-        lon_min, lat_min, lon_max, lat_max
+def dist_to_coast(coord_lat, lon=None, highres=False, signed=False):
+    """Read interpolated (signed) distance to coast (in m) from NASA data
 
-    Returns
-    -------
-    zones : list of pairs (zone_epsg, zone_wgs_bounds)
-        EPSG code and bounding box in WGS coordinates.
-    """
-    lon_min, lat_min, lon_max, lat_max = wgs_bounds
-    lon_min, lon_max = max(-179.99, lon_min), min(179.99, lon_max)
-    utm_min, utm_max = [math.floor((l + 180) / 6) for l in [lon_min, lon_max]]
-    zones = []
-    for utm in range(utm_min, utm_max + 1):
-        epsg = 32601 + utm
-        bounds = (-180 + 6 * utm, 0, -180 + 6 * (utm + 1), 90)
-        if lat_max >= 0:
-            zones.append((epsg, bounds))
-        if lat_min < 0:
-            bounds = (bounds[0], -90, bounds[2], 0)
-            zones.append((epsg + 100, bounds))
-    return zones
-
-def dist_to_coast(coord_lat, lon=None, signed=False):
-    """Compute (signed) distance to coast from input points in meters.
+    Note: The NASA raster file is 300 MB and will be downloaded on first run!
 
     Parameters
     ----------
-    coord_lat : GeoDataFrame or np.array or float
-        One of the following:
-            * GeoDataFrame with geometry column in epsg:4326
-            * np.array with two columns, first for latitude of each point
-              and second with longitude in epsg:4326
-            * np.array with one dimension containing latitudes in epsg:4326
-            * float with a latitude value in epsg:4326
-    lon : np.array or float, optional
+    coord_lat : GeoDataFrame or np.ndarray or float
         One of the following:
-            * np.array with one dimension containing longitudes in epsg:4326
-            * float with a longitude value in epsg:4326
+
+        * GeoDataFrame with geometry column in epsg:4326
+        * np.array with two columns, first for latitude of each point
+          and second with longitude in epsg:4326
+        * np.array with one dimension containing latitudes in epsg:4326
+        * float with a latitude value in epsg:4326
+
+    lon : np.ndarray or float, optional
+        If given, one of the following:
+
+        * np.array with one dimension containing longitudes in epsg:4326
+        * float with a longitude value in epsg:4326
+
+    highres : bool, optional
+        Use full resolution of NASA data (much slower). Default: False.
     signed : bool
         If True, distance is signed with positive values off shore and negative values on land.
         Default: False
@@ -557,53 +508,21 @@ def dist_to_coast(coord_lat, lon=None, signed=False):
     dist : np.array
         (Signed) distance to coast in meters.
     """
-    if isinstance(coord_lat, (gpd.GeoDataFrame, gpd.GeoSeries)):
-        if not equal_crs(coord_lat.crs, NE_CRS):
-            raise ValueError('Input CRS is not %s' % str(NE_CRS))
-        geom = coord_lat
-    else:
-        if lon is None:
-            if isinstance(coord_lat, np.ndarray) and coord_lat.shape[1] == 2:
-                lat, lon = coord_lat[:, 0], coord_lat[:, 1]
-            else:
-                raise ValueError('Missing longitude values.')
+    if lon is None:
+        if isinstance(coord_lat, (gpd.GeoDataFrame, gpd.GeoSeries)):
+            if not equal_crs(coord_lat.crs, DEF_CRS):
+                raise ValueError('Input CRS is not %s' % str(DEF_CRS))
+            geom = coord_lat if isinstance(coord_lat, gpd.GeoSeries) else coord_lat["geometry"]
+            lon, lat = geom.x.values, geom.y.values
+        elif isinstance(coord_lat, np.ndarray) and coord_lat.shape[1] == 2:
+            lat, lon = coord_lat[:, 0], coord_lat[:, 1]
         else:
-            lat, lon = [np.asarray(v).reshape(-1) for v in [coord_lat, lon]]
-            if lat.size != lon.size:
-                raise ValueError('Mismatching input coordinates size: %s != %s'
-                                 % (lat.size, lon.size))
-        geom = gpd.GeoDataFrame(geometry=gpd.points_from_xy(lon, lat), crs=NE_CRS)
-
-    pad = 20
-    bounds = (geom.total_bounds[0] - pad, geom.total_bounds[1] - pad,
-              geom.total_bounds[2] + pad, geom.total_bounds[3] + pad)
-    coast = get_coastlines(bounds, 10).geometry
-    coast = gpd.GeoDataFrame(geometry=coast, crs=NE_CRS)
-    dist = np.empty(geom.shape[0])
-    zones = utm_zones(geom.geometry.total_bounds)
-    for izone, (epsg, bounds) in enumerate(zones):
-        to_crs = f"epsg:{epsg}"
-        zone_mask = (
-            (bounds[1] <= geom.geometry.y)
-            & (geom.geometry.y <= bounds[3])
-            & (bounds[0] <= geom.geometry.x)
-            & (geom.geometry.x <= bounds[2])
-        )
-        if np.count_nonzero(zone_mask) == 0:
-            continue
-        LOGGER.info("dist_to_coast: UTM %d (%d/%d)",
-                    epsg, izone + 1, len(zones))
-        bounds = geom[zone_mask].total_bounds
-        bounds = (bounds[0] - pad, bounds[1] - pad,
-                  bounds[2] + pad, bounds[3] + pad)
-        coast_mask = coast.envelope.intersects(box(*bounds))
-        utm_coast = coast[coast_mask].geometry.unary_union
-        utm_coast = gpd.GeoDataFrame(geometry=[utm_coast], crs=NE_CRS)
-        utm_coast = utm_coast.to_crs(to_crs).geometry[0]
-        dist[zone_mask] = geom[zone_mask].to_crs(to_crs).distance(utm_coast)
-    if signed:
-        dist[coord_on_land(geom.geometry.y, geom.geometry.x)] *= -1
-    return dist
+            raise ValueError('Missing longitude values.')
+    else:
+        lat, lon = [np.asarray(v).reshape(-1) for v in [coord_lat, lon]]
+        if lat.size != lon.size:
+            raise ValueError(f'Mismatching input coordinates size: {lat.size} != {lon.size}')
+    return dist_to_coast_nasa(lat, lon, highres=highres, signed=signed)
 
 def _get_dist_to_coast_nasa_tif():
     """Get the path to the NASA raster file for distance to coast.

climada/util/test/test_coordinates.py

@@ -969,27 +969,6 @@ def test_nat_earth_resolution_fail(self):
         with self.assertRaises(ValueError):
             u_coord.nat_earth_resolution(111)
 
-    def test_get_coastlines_all_pass(self):
-        """Check get_coastlines function over whole earth"""
-        coast = u_coord.get_coastlines(resolution=110)
-        tot_bounds = coast.total_bounds
-        self.assertEqual((134, 1), coast.shape)
-        self.assertAlmostEqual(tot_bounds[0], -180)
-        self.assertAlmostEqual(tot_bounds[1], -85.60903777)
-        self.assertAlmostEqual(tot_bounds[2], 180.00000044)
-        self.assertAlmostEqual(tot_bounds[3], 83.64513)
-
-    def test_get_coastlines_pass(self):
-        """Check get_coastlines function in defined extent"""
-        bounds = (-100, -55, -20, 35)
-        coast = u_coord.get_coastlines(bounds, resolution=110)
-        ex_box = box(bounds[0], bounds[1], bounds[2], bounds[3])
-        self.assertEqual(coast.shape[0], 14)
-        self.assertTrue(coast.total_bounds[2] < 0)
-        for _row, line in coast.iterrows():
-            if not ex_box.intersects(line.geometry):
-                self.assertEqual(1, 0)
-
     def test_get_land_geometry_country_pass(self):
         """get_land_geometry with selected countries."""
         iso_countries = ['DEU', 'VNM']
@@ -1041,26 +1020,6 @@ def test_on_land_pass(self):
         self.assertEqual(res.size, lat.size)
         np.testing.assert_array_equal(res[:3], [True, False, True])
 
-    def test_dist_to_coast(self):
-        """Test point in coast and point not in coast"""
-        points = np.array([
-            # Caribbean Sea:
-            [13.208333333333329, -59.625000000000014],
-            # South America:
-            [-12.497529, -58.849505],
-            # Very close to coast of Somalia:
-            [1.96768, 45.23219],
-        ])
-        dists = [2594.2071059573445, 1382985.2459744606, 0.088222234]
-        for d, p in zip(dists, points):
-            res = u_coord.dist_to_coast(*p)
-            self.assertAlmostEqual(d, res[0])
-
-        # All at once requires more than one UTM
-        res = u_coord.dist_to_coast(points)
-        for d, r in zip(dists, res):
-            self.assertAlmostEqual(d, r)
-
     def test_dist_to_coast_nasa(self):
         """Test point in coast and point not in coast"""
         points = np.array([