Use geodesic conversion of meters, avoid empty polygon

fmtm_splitter/splitter.py

@@ -20,18 +20,19 @@
 import argparse
 import json
 import logging
+import math
 import sys
 from io import BytesIO
 from pathlib import Path
 from textwrap import dedent
-from typing import Optional, Union
+from typing import Optional, Tuple, Union
 
 import geojson
 import numpy as np
 from geojson import Feature, FeatureCollection, GeoJSON
 from osm_rawdata.postgres import PostgresClient
 from psycopg2.extensions import connection
-from shapely.geometry import Polygon, shape
+from shapely.geometry import Polygon, box, shape
 from shapely.geometry.geo import mapping
 from shapely.ops import unary_union
 
@@ -151,6 +152,46 @@ def geojson_to_shapely_polygon(
 
         return shape(features[0].get("geometry"))
 
+    def meters_to_degrees(
+        self, meters: float, reference_lat: float
+    ) -> Tuple[float, float]:
+        """Converts meters to degrees at a given latitude using WGS84 ellipsoidal calculations.
+
+        Args:
+            meters (float): The distance in meters to convert.
+            reference_lat (float): The latitude at which to perform the conversion (in degrees).
+
+        Returns:
+            Tuple[float, float]: Degree values for latitude and longitude.
+        """
+        # INFO:
+        # The geodesic distance is the shortest distance on the surface
+        # of an ellipsoidal model of the earth
+
+        lat_rad = math.radians(reference_lat)
+
+        # Using WGS84 parameters
+        a = 6378137.0  # Semi-major axis in meters
+        f = 1 / 298.257223563  # Flattening factor
+
+        # Applying formula
+        e2 = (2 * f) - (f**2)  # Eccentricity squared
+        N = a / math.sqrt(
+            1 - e2 * math.sin(lat_rad) ** 2
+        )  # Radius of curvature in the prime vertical
+        M = (
+            a * (1 - e2) / (1 - e2 * math.sin(lat_rad) ** 2) ** (3 / 2)
+        )  # Radius of curvature in the meridian
+
+        lat_deg_change = meters / M  # Latitude change in degrees
+        lon_deg_change = meters / (N * math.cos(lat_rad))  # Longitude change in degrees
+
+        # Convert changes to degrees by dividing by radians to degrees
+        lat_deg_change = math.degrees(lat_deg_change)
+        lon_deg_change = math.degrees(lon_deg_change)
+
+        return lat_deg_change, lon_deg_change
+
     def splitBySquare(  # noqa: N802
         self,
         meters: int,
@@ -170,14 +211,13 @@ def splitBySquare(  # noqa: N802
 
         xmin, ymin, xmax, ymax = self.aoi.bounds
 
-        # 1 meters is this factor in degrees
-        meter = 0.0000114
-        length = float(meters) * meter
-        width = float(meters) * meter
+        reference_lat = (ymin + ymax) / 2
+        length_deg, width_deg = self.meters_to_degrees(meters, reference_lat)
+
+        # Create grid columns and rows based on the AOI bounds
+        cols = np.arange(xmin, xmax + width_deg, width_deg)
+        rows = np.arange(ymin, ymax + length_deg, length_deg)
 
-        cols = list(np.arange(xmin, xmax + width, width))
-        rows = list(np.arange(ymin, ymax + length, length))
-        polygons = []
         extract_geoms = []
         if extract_geojson:
             features = (
@@ -187,14 +227,18 @@ def splitBySquare(  # noqa: N802
             )
             extract_geoms = [shape(feature["geometry"]) for feature in features]
 
+        # Generate grid polygons and clip them by AOI
+        polygons = []
         for x in cols[:-1]:
             for y in rows[:-1]:
-                grid_polygon = Polygon(
-                    [(x, y), (x + width, y), (x + width, y + length), (x, y + length)]
-                )
+                grid_polygon = box(x, y, x + width_deg, y + length_deg)
                 clipped_polygon = grid_polygon.intersection(self.aoi)
+
+                if clipped_polygon.is_empty:
+                    continue
+
+                # Check intersection with extract geometries if available
                 if extract_geoms:
-                    # Check if any extract geometry is within the clipped grid
                     if any(geom.within(clipped_polygon) for geom in extract_geoms):
                         polygons.append(clipped_polygon)
                 else: