From 6f0b252181f4c4b23c80703a577c6705736bc6b7 Mon Sep 17 00:00:00 2001 From: Sujan Adhikari Date: Mon, 28 Oct 2024 14:07:52 +0545 Subject: [PATCH 1/6] feat: use geodetic conversion of meters to degrees --- fmtm_splitter/splitter.py | 60 ++++++++++++++++++++++++++++++++------- 1 file changed, 49 insertions(+), 11 deletions(-) diff --git a/fmtm_splitter/splitter.py b/fmtm_splitter/splitter.py index c2497a8..31071ad 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, shape, box from shapely.geometry.geo import mapping from shapely.ops import unary_union @@ -151,6 +152,45 @@ 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,13 +210,12 @@ 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(reference_lat, meters) - cols = list(np.arange(xmin, xmax + width, width)) - rows = list(np.arange(ymin, ymax + length, length)) + # 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) polygons = [] extract_geoms = [] if extract_geojson: @@ -189,10 +228,9 @@ def splitBySquare( # noqa: N802 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 extract_geoms: # Check if any extract geometry is within the clipped grid if any(geom.within(clipped_polygon) for geom in extract_geoms): From 1bc6610ac6e4084afc905b9bca5c8de320f0afbb Mon Sep 17 00:00:00 2001 From: Sujan Adhikari Date: Mon, 28 Oct 2024 14:10:00 +0545 Subject: [PATCH 2/6] fix: avoid appending empty clipped polygons --- fmtm_splitter/splitter.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/fmtm_splitter/splitter.py b/fmtm_splitter/splitter.py index 31071ad..5b07bdb 100755 --- a/fmtm_splitter/splitter.py +++ b/fmtm_splitter/splitter.py @@ -211,12 +211,12 @@ def splitBySquare( # noqa: N802 xmin, ymin, xmax, ymax = self.aoi.bounds reference_lat = (ymin + ymax) / 2 - length_deg, width_deg = self.meters_to_degrees(reference_lat, meters) + 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) - polygons = [] + extract_geoms = [] if extract_geojson: features = ( @@ -226,13 +226,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 = 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: From 0537537dc7394dd53e42f864c42c7ba585e28bbc Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 28 Oct 2024 08:42:30 +0000 Subject: [PATCH 3/6] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- fmtm_splitter/splitter.py | 33 +++++++++++++++++---------------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/fmtm_splitter/splitter.py b/fmtm_splitter/splitter.py index 5b07bdb..e765aa1 100755 --- a/fmtm_splitter/splitter.py +++ b/fmtm_splitter/splitter.py @@ -32,7 +32,7 @@ from geojson import Feature, FeatureCollection, GeoJSON from osm_rawdata.postgres import PostgresClient from psycopg2.extensions import connection -from shapely.geometry import Polygon, shape, box +from shapely.geometry import Polygon, box, shape from shapely.geometry.geo import mapping from shapely.ops import unary_union @@ -153,17 +153,14 @@ 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. - + 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. """ @@ -174,21 +171,25 @@ def meters_to_degrees( lat_rad = math.radians(reference_lat) # Using WGS84 parameters - a = 6378137.0 # Semi-major axis in meters + 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 - + 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 From 2bdd38f902769153ab1afa0ea2404525a791e284 Mon Sep 17 00:00:00 2001 From: Sujan Adhikari Date: Tue, 29 Oct 2024 15:30:49 +0545 Subject: [PATCH 4/6] fix: add type int in argument --- fmtm_splitter/splitter.py | 1 + 1 file changed, 1 insertion(+) diff --git a/fmtm_splitter/splitter.py b/fmtm_splitter/splitter.py index 5b07bdb..331ce6f 100755 --- a/fmtm_splitter/splitter.py +++ b/fmtm_splitter/splitter.py @@ -748,6 +748,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( From bd5bdb787827b618aca6814401a9ee36d0847b4e Mon Sep 17 00:00:00 2001 From: Sujan Adhikari Date: Tue, 29 Oct 2024 15:31:34 +0545 Subject: [PATCH 5/6] fix: updated the tasks count in test cases of split by square --- tests/test_splitter.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) 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(): From 48d0e18dab65cfbbde93722fb5d459f548cefd72 Mon Sep 17 00:00:00 2001 From: Sujan Adhikari Date: Tue, 29 Oct 2024 15:45:26 +0545 Subject: [PATCH 6/6] fix: precommit errors --- fmtm_splitter/splitter.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/fmtm_splitter/splitter.py b/fmtm_splitter/splitter.py index aa800e8..dcea438 100755 --- a/fmtm_splitter/splitter.py +++ b/fmtm_splitter/splitter.py @@ -155,11 +155,14 @@ def geojson_to_shapely_polygon( 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. + """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). + reference_lat (float): The latitude at which to , + perform the conversion (in degrees). Returns: Tuple[float, float]: Degree values for latitude and longitude. @@ -176,15 +179,15 @@ def meters_to_degrees( # Applying formula e2 = (2 * f) - (f**2) # Eccentricity squared - N = a / math.sqrt( + n = a / math.sqrt( 1 - e2 * math.sin(lat_rad) ** 2 ) # Radius of curvature in the prime vertical - M = ( + 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 + 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)