From ff50a10ec19ebd5d94bd7fa44d227da702074b62 Mon Sep 17 00:00:00 2001
From: josidd <joseph.siddons@noc.ac.uk>
Date: Wed, 25 Sep 2024 13:46:28 +0100
Subject: [PATCH] feat: Comparison between quad/octtree and ellipse on surface
of earth
---
GeoSpatialTools/octtree.py | 126 +++++++++++++++++++++++++++++++++++-
GeoSpatialTools/quadtree.py | 102 ++++++++++++++++++++++++++++-
2 files changed, 225 insertions(+), 3 deletions(-)
diff --git a/GeoSpatialTools/octtree.py b/GeoSpatialTools/octtree.py
index 05f0904..28cee15 100644
--- a/GeoSpatialTools/octtree.py
+++ b/GeoSpatialTools/octtree.py
@@ -1,5 +1,7 @@
from datetime import datetime, timedelta
from .distance_metrics import haversine
+from .distance_metrics import haversine, destination
+from math import degrees, sqrt
class SpaceTimeRecord:
@@ -182,6 +184,11 @@ class SpaceTimeRectangle:
-------
bool : True if the point is <= dist + max(dist(centre, corners))
"""
+ if (
+ point.datetime - t_dist > self.datetime + self.dt / 2
+ or point.datetime + t_dist < self.datetime - self.dt / 2
+ ):
+ return False
# QUESTION: Is this sufficient? Possibly it is overkill
corner_dist = max(
haversine(
@@ -207,8 +214,96 @@ class SpaceTimeRectangle:
return (
haversine(self.lon, self.lat, point.lon, point.lat)
<= dist + corner_dist
- and point.datetime - t_dist <= self.datetime + self.dt / 2
- and point.datetime + t_dist >= self.datetime - self.dt / 2
+ )
+
+
+class SpaceTimeEllipse:
+ """A simple Ellipse Class for an ellipse on the surface of a sphere."""
+
+ def __init__(
+ self,
+ a: float,
+ b: float,
+ lon: float,
+ lat: float,
+ theta: float,
+ datetime: datetime,
+ dt: timedelta,
+ ) -> None:
+ self.a = a
+ self.b = b
+ self.lon = lon
+ self.lat = lat
+ self.datetime = datetime
+ self.dt = dt
+ # theta is anti-clockwise angle from horizontal in radians
+ self.theta = theta
+ # bearing is angle clockwise from north in degrees
+ self.bearing = (90 - degrees(self.theta)) % 360
+
+ a2 = self.a * self.a
+ b2 = self.b * self.b
+
+ self.c = sqrt(a2 - b2)
+ self.p1_lon, self.p1_lat = destination(
+ self.lon,
+ self.lat,
+ self.bearing,
+ self.c,
+ )
+ self.p2_lon, self.p2_lat = destination(
+ self.lon,
+ self.lat,
+ (180 - self.bearing) % 360,
+ self.c,
+ )
+
+ def contains(self, point: SpaceTimeRecord) -> bool:
+ """Test if a point is contained within the Ellipse"""
+ return (
+ (
+ haversine(self.p1_lon, self.p1_lat, point.lon, point.lat)
+ + haversine(self.p2_lon, self.p2_lat, point.lon, point.lat)
+ )
+ <= 2 * self.a
+ and point.datetime <= self.datetime + self.dt / 2
+ and point.datetime >= self.datetime - self.dt / 2
+ )
+
+ def nearby_rect(self, rect: SpaceTimeRectangle) -> bool:
+ """Test if a rectangle is near to the Ellipse"""
+ if (
+ rect.datetime - rect.dt / 2 > self.datetime + self.dt / 2
+ or rect.datetime + rect.dt / 2 < self.datetime - self.dt / 2
+ ):
+ return False
+ # TODO: Check corners, and 0 lat
+ corner_dist = max(
+ haversine(
+ rect.lon,
+ rect.lat,
+ rect.lon + rect.lon_range / 2,
+ rect.lat + rect.lat_range / 2,
+ ),
+ haversine(
+ rect.lon,
+ rect.lat,
+ rect.lon + rect.lon_range / 2,
+ rect.lat - rect.lat_range / 2,
+ ),
+ )
+ if (rect.lat + rect.lat_range / 2) * (
+ rect.lat - rect.lat_range / 2
+ ) < 0:
+ corner_dist = max(
+ corner_dist,
+ haversine(rect.lon, rect.lat, rect.lon + rect.lon_range / 2, 0),
+ )
+ return (
+ haversine(self.p1_lon, self.p1_lat, rect.lon, rect.lat)
+ <= corner_dist + self.a
+ or haversine(self.p2_lon, self.p2_lat, rect.lon, rect.lat)
+ <= corner_dist + self.a
)
@@ -459,6 +554,33 @@ class OctTree:
return points
+ def query_ellipse(
+ self,
+ ellipse: SpaceTimeEllipse,
+ points: SpaceTimeRecords | None = None,
+ ) -> SpaceTimeRecords:
+ """Get points that fall in an ellipse."""
+ if not points:
+ points = SpaceTimeRecords()
+ if not ellipse.nearby_rect(self.boundary):
+ return points
+
+ for point in self.points:
+ if ellipse.contains(point):
+ points.append(point)
+
+ if self.divided:
+ points = self.northwestfwd.query_ellipse(ellipse, points)
+ points = self.northeastfwd.query_ellipse(ellipse, points)
+ points = self.southwestfwd.query_ellipse(ellipse, points)
+ points = self.southeastfwd.query_ellipse(ellipse, points)
+ points = self.northwestback.query_ellipse(ellipse, points)
+ points = self.northeastback.query_ellipse(ellipse, points)
+ points = self.southwestback.query_ellipse(ellipse, points)
+ points = self.southeastback.query_ellipse(ellipse, points)
+
+ return points
+
def nearby_points(
self,
point: SpaceTimeRecord,
diff --git a/GeoSpatialTools/quadtree.py b/GeoSpatialTools/quadtree.py
index 474dac9..b5f2088 100644
--- a/GeoSpatialTools/quadtree.py
+++ b/GeoSpatialTools/quadtree.py
@@ -4,7 +4,8 @@ for detecting nearby records for example
"""
from datetime import datetime
-from .distance_metrics import haversine
+from .distance_metrics import haversine, destination
+from math import degrees, sqrt
class Record:
@@ -144,6 +145,82 @@ class Rectangle:
)
+class Ellipse:
+ """A simple Ellipse Class for an ellipse on the surface of a sphere."""
+
+ def __init__(
+ self,
+ a: float,
+ b: float,
+ lon: float,
+ lat: float,
+ theta: float,
+ ) -> None:
+ self.a = a
+ self.b = b
+ self.lon = lon
+ self.lat = lat
+ # theta is anti-clockwise angle from horizontal in radians
+ self.theta = theta
+ # bearing is angle clockwise from north in degrees
+ self.bearing = (90 - degrees(self.theta)) % 360
+
+ a2 = self.a * self.a
+ b2 = self.b * self.b
+
+ self.c = sqrt(a2 - b2)
+ self.p1_lon, self.p1_lat = destination(
+ self.lon,
+ self.lat,
+ self.bearing,
+ self.c,
+ )
+ self.p2_lon, self.p2_lat = destination(
+ self.lon,
+ self.lat,
+ (180 - self.bearing) % 360,
+ self.c,
+ )
+
+ def contains(self, point: Record) -> bool:
+ """Test if a point is contained within the Ellipse"""
+ return (
+ haversine(self.p1_lon, self.p1_lat, point.lon, point.lat)
+ + haversine(self.p2_lon, self.p2_lat, point.lon, point.lat)
+ ) <= 2 * self.a
+
+ def nearby_rect(self, rect: Rectangle) -> bool:
+ """Test if a rectangle is near to the Ellipse"""
+ # TODO: Check corners, and 0 lat
+ corner_dist = max(
+ haversine(
+ rect.lon,
+ rect.lat,
+ rect.lon + rect.lon_range / 2,
+ rect.lat + rect.lat_range / 2,
+ ),
+ haversine(
+ rect.lon,
+ rect.lat,
+ rect.lon + rect.lon_range / 2,
+ rect.lat - rect.lat_range / 2,
+ ),
+ )
+ if (rect.lat + rect.lat_range / 2) * (
+ rect.lat - rect.lat_range / 2
+ ) < 0:
+ corner_dist = max(
+ corner_dist,
+ haversine(rect.lon, rect.lat, rect.lon + rect.lon_range / 2, 0),
+ )
+ return (
+ haversine(self.p1_lon, self.p1_lat, rect.lon, rect.lat)
+ <= corner_dist + self.a
+ or haversine(self.p2_lon, self.p2_lat, rect.lon, rect.lat)
+ <= corner_dist + self.a
+ )
+
+
class QuadTree:
"""
A Simple QuadTree class for PyCOADS
@@ -288,6 +365,29 @@ class QuadTree:
return points
+ def query_ellipse(
+ self,
+ ellipse: Ellipse,
+ points: list[Record] | None = None,
+ ) -> list[Record]:
+ """Get points that fall in an ellipse."""
+ if not points:
+ points = list()
+ if not ellipse.nearby_rect(self.boundary):
+ return points
+
+ for point in self.points:
+ if ellipse.contains(point):
+ points.append(point)
+
+ if self.divided:
+ points = self.northwest.query_ellipse(ellipse, points)
+ points = self.northeast.query_ellipse(ellipse, points)
+ points = self.southwest.query_ellipse(ellipse, points)
+ points = self.southeast.query_ellipse(ellipse, points)
+
+ return points
+
def nearby_points(
self,
point: Record,
--
GitLab