diff --git a/fmtm_splitter/splitter.py b/fmtm_splitter/splitter.py index c2497a8..dcea438 100755 --- a/fmtm_splitter/splitter.py +++ b/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( diff --git a/tests/test_splitter.py b/tests/test_splitter.py index 589d87d..cc88db5 100644 --- a/tests/test_splitter.py +++ b/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():