diff --git a/pyproject.toml b/pyproject.toml index 7cd91ae..99ab3a8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -16,9 +16,11 @@ requires-python = ">=3.8" keywords = ["WIS2.0", "CAP", "XML", "GeoJSON", "convert"] license = {file = "LICENSE"} dependencies = [ - "xmltodict>=0.13.0", + "click>=8.1.7", "geojson>=3.1.0", - "click>=8.1.7" + "pyproj>=3.6.0", + "shapely>=2.0.0", + "xmltodict>=0.13.0", ] dynamic = ["version"] diff --git a/src/cap2geojson/__init__.py b/src/cap2geojson/__init__.py index c2bfd33..115245f 100644 --- a/src/cap2geojson/__init__.py +++ b/src/cap2geojson/__init__.py @@ -23,7 +23,7 @@ from .convert import to_geojson -__version__ = "0.1.0-dev3" +__version__ = "0.1.0-dev4" logging.basicConfig( level=logging.INFO, diff --git a/src/cap2geojson/convert.py b/src/cap2geojson/convert.py index a9af932..9efda64 100644 --- a/src/cap2geojson/convert.py +++ b/src/cap2geojson/convert.py @@ -1,10 +1,14 @@ -import geojson import json import logging -import math -from typing import Generator, Union import re +from typing import Generator, Union + +import geojson import xmltodict +from pyproj import Transformer +from shapely.geometry import Point +from shapely.geometry.polygon import orient +from shapely.ops import transform as transform_bufr logger = logging.getLogger(__name__) @@ -61,31 +65,56 @@ def get_area_desc(area: Union[dict, list]) -> str: def get_circle_coords( - x_centre: float, y_centre: float, radius: float, n_points: int + x_centre: float, y_centre: float, radius: float ) -> Generator[list, None, None]: """ - Estimate the n+1 coordinates of a circle with a given centre and radius. + Estimate the n+1 coordinates of a circle with a given centre and radius + using the azimuthal equidistant projection. Args: x_centre (float): The longitude of the circle's centre. y_centre (float): The latitude of the circle's centre. - radius (float): The radius of the circle. - n_points (int): The number of edges in the n-gon to approximate - the circle. + radius (float): The radius of the circle in kilometres. Returns: Generator: Yield the n+1 estimated coordinates of the circle. """ - # Generate thetas for the n-gon - thetas = [i / n_points * math.tau for i in range(n_points)] - # Add theta = 0 to the end to ensure the circle is closed - thetas.append(0) - - for theta in thetas: - x = radius * math.cos(theta) + x_centre - y = radius * math.sin(theta) + y_centre - # Round to 5 decimal places to prevent excessive precision - yield [round(x, 5), round(y, 5)] + # Validate latitude and longitude + if not (-90 <= y_centre <= 90): + raise ValueError( + f"Invalid latitude value: {y_centre}. Must be between -90 and 90 degrees." # noqa + ) + if not (-180 <= x_centre <= 180): + raise ValueError( + f"Invalid longitude value: {x_centre}. Must be between -180 and 180 degrees." # noqa + ) + + # Create local azimuthal equidistant projection + local_azimuthal_projection = ( + f"+proj=aeqd +R=6371000 +units=m +lat_0={y_centre} +lon_0={x_centre}" + ) + wgs84_to_aeqd = Transformer.from_proj( + "+proj=longlat +datum=WGS84 +no_defs", local_azimuthal_projection + ) + aeqd_to_wgs84 = Transformer.from_proj( + local_azimuthal_projection, "+proj=longlat +datum=WGS84 +no_defs" + ) + + # Transform the center point to the local projection + point_transformed = Point(wgs84_to_aeqd.transform(x_centre, y_centre)) + + # Create a buffer around the transformed point, with radius in metres + buffer = point_transformed.buffer(radius*1000) + + # Transform the buffer back to WGS84 coordinates + circle = transform_bufr(aeqd_to_wgs84.transform, buffer) + + # Ensure the coordinates follow the right-hand rule (counter-clockwise) + circle = orient(circle, sign=1.0) + + # Extract the coordinates from the transformed buffer + for coord in circle.exterior.coords: + yield [round(coord[0], 5), round(coord[1], 5)] def ensure_counter_clockwise(coords: list) -> list: @@ -132,7 +161,7 @@ def get_polygon_coordinates(single_area: dict) -> list: radius = float(radius) y_centre, x_centre = map(float, centre.split(",")) # Estimate the circle coordinates with n=100 points - return list(get_circle_coords(x_centre, y_centre, radius, 100)) + return list(get_circle_coords(x_centre, y_centre, radius)) if "polygon" in single_area: # Takes form "y,x y,x y,x". So, split on whitespace, then comma, and diff --git a/tests/output/circle_estimation.json b/tests/output/circle_estimation.json index 6b1f7da..714de7d 100644 --- a/tests/output/circle_estimation.json +++ b/tests/output/circle_estimation.json @@ -1,406 +1,262 @@ [ [ - 12.0, + 5.06304, 3.0 ], [ - 11.98619, - 3.43953 + 5.06274, + 3.00617 ], [ - 11.9448, - 3.87733 + 5.06183, + 3.01228 ], [ - 11.87601, - 4.31167 + 5.06033, + 3.01827 ], [ - 11.78008, - 4.74083 + 5.05824, + 3.02409 ], [ - 11.6574, - 5.16312 + 5.0556, + 3.02967 ], [ - 11.50844, - 5.57687 + 5.05242, + 3.03497 ], [ - 11.33379, - 5.98046 + 5.04873, + 3.03994 ], [ - 11.13415, - 6.37228 + 5.04458, + 3.04451 ], [ - 10.9103, - 6.75079 + 5.03999, + 3.04866 ], [ - 10.66312, - 7.1145 + 5.03502, + 3.05234 ], [ - 10.39359, - 7.46197 + 5.02972, + 3.05552 ], [ - 10.10278, - 7.79183 + 5.02413, + 3.05816 ], [ - 9.79183, - 8.10278 + 5.0183, + 3.06024 ], [ - 9.46197, - 8.39359 + 5.0123, + 3.06174 ], [ - 9.1145, - 8.66312 - ], - [ - 8.75079, - 8.9103 - ], - [ - 8.37228, - 9.13415 - ], - [ - 7.98046, - 9.33379 - ], - [ - 7.57687, - 9.50844 - ], - [ - 7.16312, - 9.6574 - ], - [ - 6.74083, - 9.78008 - ], - [ - 6.31167, - 9.87601 - ], - [ - 5.87733, - 9.9448 - ], - [ - 5.43953, - 9.98619 + 5.00618, + 3.06265 ], [ 5.0, - 10.0 - ], - [ - 4.56047, - 9.98619 - ], - [ - 4.12267, - 9.9448 - ], - [ - 3.68833, - 9.87601 - ], - [ - 3.25917, - 9.78008 - ], - [ - 2.83688, - 9.6574 - ], - [ - 2.42313, - 9.50844 - ], - [ - 2.01954, - 9.33379 - ], - [ - 1.62772, - 9.13415 - ], - [ - 1.24921, - 8.9103 + 3.06295 ], [ - 0.8855, - 8.66312 + 4.99382, + 3.06265 ], [ - 0.53803, - 8.39359 + 4.9877, + 3.06174 ], [ - 0.20817, - 8.10278 + 4.9817, + 3.06024 ], [ - -0.10278, - 7.79183 + 4.97587, + 3.05816 ], [ - -0.39359, - 7.46197 + 4.97028, + 3.05552 ], [ - -0.66312, - 7.1145 + 4.96498, + 3.05234 ], [ - -0.9103, - 6.75079 + 4.96001, + 3.04866 ], [ - -1.13415, - 6.37228 + 4.95542, + 3.04451 ], [ - -1.33379, - 5.98046 + 4.95127, + 3.03994 ], [ - -1.50844, - 5.57687 + 4.94758, + 3.03497 ], [ - -1.6574, - 5.16312 + 4.9444, + 3.02967 ], [ - -1.78008, - 4.74083 + 4.94176, + 3.02409 ], [ - -1.87601, - 4.31167 + 4.93967, + 3.01827 ], [ - -1.9448, - 3.87733 + 4.93817, + 3.01228 ], [ - -1.98619, - 3.43953 + 4.93726, + 3.00617 ], [ - -2.0, + 4.93696, 3.0 ], [ - -1.98619, - 2.56047 + 4.93726, + 2.99383 ], [ - -1.9448, - 2.12267 + 4.93817, + 2.98772 ], [ - -1.87601, - 1.68833 + 4.93968, + 2.98172 ], [ - -1.78008, - 1.25917 + 4.94176, + 2.97591 ], [ - -1.6574, - 0.83688 + 4.94441, + 2.97032 ], [ - -1.50844, - 0.42313 + 4.94759, + 2.96502 ], [ - -1.33379, - 0.01954 + 4.95127, + 2.96006 ], [ - -1.13415, - -0.37228 + 4.95543, + 2.95548 ], [ - -0.9103, - -0.75079 + 4.96001, + 2.95134 ], [ - -0.66312, - -1.1145 + 4.96498, + 2.94766 ], [ - -0.39359, - -1.46197 + 4.97029, + 2.94448 ], [ - -0.10278, - -1.79183 + 4.97588, + 2.94184 ], [ - 0.20817, - -2.10278 + 4.9817, + 2.93976 ], [ - 0.53803, - -2.39359 + 4.9877, + 2.93826 ], [ - 0.8855, - -2.66312 - ], - [ - 1.24921, - -2.9103 - ], - [ - 1.62772, - -3.13415 - ], - [ - 2.01954, - -3.33379 - ], - [ - 2.42313, - -3.50844 - ], - [ - 2.83688, - -3.6574 - ], - [ - 3.25917, - -3.78008 - ], - [ - 3.68833, - -3.87601 - ], - [ - 4.12267, - -3.9448 - ], - [ - 4.56047, - -3.98619 + 4.99382, + 2.93735 ], [ 5.0, - -4.0 - ], - [ - 5.43953, - -3.98619 - ], - [ - 5.87733, - -3.9448 - ], - [ - 6.31167, - -3.87601 - ], - [ - 6.74083, - -3.78008 - ], - [ - 7.16312, - -3.6574 - ], - [ - 7.57687, - -3.50844 - ], - [ - 7.98046, - -3.33379 - ], - [ - 8.37228, - -3.13415 - ], - [ - 8.75079, - -2.9103 + 2.93705 ], [ - 9.1145, - -2.66312 + 5.00618, + 2.93735 ], [ - 9.46197, - -2.39359 + 5.0123, + 2.93826 ], [ - 9.79183, - -2.10278 + 5.0183, + 2.93976 ], [ - 10.10278, - -1.79183 + 5.02412, + 2.94184 ], [ - 10.39359, - -1.46197 + 5.02971, + 2.94448 ], [ - 10.66312, - -1.1145 + 5.03502, + 2.94766 ], [ - 10.9103, - -0.75079 + 5.03999, + 2.95134 ], [ - 11.13415, - -0.37228 + 5.04457, + 2.95548 ], [ - 11.33379, - 0.01954 + 5.04873, + 2.96006 ], [ - 11.50844, - 0.42313 + 5.05241, + 2.96502 ], [ - 11.6574, - 0.83688 + 5.05559, + 2.97032 ], [ - 11.78008, - 1.25917 + 5.05824, + 2.97591 ], [ - 11.87601, - 1.68833 + 5.06032, + 2.98172 ], [ - 11.9448, - 2.12267 + 5.06183, + 2.98772 ], [ - 11.98619, - 2.56047 + 5.06274, + 2.99383 ], [ - 12.0, + 5.06304, 3.0 ] ] diff --git a/tests/test_cap2geojson.py b/tests/test_cap2geojson.py index d63c5ed..7466062 100644 --- a/tests/test_cap2geojson.py +++ b/tests/test_cap2geojson.py @@ -21,11 +21,11 @@ import json import logging -import pytest import time +import pytest + from cap2geojson.convert import ( - get_circle_coords, ensure_counter_clockwise, get_polygon_coordinates, preprocess_alert, @@ -53,33 +53,6 @@ def test_to_geojson(sc_alert): assert actual == expected -@pytest.fixture -def circle(): - return { - "x_centre": 5.0, - "y_centre": 3.0, - "radius": 7.0, - } - - -def test_circle_coords(circle): - assert list(get_circle_coords( - circle["x_centre"], circle["y_centre"], circle["radius"], 10 - )) == [ - [12.0, 3.0], - [10.66312, 7.1145], - [7.16312, 9.6574], - [2.83688, 9.6574], - [-0.66312, 7.1145], - [-2.0, 3.0], - [-0.66312, -1.1145], - [2.83688, -3.6574], - [7.16312, -3.6574], - [10.66312, -1.1145], - [12.0, 3.0], - ] - - @pytest.fixture def circle_area(): return {