move grid area function to util

Author: NicolasColombi

climada/hazard/tc_tracks.py

@@ -48,7 +48,6 @@
 from matplotlib.collections import LineCollection
 from matplotlib.colors import BoundaryNorm, ListedColormap
 from matplotlib.lines import Line2D
-from scipy.sparse import csr_matrix
 from shapely.geometry import LineString, MultiLineString, Point
 from sklearn.metrics import DistanceMetric
 from tqdm import tqdm
@@ -2957,74 +2956,7 @@ def compute_track_density(
         hist_count += hist_new
 
     if density:
-        grid_area, _ = compute_grid_cell_area(res=res)
+        grid_area, _ = u_coord.compute_grid_cell_area(res=res)
         hist_count = hist_count / grid_area
 
     return hist_count, lat_bins, lon_bins
-
-
-def compute_grid_cell_area(res) -> tuple[np.ndarray]:
-    """
-    This function computes the area of each grid cell on a sphere (Earth), using latitude and
-    longitude bins based on the given resolution. The area is computed using the spherical cap
-    approximation for each grid cell. The function return a 2D matrix with the corresponding area.
-
-    The formula used to compute the area of each grid cell is derived from the integral of the
-    surface area of a spherical cap between squared latitude and longitude bins:
-
-    A = R**2 * (Δλ) * (sin(ϕ2) - sin(ϕ1))
-
-    Where:
-    - R is the radius of the Earth (in km).
-    - Δλ is the width of the grid cell in terms of longitude (in radians).
-    - sin(ϕ2) - sin(ϕ1) is the difference in the sine of the latitudes, which
-    accounts for the varying horizontal distance between longitudinal lines at different latitudes.
-
-    This formula is the direct integration over λ and ϕ2 of:
-
-    A = R**2 * Δλ * Δϕ * cos(ϕ1)
-
-    which approximate the grid cell area as a square computed by the multiplication of two
-    arc legths, with the logitudal arc length ajdusted by latitude.
-
-    Parameters:
-    ----------
-    res: int
-        The resolution of the grid in degrees. The grid will have cells of size `res x res` in
-        terms of latitude and longitude.
-
-    Returns:
-    -------
-    grid_area: np.ndarray
-        A 2D array of shape `(len(lat_bins)-1, len(lon_bins)-1)` containing the area of each grid
-        cell, in square kilometers. Each entry in the array corresponds to the area of the
-        corresponding grid cell on Earth.
-
-    Example:
-    --------
-    >>> res = 10  #10° resolution
-    >>> area = compute_grid_cell_area(res = res)
-    >>> print(area.shape) # (180/10, 360/10)
-    (18, 36)
-    """
-
-    lat_bins = np.linspace(-90, 90, int(180 / res))  # lat bins
-    lon_bins = np.linspace(-180, 180, int(360 / res))
-
-    R = 6371  # Earth's radius [km]
-    # Convert to radians
-    lat_bin_edges = np.radians(lat_bins)
-    lon_res_rad = np.radians(res)
-
-    # Compute area
-    areas = (
-        R**2
-        * lon_res_rad
-        * np.abs(np.sin(lat_bin_edges[1:]) - np.sin(lat_bin_edges[:-1]))
-    )
-    # Broadcast to create a full 2D grid of areas
-    grid_area = np.tile(
-        areas[:, np.newaxis], (1, len(lon_bins) - 1)
-    )  # each row as same area
-
-    return grid_area, [lat_bins, lon_bins]

climada/util/coordinates.py

@@ -459,6 +459,73 @@ def get_gridcellarea(lat, resolution=0.5, unit="ha"):
     return area
 
 
+def compute_grid_cell_area(res: float) -> tuple[np.ndarray]:
+    """
+    This function computes the area of each grid cell on a sphere (Earth), using latitude and
+    longitude bins based on the given resolution. The area is computed using the spherical cap
+    approximation for each grid cell. The function return a 2D matrix with the corresponding area.
+
+    The formula used to compute the area of each grid cell is derived from the integral of the
+    surface area of a spherical cap between squared latitude and longitude bins:
+
+    A = R**2 * (Δλ) * (sin(ϕ2) - sin(ϕ1))
+
+    Where:
+    - R is the radius of the Earth (in km).
+    - Δλ is the width of the grid cell in terms of longitude (in radians).
+    - sin(ϕ2) - sin(ϕ1) is the difference in the sine of the latitudes, which
+    accounts for the varying horizontal distance between longitudinal lines at different latitudes.
+
+    This formula is the direct integration over λ and ϕ2 of:
+
+    A = R**2 * Δλ * Δϕ * cos(ϕ1)
+
+    which approximate the grid cell area as a square computed by the multiplication of two
+    arc legths, with the logitudal arc length ajdusted by latitude.
+
+    Parameters:
+    ----------
+    res: int
+        The resolution of the grid in degrees. The grid will have cells of size `res x res` in
+        terms of latitude and longitude.
+
+    Returns:
+    -------
+    grid_area: np.ndarray
+        A 2D array of shape `(len(lat_bins)-1, len(lon_bins)-1)` containing the area of each grid
+        cell, in square kilometers. Each entry in the array corresponds to the area of the
+        corresponding grid cell on Earth.
+
+    Example:
+    --------
+    >>> res = 10  #10° resolution
+    >>> area = compute_grid_cell_area(res = res)
+    >>> print(area.shape) # (180/10, 360/10)
+    (18, 36)
+    """
+
+    lat_bins = np.linspace(-90, 90, int(180 / res))  # lat bins
+    lon_bins = np.linspace(-180, 180, int(360 / res))
+
+    R = 6371  # Earth's radius [km]
+    # Convert to radians
+    lat_bin_edges = np.radians(lat_bins)
+    lon_res_rad = np.radians(res)
+
+    # Compute area
+    areas = (
+        R**2
+        * lon_res_rad
+        * np.abs(np.sin(lat_bin_edges[1:]) - np.sin(lat_bin_edges[:-1]))
+    )
+    # Broadcast to create a full 2D grid of areas
+    grid_area = np.tile(
+        areas[:, np.newaxis], (1, len(lon_bins) - 1)
+    )  # each row as same area
+
+    return grid_area, [lat_bins, lon_bins]
+
+
 def grid_is_regular(coord):
     """Return True if grid is regular. If True, returns height and width.