Skip to content

GitLab

Commit 762dfafa authored by Joseph Siddons's avatar Joseph Siddons
Browse files

refactor(Rectangle)!: Fix wrapping at -180, 180. Rewrite as dataclass. 

Adds north, east, south, west attributes
parent a10a6cf1
Hide whitespace changes
Inline Side-by-side
Showing with 224 additions and 195 deletions
GeoSpatialTools/octtree.py
View file @ 762dfafa
from datetime import datetime, timedelta
from dataclasses import dataclass
import datetime
from .distance_metrics import haversine, destination
from .utils import LatitudeError
......@@ -83,6 +83,7 @@ class SpaceTimeRecords(list[SpaceTimeRecord]):
"""List of SpaceTimeRecords"""
@dataclass
class SpaceTimeRectangle:
"""
A simple Space Time SpaceTimeRectangle class.
......@@ -104,73 +105,134 @@ class SpaceTimeRectangle:
Horizontal centre of the rectangle (longitude).
lat : float
Vertical centre of the rectangle (latitude).
datetime : datetime
datetime : datetime.datetime
Datetime centre of the rectangle.
w : float
Width of the rectangle (longitude range).
h : float
Height of the rectangle (latitude range).
dt : timedelta
dt : datetime.timedelta
time extent of the rectangle.
"""
def __init__(
self,
lon: float,
lat: float,
datetime: datetime,
lon_range: float,
lat_range: float,
dt: timedelta,
) -> None:
self.lon = lon
self.lat = lat
self.lon_range = lon_range
self.lat_range = lat_range
self.datetime = datetime
self.dt = dt
lon: float
lat: float
date: datetime.datetime
lon_range: float
lat_range: float
dt: datetime.timedelta
def __post_init__(self):
if self.lon > 180:
self.lon -= 360
if self.lat > 90 or self.lat < -90:
raise LatitudeError(
f"Central latitude value out of range {self.lat}, "
+ "should be between -90, 90 degrees"
)
def __str__(self) -> str:
return f"SpaceTimeRectangle(x = {self.lon}, y = {self.lat}, w = {self.lon_range}, h = {self.lat_range}, t = {self.datetime}, dt = {self.dt})"
@property
def west(self) -> float:
"""Western boundary of the Rectangle"""
return (((self.lon - self.lon_range / 2) + 540) % 360) - 180
@property
def east(self) -> float:
"""Eastern boundary of the Rectangle"""
return (((self.lon + self.lon_range / 2) + 540) % 360) - 180
@property
def north(self) -> float:
"""Northern boundary of the Rectangle"""
north = self.lat + self.lat_range / 2
if north > 90:
raise LatitudeError(
"Rectangle crosses north pole - Use two Rectangles"
)
return north
def __eq__(self, other: object) -> bool:
return (
isinstance(other, SpaceTimeRectangle)
and self.lon == other.lon
and self.lat == other.lat
and self.lon_range == other.lon_range
and self.lat_range == other.lat_range
and self.datetime == other.datetime
and self.dt == other.dt
@property
def south(self) -> float:
"""Southern boundary of the Rectangle"""
south = self.lat - self.lat_range / 2
if south < -90:
raise LatitudeError(
"Rectangle crosses south pole - Use two Rectangles"
)
return south
@property
def start(self) -> datetime.datetime:
"""Start date of the Rectangle"""
return self.date - self.dt / 2
@property
def end(self) -> datetime.datetime:
"""End date of the Rectangle"""
return self.date + self.dt / 2
@property
def edge_dist(self) -> float:
"""Approximate maximum distance from the centre to an edge"""
corner_dist = max(
haversine(self.lon, self.lat, self.east, self.north),
haversine(self.lon, self.lat, self.east, self.south),
)
if self.east * self.west < 0:
corner_dist = max(
corner_dist,
haversine(self.lon, self.lat, self.east, 0),
)
return corner_dist
def _test_east_west(self, lon: float) -> bool:
if self.lon_range >= 360:
# Rectangle encircles earth
return True
if self.east > self.lon and self.west < self.lon:
return lon <= self.east and lon >= self.west
if self.east < self.lon:
return not (lon > self.east and lon < self.west)
if self.west > self.lon:
return not (lon < self.east and lon > self.west)
return False
def _test_north_south(self, lat: float) -> bool:
return lat <= self.north and lat >= self.south
def contains(self, point: SpaceTimeRecord) -> bool:
"""Test if a point is contained within the SpaceTimeRectangle"""
return (
point.lon <= self.lon + self.lon_range / 2
and point.lon >= self.lon - self.lon_range / 2
and point.lat <= self.lat + self.lat_range / 2
and point.lat >= self.lat - self.lat_range / 2
and point.datetime <= self.datetime + self.dt / 2
and point.datetime >= self.datetime - self.dt / 2
if point.datetime > self.end or point.datetime < self.start:
return False
return self._test_north_south(point.lat) and self._test_east_west(
point.lon
)
def intersects(self, other: object) -> bool:
"""Test if another Rectangle object intersects this Rectangle"""
return isinstance(other, SpaceTimeRectangle) and not (
self.lon - self.lon_range / 2 > other.lon + other.lon_range / 2
or self.lon + self.lon_range / 2 < other.lon - other.lon_range / 2
or self.lat - self.lat_range / 2 > other.lat + other.lat_range / 2
or self.lat + self.lat_range / 2 < other.lat - other.lat_range / 2
or self.datetime - self.dt / 2 > other.datetime + other.dt / 2
or self.datetime + self.dt / 2 < other.datetime - other.dt / 2
if not isinstance(other, SpaceTimeRectangle):
raise TypeError(
f"other must be a Rectangle class, got {type(other)}"
)
if other.end < self.start or other.start > self.end:
# Not in the same time range
return False
if other.south > self.north:
# Other is fully north of self
return False
if other.north < self.south:
# Other is fully south of self
return False
# Handle east / west edges
return self._test_east_west(other.west) or self._test_east_west(
other.east
)
def nearby(
self,
point: SpaceTimeRecord,
dist: float,
t_dist: timedelta,
t_dist: datetime.timedelta,
) -> bool:
"""
Check if point is nearby the Rectangle
......@@ -192,42 +254,21 @@ class SpaceTimeRectangle:
----------
point : SpaceTimeRecord
dist : float,
t_dist : timedelta
t_dist : datetime.timedelta
Returns
-------
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
point.datetime - t_dist > self.date + self.dt / 2
or point.datetime + t_dist < self.date - self.dt / 2
):
return False
# QUESTION: Is this sufficient? Possibly it is overkill
corner_dist = max(
haversine(
self.lon,
self.lat,
self.lon + self.lon_range / 2,
self.lat + self.lat_range / 2,
),
haversine(
self.lon,
self.lat,
self.lon + self.lon_range / 2,
self.lat - self.lat_range / 2,
),
)
if (self.lat + self.lat_range / 2) * (
self.lat - self.lat_range / 2
) < 0:
corner_dist = max(
corner_dist,
haversine(self.lon, self.lat, self.lon + self.lon_range / 2, 0),
)
return (
haversine(self.lon, self.lat, point.lon, point.lat)
<= dist + corner_dist
<= dist + self.edge_dist
)
......@@ -241,7 +282,7 @@ class SpaceTimeEllipse:
Horizontal centre of the ellipse
lat : float
Vertical centre of the ellipse
datetime : datetime
datetime : datetime.datetime
Datetime centre of the ellipse.
a : float
Length of the semi-major axis
......@@ -249,7 +290,7 @@ class SpaceTimeEllipse:
Length of the semi-minor axis
theta : float
Angle of the semi-major axis from horizontal anti-clockwise in radians
dt : timedelta
dt : datetime.timedelta
(full) time extent of the ellipse.
"""
......@@ -257,11 +298,11 @@ class SpaceTimeEllipse:
self,
lon: float,
lat: float,
datetime: datetime,
datetime: datetime.datetime,
a: float,
b: float,
theta: float,
dt: timedelta,
dt: datetime.timedelta,
) -> None:
self.a = a
self.b = b
......@@ -290,53 +331,28 @@ class SpaceTimeEllipse:
(self.bearing - 180) % 360,
self.c,
)
self.start = self.datetime - self.dt / 2
self.end = self.datetime + self.dt / 2
def contains(self, point: SpaceTimeRecord) -> bool:
"""Test if a point is contained within the Ellipse"""
if point.datetime > self.end or point.datetime < self.start:
return False
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
)
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: 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
):
if rect.start > self.end or rect.end < self.start:
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
<= rect.edge_dist + self.a
and haversine(self.p2_lon, self.p2_lat, rect.lon, rect.lat)
<= corner_dist + self.a
<= rect.edge_dist + self.a
)
......@@ -419,7 +435,7 @@ class OctTree:
SpaceTimeRectangle(
self.boundary.lon - self.boundary.lon_range / 4,
self.boundary.lat + self.boundary.lat_range / 4,
self.boundary.datetime + self.boundary.dt / 4,
self.boundary.date + self.boundary.dt / 4,
self.boundary.lon_range / 2,
self.boundary.lat_range / 2,
self.boundary.dt / 2,
......@@ -432,7 +448,7 @@ class OctTree:
SpaceTimeRectangle(
self.boundary.lon + self.boundary.lon_range / 4,
self.boundary.lat + self.boundary.lat_range / 4,
self.boundary.datetime + self.boundary.dt / 4,
self.boundary.date + self.boundary.dt / 4,
self.boundary.lon_range / 2,
self.boundary.lat_range / 2,
self.boundary.dt / 2,
......@@ -445,7 +461,7 @@ class OctTree:
SpaceTimeRectangle(
self.boundary.lon - self.boundary.lon_range / 4,
self.boundary.lat - self.boundary.lat_range / 4,
self.boundary.datetime + self.boundary.dt / 4,
self.boundary.date + self.boundary.dt / 4,
self.boundary.lon_range / 2,
self.boundary.lat_range / 2,
self.boundary.dt / 2,
......@@ -458,7 +474,7 @@ class OctTree:
SpaceTimeRectangle(
self.boundary.lon + self.boundary.lon_range / 4,
self.boundary.lat - self.boundary.lat_range / 4,
self.boundary.datetime + self.boundary.dt / 4,
self.boundary.date + self.boundary.dt / 4,
self.boundary.lon_range / 2,
self.boundary.lat_range / 2,
self.boundary.dt / 2,
......@@ -471,7 +487,7 @@ class OctTree:
SpaceTimeRectangle(
self.boundary.lon - self.boundary.lon_range / 4,
self.boundary.lat + self.boundary.lat_range / 4,
self.boundary.datetime - self.boundary.dt / 4,
self.boundary.date - self.boundary.dt / 4,
self.boundary.lon_range / 2,
self.boundary.lat_range / 2,
self.boundary.dt / 2,
......@@ -484,7 +500,7 @@ class OctTree:
SpaceTimeRectangle(
self.boundary.lon + self.boundary.lon_range / 4,
self.boundary.lat + self.boundary.lat_range / 4,
self.boundary.datetime - self.boundary.dt / 4,
self.boundary.date - self.boundary.dt / 4,
self.boundary.lon_range / 2,
self.boundary.lat_range / 2,
self.boundary.dt / 2,
......@@ -497,7 +513,7 @@ class OctTree:
SpaceTimeRectangle(
self.boundary.lon - self.boundary.lon_range / 4,
self.boundary.lat - self.boundary.lat_range / 4,
self.boundary.datetime - self.boundary.dt / 4,
self.boundary.date - self.boundary.dt / 4,
self.boundary.lon_range / 2,
self.boundary.lat_range / 2,
self.boundary.dt / 2,
......@@ -510,7 +526,7 @@ class OctTree:
SpaceTimeRectangle(
self.boundary.lon + self.boundary.lon_range / 4,
self.boundary.lat - self.boundary.lat_range / 4,
self.boundary.datetime - self.boundary.dt / 4,
self.boundary.date - self.boundary.dt / 4,
self.boundary.lon_range / 2,
self.boundary.lat_range / 2,
self.boundary.dt / 2,
......@@ -522,7 +538,7 @@ class OctTree:
self.divided = True
def _datetime_is_numeric(self) -> bool:
return not isinstance(self.boundary.datetime, datetime)
return not isinstance(self.boundary.date, datetime)
def insert(self, point: SpaceTimeRecord) -> bool:
"""
......@@ -618,7 +634,7 @@ class OctTree:
self,
point: SpaceTimeRecord,
dist: float,
t_dist: timedelta,
t_dist: datetime.timedelta,
points: SpaceTimeRecords | None = None,
) -> SpaceTimeRecords:
"""
......@@ -637,7 +653,7 @@ class OctTree:
The distance for comparison. Note that Haversine distance is used
as the distance metric as the query SpaceTimeRecord and OctTree are
assumed to lie on the surface of Earth.
t_dist : timedelta
t_dist : datetime.timedelta
Max time gap between SpaceTimeRecords within the OctTree and the
query SpaceTimeRecord. Can be numeric if the OctTree boundaries,
SpaceTimeRecords, and query SpaceTimeRecord have numeric datetime
......
GeoSpatialTools/quadtree.py
View file @ 762dfafa
......@@ -3,6 +3,7 @@ Constuctors for QuadTree classes that can decrease the number of comparisons
for detecting nearby records for example
"""
from dataclasses import dataclass
from datetime import datetime
from .distance_metrics import haversine, destination
from .utils import LatitudeError
......@@ -82,6 +83,7 @@ class Record:
return haversine(self.lon, self.lat, other.lon, other.lat)
@dataclass
class Rectangle:
"""
A simple Rectangle class
......@@ -98,46 +100,100 @@ class Rectangle:
Height of the rectangle
"""
def __init__(
self,
lon: float,
lat: float,
lon_range: float,
lat_range: float,
) -> None:
self.lon = lon
self.lat = lat
self.lon_range = lon_range
self.lat_range = lat_range
lon: float
lat: float
lon_range: float
lat_range: float
def __str__(self) -> str:
return f"Rectangle(x = {self.lon}, y = {self.lat}, w = {self.lon_range}, h = {self.lat_range})"
def __post_init__(self):
if self.lon > 180:
self.lon -= 360
if self.lat > 90 or self.lat < -90:
raise LatitudeError(
f"Central latitude value out of range {self.lat}, "
+ "should be between -90, 90 degrees"
)
def __eq__(self, other: object) -> bool:
return (
isinstance(other, Rectangle)
and self.lon == other.lon
and self.lat == other.lat
and self.lon_range == other.lon_range
and self.lat_range == other.lat_range
@property
def west(self) -> float:
"""Western boundary of the Rectangle"""
return (((self.lon - self.lon_range / 2) + 540) % 360) - 180
@property
def east(self) -> float:
"""Eastern boundary of the Rectangle"""
return (((self.lon + self.lon_range / 2) + 540) % 360) - 180
@property
def north(self) -> float:
"""Northern boundary of the Rectangle"""
north = self.lat + self.lat_range / 2
if north > 90:
raise LatitudeError(
"Rectangle crosses north pole - Use two Rectangles"
)
return north
@property
def south(self) -> float:
"""Southern boundary of the Rectangle"""
south = self.lat - self.lat_range / 2
if south < -90:
raise LatitudeError(
"Rectangle crosses south pole - Use two Rectangles"
)
return south
@property
def edge_dist(self) -> float:
"""Approximate maximum distance from the centre to an edge"""
corner_dist = max(
haversine(self.lon, self.lat, self.east, self.north),
haversine(self.lon, self.lat, self.east, self.south),
)
if self.east * self.west < 0:
corner_dist = max(
corner_dist,
haversine(self.lon, self.lat, self.east, 0),
)
return corner_dist
def _test_east_west(self, lon: float) -> bool:
if self.lon_range >= 360:
# Rectangle encircles earth
return True
if self.east > self.lon and self.west < self.lon:
return lon <= self.east and lon >= self.west
if self.east < self.lon:
return not (lon > self.east and lon < self.west)
if self.west > self.lon:
return not (lon < self.east and lon > self.west)
return False
def _test_north_south(self, lat: float) -> bool:
return lat <= self.north and lat >= self.south
def contains(self, point: Record) -> bool:
"""Test if a point is contained within the Rectangle"""
return (
point.lon <= self.lon + self.lon_range / 2
and point.lon >= self.lon - self.lon_range / 2
and point.lat <= self.lat + self.lat_range / 2
and point.lat >= self.lat - self.lat_range / 2
return self._test_north_south(point.lat) and self._test_east_west(
point.lon
)
def intersects(self, other: object) -> bool:
"""Test if another Rectangle object intersects this Rectangle"""
return isinstance(other, Rectangle) and not (
self.lon - self.lon_range / 2 > other.lon + other.lon_range / 2
or self.lon + self.lon_range / 2 < other.lon - other.lon_range / 2
or self.lat - self.lat_range / 2 > other.lat + other.lat_range / 2
or self.lat + self.lat_range / 2 < other.lat - other.lat_range / 2
if not isinstance(other, Rectangle):
raise TypeError(
f"other must be a Rectangle class, got {type(other)}"
)
if other.south > self.north:
# Other is fully north of self
return False
if other.north < self.south:
# Other is fully south of self
return False
# Handle east / west edges
return self._test_east_west(other.west) or self._test_east_west(
other.east
)
def nearby(
......@@ -147,30 +203,9 @@ class Rectangle:
) -> bool:
"""Check if point is nearby the Rectangle"""
# QUESTION: Is this sufficient? Possibly it is overkill
corner_dist = max(
haversine(
self.lon,
self.lat,
self.lon + self.lon_range / 2,
self.lat + self.lat_range / 2,
),
haversine(
self.lon,
self.lat,
self.lon + self.lon_range / 2,
self.lat - self.lat_range / 2,
),
)
if (self.lat + self.lat_range / 2) * (
self.lat - self.lat_range / 2
) < 0:
corner_dist = max(
corner_dist,
haversine(self.lon, self.lat, self.lon + self.lon_range / 2, 0),
)
return (
haversine(self.lon, self.lat, point.lon, point.lat)
<= dist + corner_dist
<= dist + self.edge_dist
)
......@@ -235,33 +270,11 @@ class Ellipse:
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
<= rect.edge_dist + self.a
and haversine(self.p2_lon, self.p2_lat, rect.lon, rect.lat)
<= corner_dist + self.a
<= rect.edge_dist + self.a
)
......
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