Merge pull request #60 from hotosm/fix/square-split

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,49 @@ 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 +214,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 +230,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:
@@ -705,6 +752,7 @@ def main(args_list: list[str] | None = None):
         "--meters",
         nargs="?",
         const=50,
+        type=int,
         help="Size in meters if using square splitting",
     )
     parser.add_argument(

tests/test_splitter.py

@@ -66,11 +66,11 @@ def test_split_by_square_with_dict(aoi_json, extract_json):
     features = split_by_square(
         aoi_json.get("features")[0], meters=50, osm_extract=extract_json
     )
-    assert len(features.get("features")) == 50
+    assert len(features.get("features")) == 60
     features = split_by_square(
         aoi_json.get("features")[0].get("geometry"), meters=50, osm_extract=extract_json
     )
-    assert len(features.get("features")) == 50
+    assert len(features.get("features")) == 60
 
 
 def test_split_by_square_with_str(aoi_json, extract_json):
@@ -79,21 +79,21 @@ def test_split_by_square_with_str(aoi_json, extract_json):
     features = split_by_square(
         geojson.dumps(aoi_json.get("features")[0]), meters=50, osm_extract=extract_json
     )
-    assert len(features.get("features")) == 50
+    assert len(features.get("features")) == 60
     # JSON Dumps
     features = split_by_square(
         json.dumps(aoi_json.get("features")[0].get("geometry")),
         meters=50,
         osm_extract=extract_json,
     )
-    assert len(features.get("features")) == 50
+    assert len(features.get("features")) == 60
     # File
     features = split_by_square(
         "tests/testdata/kathmandu.geojson",
         meters=100,
         osm_extract="tests/testdata/kathmandu_extract.geojson",
     )
-    assert len(features.get("features")) == 15
+    assert len(features.get("features")) == 20
 
 
 def test_split_by_square_with_file_output():
@@ -108,11 +108,11 @@ def test_split_by_square_with_file_output():
         meters=50,
         outfile=str(outfile),
     )
-    assert len(features.get("features")) == 50
+    assert len(features.get("features")) == 60
     # Also check output file
     with open(outfile, "r") as jsonfile:
         output_geojson = geojson.load(jsonfile)
-    assert len(output_geojson.get("features")) == 50
+    assert len(output_geojson.get("features")) == 60
 
 
 def test_split_by_square_with_multigeom_input(aoi_multi_json, extract_json):
@@ -125,7 +125,7 @@ def test_split_by_square_with_multigeom_input(aoi_multi_json, extract_json):
         osm_extract=extract_json,
         outfile=str(outfile),
     )
-    assert len(features.get("features", [])) == 50
+    assert len(features.get("features", [])) == 80
     for index in [0, 1, 2, 3]:
         assert Path(f"{Path(outfile).stem}_{index}.geojson)").exists()
 
@@ -231,7 +231,7 @@ def test_split_by_square_cli():
     with open(outfile, "r") as jsonfile:
         output_geojson = geojson.load(jsonfile)
 
-    assert len(output_geojson.get("features")) == 15
+    assert len(output_geojson.get("features")) == 20
 
 
 def test_split_by_features_cli():