Skip to content

GitLab

Commit ff50a10e authored by Joseph Siddons's avatar Joseph Siddons
Browse files

feat: Comparison between quad/octtree and ellipse on surface of earth

parent 1499ecd5
Hide whitespace changes
Inline Side-by-side
Showing with 225 additions and 3 deletions
GeoSpatialTools/octtree.py
View file @ ff50a10e
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,
......
GeoSpatialTools/quadtree.py
View file @ ff50a10e
......@@ -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,
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment