Skip to content

GitLab

Commit 77fb3001 authored by Joseph Siddons's avatar Joseph Siddons
Browse files

Merge branch 'rect_lon_wrap' into 'main' 

Fix Rectangle query for QuadTree and OctTree classes

See merge request nocsurfaceprocesses/geospatialtools!8
parents 85e6d8af 6b1e0751
Hide whitespace changes
Inline Side-by-side
Showing with 512 additions and 297 deletions
GeoSpatialTools/kdtree.py
View file @ 77fb3001
...@@ -49,12 +49,6 @@ class KDTree: ...@@ -49,12 +49,6 @@ class KDTree:
points.sort(key=lambda p: getattr(p, self.variable)) points.sort(key=lambda p: getattr(p, self.variable))
split_index = n_points // 2 split_index = n_points // 2
self.partition_value = getattr(points[split_index - 1], self.variable) self.partition_value = getattr(points[split_index - 1], self.variable)
# while (
# split_index < n_points
# and getattr(points[split_index], self.variable)
# == self.partition_value
# ):
# split_index += 1
self.split = True self.split = True
...@@ -119,9 +113,9 @@ class KDTree: ...@@ -119,9 +113,9 @@ class KDTree:
def query(self, point) -> tuple[list[Record], float]: def query(self, point) -> tuple[list[Record], float]:
"""Find the nearest Record within the KDTree to a query Record""" """Find the nearest Record within the KDTree to a query Record"""
if point.lon < 0: if point.lon < 0:
point2 = Record(point.lon + 360, point.lat) point2 = Record(point.lon + 360, point.lat, fix_lon=False)
else: else:
point2 = Record(point.lon - 360, point.lat) point2 = Record(point.lon - 360, point.lat, fix_lon=False)
r1, d1 = self._query(point) r1, d1 = self._query(point)
r2, d2 = self._query(point2) r2, d2 = self._query(point2)
......
GeoSpatialTools/neighbours.py
View file @ 77fb3001
...@@ -11,10 +11,14 @@ Numeric = TypeVar("Numeric", int, float, datetime, date) ...@@ -11,10 +11,14 @@ Numeric = TypeVar("Numeric", int, float, datetime, date)
class SortedWarning(Warning): class SortedWarning(Warning):
"""Warning class for Sortedness"""
pass pass
class SortedError(Exception): class SortedError(Exception):
"""Error class for Sortedness"""
pass pass
......
GeoSpatialTools/octtree.py
View file @ 77fb3001
from datetime import datetime, timedelta from dataclasses import dataclass
from .distance_metrics import haversine import datetime
from .distance_metrics import haversine, destination from .distance_metrics import haversine, destination
from .utils import LatitudeError
from math import degrees, sqrt from math import degrees, sqrt
...@@ -25,12 +26,14 @@ class SpaceTimeRecord: ...@@ -25,12 +26,14 @@ class SpaceTimeRecord:
Horizontal coordinate (longitude). Horizontal coordinate (longitude).
lat : float lat : float
Vertical coordinate (latitude). Vertical coordinate (latitude).
datetime : datetime datetime : datetime.datetime
Datetime of the record. Can also be a numeric value such as pentad. Datetime of the record. Can also be a numeric value such as pentad.
Comparisons between Records with datetime and Records with numeric Comparisons between Records with datetime and Records with numeric
datetime will fail. datetime will fail.
uid : str | None uid : str | None
Unique Identifier. Unique Identifier.
fix_lon : bool
Force longitude to -180, 180
**data **data
Additional data passed to the SpaceTimeRecord for use by other functions Additional data passed to the SpaceTimeRecord for use by other functions
or classes. or classes.
...@@ -40,11 +43,19 @@ class SpaceTimeRecord: ...@@ -40,11 +43,19 @@ class SpaceTimeRecord:
self, self,
lon: float, lon: float,
lat: float, lat: float,
datetime: datetime, datetime: datetime.datetime,
uid: str | None = None, uid: str | None = None,
fix_lon: bool = True,
**data, **data,
) -> None: ) -> None:
self.lon = lon self.lon = lon
if fix_lon:
# Move lon to -180, 180
self.lon = ((self.lon + 540) % 360) - 180
if lat < -90 or lat > 90:
raise LatitudeError(
"Expected latitude value to be between -90 and 90 degrees"
)
self.lat = lat self.lat = lat
self.datetime = datetime self.datetime = datetime
self.uid = uid self.uid = uid
...@@ -53,12 +64,15 @@ class SpaceTimeRecord: ...@@ -53,12 +64,15 @@ class SpaceTimeRecord:
return None return None
def __str__(self) -> str: def __str__(self) -> str:
return f"Record(x = {self.lon}, y = {self.lat}, datetime = {self.datetime}, uid = {self.uid})" return f"SpaceTimeRecord(x = {self.lon}, y = {self.lat}, datetime = {self.datetime}, uid = {self.uid})"
def __eq__(self, other: object) -> bool: def __eq__(self, other: object) -> bool:
if not isinstance(other, SpaceTimeRecord):
return False
if self.uid and other.uid:
return self.uid == other.uid
return ( return (
isinstance(other, SpaceTimeRecord) self.lon == other.lon
and self.lon == other.lon
and self.lat == other.lat and self.lat == other.lat
and self.datetime == other.datetime and self.datetime == other.datetime
and (not (self.uid or other.uid) or self.uid == other.uid) and (not (self.uid or other.uid) or self.uid == other.uid)
...@@ -69,6 +83,7 @@ class SpaceTimeRecords(list[SpaceTimeRecord]): ...@@ -69,6 +83,7 @@ class SpaceTimeRecords(list[SpaceTimeRecord]):
"""List of SpaceTimeRecords""" """List of SpaceTimeRecords"""
@dataclass
class SpaceTimeRectangle: class SpaceTimeRectangle:
""" """
A simple Space Time SpaceTimeRectangle class. A simple Space Time SpaceTimeRectangle class.
...@@ -90,73 +105,134 @@ class SpaceTimeRectangle: ...@@ -90,73 +105,134 @@ class SpaceTimeRectangle:
Horizontal centre of the rectangle (longitude). Horizontal centre of the rectangle (longitude).
lat : float lat : float
Vertical centre of the rectangle (latitude). Vertical centre of the rectangle (latitude).
datetime : datetime datetime : datetime.datetime
Datetime centre of the rectangle. Datetime centre of the rectangle.
w : float w : float
Width of the rectangle (longitude range). Width of the rectangle (longitude range).
h : float h : float
Height of the rectangle (latitude range). Height of the rectangle (latitude range).
dt : timedelta dt : datetime.timedelta
time extent of the rectangle. time extent of the rectangle.
""" """
def __init__( lon: float
self, lat: float
lon: float, date: datetime.datetime
lat: float, lon_range: float
datetime: datetime, lat_range: float
lon_range: float, dt: datetime.timedelta
lat_range: float,
dt: timedelta, def __post_init__(self):
) -> None: if self.lon > 180:
self.lon = lon self.lon -= 360
self.lat = lat if self.lat > 90 or self.lat < -90:
self.lon_range = lon_range raise LatitudeError(
self.lat_range = lat_range f"Central latitude value out of range {self.lat}, "
self.datetime = datetime + "should be between -90, 90 degrees"
self.dt = dt )
def __str__(self) -> str: @property
return f"SpaceTimeRectangle(x = {self.lon}, y = {self.lat}, w = {self.lon_range}, h = {self.lat_range}, t = {self.datetime}, dt = {self.dt})" 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
def __eq__(self, other: object) -> bool: @property
return ( def start(self) -> datetime.datetime:
isinstance(other, SpaceTimeRectangle) """Start date of the Rectangle"""
and self.lon == other.lon return self.date - self.dt / 2
and self.lat == other.lat
and self.lon_range == other.lon_range @property
and self.lat_range == other.lat_range def end(self) -> datetime.datetime:
and self.datetime == other.datetime """End date of the Rectangle"""
and self.dt == other.dt 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.north * self.south < 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: def contains(self, point: SpaceTimeRecord) -> bool:
"""Test if a point is contained within the SpaceTimeRectangle""" """Test if a point is contained within the SpaceTimeRectangle"""
return ( if point.datetime > self.end or point.datetime < self.start:
point.lon <= self.lon + self.lon_range / 2 return False
and point.lon >= self.lon - self.lon_range / 2 return self._test_north_south(point.lat) and self._test_east_west(
and point.lat <= self.lat + self.lat_range / 2 point.lon
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
) )
def intersects(self, other: object) -> bool: def intersects(self, other: object) -> bool:
"""Test if another Rectangle object intersects this Rectangle""" """Test if another Rectangle object intersects this Rectangle"""
return isinstance(other, SpaceTimeRectangle) and not ( if not isinstance(other, SpaceTimeRectangle):
self.lon - self.lon_range / 2 > other.lon + other.lon_range / 2 raise TypeError(
or self.lon + self.lon_range / 2 < other.lon - other.lon_range / 2 f"other must be a Rectangle class, got {type(other)}"
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 other.end < self.start or other.start > self.end:
or self.datetime - self.dt / 2 > other.datetime + other.dt / 2 # Not in the same time range
or self.datetime + self.dt / 2 < other.datetime - other.dt / 2 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( def nearby(
self, self,
point: SpaceTimeRecord, point: SpaceTimeRecord,
dist: float, dist: float,
t_dist: timedelta, t_dist: datetime.timedelta,
) -> bool: ) -> bool:
""" """
Check if point is nearby the Rectangle Check if point is nearby the Rectangle
...@@ -178,42 +254,21 @@ class SpaceTimeRectangle: ...@@ -178,42 +254,21 @@ class SpaceTimeRectangle:
---------- ----------
point : SpaceTimeRecord point : SpaceTimeRecord
dist : float, dist : float,
t_dist : timedelta t_dist : datetime.timedelta
Returns Returns
------- -------
bool : True if the point is <= dist + max(dist(centre, corners)) bool : True if the point is <= dist + max(dist(centre, corners))
""" """
if ( if (
point.datetime - t_dist > self.datetime + self.dt / 2 point.datetime - t_dist > self.date + self.dt / 2
or point.datetime + t_dist < self.datetime - self.dt / 2 or point.datetime + t_dist < self.date - self.dt / 2
): ):
return False return False
# QUESTION: Is this sufficient? Possibly it is overkill # 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 ( return (
haversine(self.lon, self.lat, point.lon, point.lat) haversine(self.lon, self.lat, point.lon, point.lat)
<= dist + corner_dist <= dist + self.edge_dist
) )
...@@ -227,7 +282,7 @@ class SpaceTimeEllipse: ...@@ -227,7 +282,7 @@ class SpaceTimeEllipse:
Horizontal centre of the ellipse Horizontal centre of the ellipse
lat : float lat : float
Vertical centre of the ellipse Vertical centre of the ellipse
datetime : datetime datetime : datetime.datetime
Datetime centre of the ellipse. Datetime centre of the ellipse.
a : float a : float
Length of the semi-major axis Length of the semi-major axis
...@@ -235,7 +290,7 @@ class SpaceTimeEllipse: ...@@ -235,7 +290,7 @@ class SpaceTimeEllipse:
Length of the semi-minor axis Length of the semi-minor axis
theta : float theta : float
Angle of the semi-major axis from horizontal anti-clockwise in radians Angle of the semi-major axis from horizontal anti-clockwise in radians
dt : timedelta dt : datetime.timedelta
(full) time extent of the ellipse. (full) time extent of the ellipse.
""" """
...@@ -243,11 +298,11 @@ class SpaceTimeEllipse: ...@@ -243,11 +298,11 @@ class SpaceTimeEllipse:
self, self,
lon: float, lon: float,
lat: float, lat: float,
datetime: datetime, datetime: datetime.datetime,
a: float, a: float,
b: float, b: float,
theta: float, theta: float,
dt: timedelta, dt: datetime.timedelta,
) -> None: ) -> None:
self.a = a self.a = a
self.b = b self.b = b
...@@ -276,53 +331,28 @@ class SpaceTimeEllipse: ...@@ -276,53 +331,28 @@ class SpaceTimeEllipse:
(self.bearing - 180) % 360, (self.bearing - 180) % 360,
self.c, self.c,
) )
self.start = self.datetime - self.dt / 2
self.end = self.datetime + self.dt / 2
def contains(self, point: SpaceTimeRecord) -> bool: def contains(self, point: SpaceTimeRecord) -> bool:
"""Test if a point is contained within the Ellipse""" """Test if a point is contained within the Ellipse"""
if point.datetime > self.end or point.datetime < self.start:
return False
return ( return (
( haversine(self.p1_lon, self.p1_lat, point.lon, point.lat)
haversine(self.p1_lon, self.p1_lat, point.lon, point.lat) + haversine(self.p2_lon, self.p2_lat, point.lon, point.lat)
+ haversine(self.p2_lon, self.p2_lat, point.lon, point.lat) ) <= 2 * self.a
)
<= 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: def nearby_rect(self, rect: SpaceTimeRectangle) -> bool:
"""Test if a rectangle is near to the Ellipse""" """Test if a rectangle is near to the Ellipse"""
if ( if rect.start > self.end or rect.end < self.start:
rect.datetime - rect.dt / 2 > self.datetime + self.dt / 2
or rect.datetime + rect.dt / 2 < self.datetime - self.dt / 2
):
return False return False
# TODO: Check corners, and 0 lat # 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 ( return (
haversine(self.p1_lon, self.p1_lat, rect.lon, rect.lat) 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) and haversine(self.p2_lon, self.p2_lat, rect.lon, rect.lat)
<= corner_dist + self.a <= rect.edge_dist + self.a
) )
...@@ -405,7 +435,7 @@ class OctTree: ...@@ -405,7 +435,7 @@ class OctTree:
SpaceTimeRectangle( SpaceTimeRectangle(
self.boundary.lon - self.boundary.lon_range / 4, self.boundary.lon - self.boundary.lon_range / 4,
self.boundary.lat + self.boundary.lat_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.lon_range / 2,
self.boundary.lat_range / 2, self.boundary.lat_range / 2,
self.boundary.dt / 2, self.boundary.dt / 2,
...@@ -418,7 +448,7 @@ class OctTree: ...@@ -418,7 +448,7 @@ class OctTree:
SpaceTimeRectangle( SpaceTimeRectangle(
self.boundary.lon + self.boundary.lon_range / 4, self.boundary.lon + self.boundary.lon_range / 4,
self.boundary.lat + self.boundary.lat_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.lon_range / 2,
self.boundary.lat_range / 2, self.boundary.lat_range / 2,
self.boundary.dt / 2, self.boundary.dt / 2,
...@@ -431,7 +461,7 @@ class OctTree: ...@@ -431,7 +461,7 @@ class OctTree:
SpaceTimeRectangle( SpaceTimeRectangle(
self.boundary.lon - self.boundary.lon_range / 4, self.boundary.lon - self.boundary.lon_range / 4,
self.boundary.lat - self.boundary.lat_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.lon_range / 2,
self.boundary.lat_range / 2, self.boundary.lat_range / 2,
self.boundary.dt / 2, self.boundary.dt / 2,
...@@ -444,7 +474,7 @@ class OctTree: ...@@ -444,7 +474,7 @@ class OctTree:
SpaceTimeRectangle( SpaceTimeRectangle(
self.boundary.lon + self.boundary.lon_range / 4, self.boundary.lon + self.boundary.lon_range / 4,
self.boundary.lat - self.boundary.lat_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.lon_range / 2,
self.boundary.lat_range / 2, self.boundary.lat_range / 2,
self.boundary.dt / 2, self.boundary.dt / 2,
...@@ -457,7 +487,7 @@ class OctTree: ...@@ -457,7 +487,7 @@ class OctTree:
SpaceTimeRectangle( SpaceTimeRectangle(
self.boundary.lon - self.boundary.lon_range / 4, self.boundary.lon - self.boundary.lon_range / 4,
self.boundary.lat + self.boundary.lat_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.lon_range / 2,
self.boundary.lat_range / 2, self.boundary.lat_range / 2,
self.boundary.dt / 2, self.boundary.dt / 2,
...@@ -470,7 +500,7 @@ class OctTree: ...@@ -470,7 +500,7 @@ class OctTree:
SpaceTimeRectangle( SpaceTimeRectangle(
self.boundary.lon + self.boundary.lon_range / 4, self.boundary.lon + self.boundary.lon_range / 4,
self.boundary.lat + self.boundary.lat_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.lon_range / 2,
self.boundary.lat_range / 2, self.boundary.lat_range / 2,
self.boundary.dt / 2, self.boundary.dt / 2,
...@@ -483,7 +513,7 @@ class OctTree: ...@@ -483,7 +513,7 @@ class OctTree:
SpaceTimeRectangle( SpaceTimeRectangle(
self.boundary.lon - self.boundary.lon_range / 4, self.boundary.lon - self.boundary.lon_range / 4,
self.boundary.lat - self.boundary.lat_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.lon_range / 2,
self.boundary.lat_range / 2, self.boundary.lat_range / 2,
self.boundary.dt / 2, self.boundary.dt / 2,
...@@ -496,7 +526,7 @@ class OctTree: ...@@ -496,7 +526,7 @@ class OctTree:
SpaceTimeRectangle( SpaceTimeRectangle(
self.boundary.lon + self.boundary.lon_range / 4, self.boundary.lon + self.boundary.lon_range / 4,
self.boundary.lat - self.boundary.lat_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.lon_range / 2,
self.boundary.lat_range / 2, self.boundary.lat_range / 2,
self.boundary.dt / 2, self.boundary.dt / 2,
...@@ -508,7 +538,7 @@ class OctTree: ...@@ -508,7 +538,7 @@ class OctTree:
self.divided = True self.divided = True
def _datetime_is_numeric(self) -> bool: 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: def insert(self, point: SpaceTimeRecord) -> bool:
""" """
...@@ -604,7 +634,7 @@ class OctTree: ...@@ -604,7 +634,7 @@ class OctTree:
self, self,
point: SpaceTimeRecord, point: SpaceTimeRecord,
dist: float, dist: float,
t_dist: timedelta, t_dist: datetime.timedelta,
points: SpaceTimeRecords | None = None, points: SpaceTimeRecords | None = None,
) -> SpaceTimeRecords: ) -> SpaceTimeRecords:
""" """
...@@ -623,7 +653,7 @@ class OctTree: ...@@ -623,7 +653,7 @@ class OctTree:
The distance for comparison. Note that Haversine distance is used The distance for comparison. Note that Haversine distance is used
as the distance metric as the query SpaceTimeRecord and OctTree are as the distance metric as the query SpaceTimeRecord and OctTree are
assumed to lie on the surface of Earth. 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 Max time gap between SpaceTimeRecords within the OctTree and the
query SpaceTimeRecord. Can be numeric if the OctTree boundaries, query SpaceTimeRecord. Can be numeric if the OctTree boundaries,
SpaceTimeRecords, and query SpaceTimeRecord have numeric datetime SpaceTimeRecords, and query SpaceTimeRecord have numeric datetime
......
GeoSpatialTools/quadtree.py
View file @ 77fb3001
...@@ -3,8 +3,10 @@ Constuctors for QuadTree classes that can decrease the number of comparisons ...@@ -3,8 +3,10 @@ Constuctors for QuadTree classes that can decrease the number of comparisons
for detecting nearby records for example for detecting nearby records for example
""" """
from dataclasses import dataclass
from datetime import datetime from datetime import datetime
from .distance_metrics import haversine, destination from .distance_metrics import haversine, destination
from .utils import LatitudeError
from math import degrees, sqrt from math import degrees, sqrt
...@@ -12,6 +14,12 @@ class Record: ...@@ -12,6 +14,12 @@ class Record:
""" """
ICOADS Record class ICOADS Record class
This is a simple instance of an ICOARDS record, it requires position data.
It can optionally include datetime, a UID, and extra data passed as
keyword arguments.
Equality is checked only on the required fields + UID if it is specified.
Parameters Parameters
---------- ----------
lon : float lon : float
...@@ -22,6 +30,11 @@ class Record: ...@@ -22,6 +30,11 @@ class Record:
Datetime of the record Datetime of the record
uid : str | None uid : str | None
Unique Identifier Unique Identifier
fix_lon : bool
Force longitude to -180, 180
**data
Additional data passed to the Record for use by other functions or
classes.
""" """
def __init__( def __init__(
...@@ -30,9 +43,17 @@ class Record: ...@@ -30,9 +43,17 @@ class Record:
lat: float, lat: float,
datetime: datetime | None = None, datetime: datetime | None = None,
uid: str | None = None, uid: str | None = None,
fix_lon: bool = True,
**data, **data,
) -> None: ) -> None:
self.lon = lon self.lon = lon
if fix_lon:
# Move lon to -180, 180
self.lon = ((self.lon + 540) % 360) - 180
if lat < -90 or lat > 90:
raise LatitudeError(
"Expected latitude value to be between -90 and 90 degrees"
)
self.lat = lat self.lat = lat
self.datetime = datetime self.datetime = datetime
self.uid = uid self.uid = uid
...@@ -44,9 +65,12 @@ class Record: ...@@ -44,9 +65,12 @@ class Record:
return f"Record(lon = {self.lon}, lat = {self.lat}, datetime = {self.datetime}, uid = {self.uid})" return f"Record(lon = {self.lon}, lat = {self.lat}, datetime = {self.datetime}, uid = {self.uid})"
def __eq__(self, other: object) -> bool: def __eq__(self, other: object) -> bool:
if not isinstance(other, Record):
return False
if self.uid and other.uid:
return self.uid == other.uid
return ( return (
isinstance(other, Record) self.lon == other.lon
and self.lon == other.lon
and self.lat == other.lat and self.lat == other.lat
and self.datetime == other.datetime and self.datetime == other.datetime
and (not (self.uid or other.uid) or self.uid == other.uid) and (not (self.uid or other.uid) or self.uid == other.uid)
...@@ -59,6 +83,7 @@ class Record: ...@@ -59,6 +83,7 @@ class Record:
return haversine(self.lon, self.lat, other.lon, other.lat) return haversine(self.lon, self.lat, other.lon, other.lat)
@dataclass
class Rectangle: class Rectangle:
""" """
A simple Rectangle class A simple Rectangle class
...@@ -75,46 +100,100 @@ class Rectangle: ...@@ -75,46 +100,100 @@ class Rectangle:
Height of the rectangle Height of the rectangle
""" """
def __init__( lon: float
self, lat: float
lon: float, lon_range: float
lat: float, lat_range: float
lon_range: float,
lat_range: float, def __post_init__(self):
) -> None: if self.lon > 180:
self.lon = lon self.lon -= 360
self.lat = lat if self.lat > 90 or self.lat < -90:
self.lon_range = lon_range raise LatitudeError(
self.lat_range = lat_range f"Central latitude value out of range {self.lat}, "
+ "should be between -90, 90 degrees"
)
def __str__(self) -> str: @property
return f"Rectangle(x = {self.lon}, y = {self.lat}, w = {self.lon_range}, h = {self.lat_range})" 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
def __eq__(self, other: object) -> bool: @property
return ( def edge_dist(self) -> float:
isinstance(other, Rectangle) """Approximate maximum distance from the centre to an edge"""
and self.lon == other.lon corner_dist = max(
and self.lat == other.lat haversine(self.lon, self.lat, self.east, self.north),
and self.lon_range == other.lon_range haversine(self.lon, self.lat, self.east, self.south),
and self.lat_range == other.lat_range
) )
if self.north * self.south < 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: def contains(self, point: Record) -> bool:
"""Test if a point is contained within the Rectangle""" """Test if a point is contained within the Rectangle"""
return ( return self._test_north_south(point.lat) and self._test_east_west(
point.lon <= self.lon + self.lon_range / 2 point.lon
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
) )
def intersects(self, other: object) -> bool: def intersects(self, other: object) -> bool:
"""Test if another Rectangle object intersects this Rectangle""" """Test if another Rectangle object intersects this Rectangle"""
return isinstance(other, Rectangle) and not ( if not isinstance(other, Rectangle):
self.lon - self.lon_range / 2 > other.lon + other.lon_range / 2 raise TypeError(
or self.lon + self.lon_range / 2 < other.lon - other.lon_range / 2 f"other must be a Rectangle class, got {type(other)}"
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 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( def nearby(
...@@ -124,30 +203,9 @@ class Rectangle: ...@@ -124,30 +203,9 @@ class Rectangle:
) -> bool: ) -> bool:
"""Check if point is nearby the Rectangle""" """Check if point is nearby the Rectangle"""
# QUESTION: Is this sufficient? Possibly it is overkill # 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 ( return (
haversine(self.lon, self.lat, point.lon, point.lat) haversine(self.lon, self.lat, point.lon, point.lat)
<= dist + corner_dist <= dist + self.edge_dist
) )
...@@ -212,33 +270,11 @@ class Ellipse: ...@@ -212,33 +270,11 @@ class Ellipse:
def nearby_rect(self, rect: Rectangle) -> bool: def nearby_rect(self, rect: Rectangle) -> bool:
"""Test if a rectangle is near to the Ellipse""" """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 ( return (
haversine(self.p1_lon, self.p1_lat, rect.lon, rect.lat) 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) and haversine(self.p2_lon, self.p2_lat, rect.lon, rect.lat)
<= corner_dist + self.a <= rect.edge_dist + self.a
) )
......
GeoSpatialTools/utils.py 0 → 100644
View file @ 77fb3001
class LatitudeError(ValueError):
"""Error for invalid Latitude Value"""
pass
notebooks/kdtree.ipynb
View file @ 77fb3001
...@@ -48,7 +48,7 @@ ...@@ -48,7 +48,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 11, "execution_count": 4,
"id": "c60b30de-f864-477a-a09a-5f1caa4d9b9a", "id": "c60b30de-f864-477a-a09a-5f1caa4d9b9a",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
...@@ -63,11 +63,11 @@ ...@@ -63,11 +63,11 @@
"│ --- ┆ --- │\n", "│ --- ┆ --- │\n",
"│ i64 ┆ i64 │\n", "│ i64 ┆ i64 │\n",
"╞══════╪═════╡\n", "╞══════╪═════╡\n",
"│ 16 ┆ -75 │\n", "│ 127 ┆ 21 │\n",
"│ 144 ┆ -77 │\n", "│ -148 ┆ 36 │\n",
"│ -173 ┆ -83 │\n", "│ -46 ┆ -15 │\n",
"│ 142 ┆ -81 │\n", "│ 104 ┆ 89 │\n",
"│ -50 ┆ -38 │\n", "│ -57 ┆ -31 │\n",
"└──────┴─────┘\n" "└──────┴─────┘\n"
] ]
} }
...@@ -90,7 +90,7 @@ ...@@ -90,7 +90,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 12, "execution_count": 5,
"id": "875f2a67-49fe-476f-add1-b1d76c6cd8f9", "id": "875f2a67-49fe-476f-add1-b1d76c6cd8f9",
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
...@@ -100,7 +100,7 @@ ...@@ -100,7 +100,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 13, "execution_count": 6,
"id": "1e883e5a-5086-4c29-aff2-d308874eae16", "id": "1e883e5a-5086-4c29-aff2-d308874eae16",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
...@@ -108,8 +108,8 @@ ...@@ -108,8 +108,8 @@
"name": "stdout", "name": "stdout",
"output_type": "stream", "output_type": "stream",
"text": [ "text": [
"CPU times: user 82 ms, sys: 4.14 ms, total: 86.1 ms\n", "CPU times: user 151 ms, sys: 360 ms, total: 511 ms\n",
"Wall time: 84.3 ms\n" "Wall time: 57.3 ms\n"
] ]
} }
], ],
...@@ -120,7 +120,7 @@ ...@@ -120,7 +120,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 14, "execution_count": 7,
"id": "69022ad1-5ec8-4a09-836c-273ef452451f", "id": "69022ad1-5ec8-4a09-836c-273ef452451f",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
...@@ -128,7 +128,7 @@ ...@@ -128,7 +128,7 @@
"name": "stdout", "name": "stdout",
"output_type": "stream", "output_type": "stream",
"text": [ "text": [
"188 μs ± 3.45 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)\n" "203 μs ± 4.56 μs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)\n"
] ]
} }
], ],
...@@ -140,7 +140,7 @@ ...@@ -140,7 +140,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 15, "execution_count": 8,
"id": "28031966-c7d0-4201-a467-37590118e851", "id": "28031966-c7d0-4201-a467-37590118e851",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
...@@ -148,7 +148,7 @@ ...@@ -148,7 +148,7 @@
"name": "stdout", "name": "stdout",
"output_type": "stream", "output_type": "stream",
"text": [ "text": [
"8.72 ms ± 74.8 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n" "8.87 ms ± 188 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
] ]
} }
], ],
...@@ -160,7 +160,7 @@ ...@@ -160,7 +160,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 16, "execution_count": 9,
"id": "0d10b2ba-57b2-475c-9d01-135363423990", "id": "0d10b2ba-57b2-475c-9d01-135363423990",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
...@@ -168,8 +168,8 @@ ...@@ -168,8 +168,8 @@
"name": "stdout", "name": "stdout",
"output_type": "stream", "output_type": "stream",
"text": [ "text": [
"CPU times: user 17.3 s, sys: 31.6 ms, total: 17.3 s\n", "CPU times: user 17.4 s, sys: 147 ms, total: 17.6 s\n",
"Wall time: 17.3 s\n" "Wall time: 17.6 s\n"
] ]
} }
], ],
...@@ -188,7 +188,7 @@ ...@@ -188,7 +188,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 17, "execution_count": 10,
"id": "a6aa6926-7fd5-4fff-bd20-7bc0305b948d", "id": "a6aa6926-7fd5-4fff-bd20-7bc0305b948d",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
...@@ -214,7 +214,7 @@ ...@@ -214,7 +214,7 @@
"└──────────┴──────────┴─────────┴────────┴────────┴─────────┴────────┴────────┘" "└──────────┴──────────┴─────────┴────────┴────────┴─────────┴────────┴────────┘"
] ]
}, },
"execution_count": 17, "execution_count": 10,
"metadata": {}, "metadata": {},
"output_type": "execute_result" "output_type": "execute_result"
} }
...@@ -259,7 +259,7 @@ ...@@ -259,7 +259,7 @@
"name": "python", "name": "python",
"nbconvert_exporter": "python", "nbconvert_exporter": "python",
"pygments_lexer": "ipython3", "pygments_lexer": "ipython3",
"version": "3.12.6" "version": "3.12.7"
} }
}, },
"nbformat": 4, "nbformat": 4,
......
notebooks/octtree.ipynb
View file @ 77fb3001
...@@ -133,8 +133,8 @@ ...@@ -133,8 +133,8 @@
"name": "stdout", "name": "stdout",
"output_type": "stream", "output_type": "stream",
"text": [ "text": [
"CPU times: user 186 ms, sys: 191 ms, total: 377 ms\n", "CPU times: user 213 ms, sys: 5.74 ms, total: 219 ms\n",
"Wall time: 118 ms\n" "Wall time: 218 ms\n"
] ]
} }
], ],
...@@ -157,105 +157,105 @@ ...@@ -157,105 +157,105 @@
"output_type": "stream", "output_type": "stream",
"text": [ "text": [
"OctTree:\n", "OctTree:\n",
"- boundary: SpaceTimeRectangle(x = 0, y = 0, w = 360, h = 180, t = 1900-01-16 00:00:00, dt = 32 days, 0:00:00)\n", "- boundary: SpaceTimeRectangle(lon=0, lat=0, date=datetime.datetime(1900, 1, 16, 0, 0), lon_range=360, lat_range=180, dt=datetime.timedelta(days=32))\n",
"- capacity: 10\n", "- capacity: 10\n",
"- depth: 0\n", "- depth: 0\n",
"- max_depth: 25\n", "- max_depth: 25\n",
"- contents:\n", "- contents:\n",
"- number of elements: 10\n", "- number of elements: 10\n",
" * Record(x = 43, y = -68, datetime = 1900-01-08 13:00:00, uid = OBiqSYcn)\n", " * SpaceTimeRecord(x = -12, y = 89, datetime = 1900-01-03 18:00:00, uid = 70Fv8aPq)\n",
" * Record(x = 97, y = -47, datetime = 1900-01-02 14:00:00, uid = w589k3Oe)\n", " * SpaceTimeRecord(x = -157, y = 21, datetime = 1900-01-29 16:00:00, uid = j2xUkgob)\n",
" * Record(x = -68, y = 44, datetime = 1900-01-30 11:00:00, uid = XAaA7McU)\n", " * SpaceTimeRecord(x = -179, y = -25, datetime = 1900-01-30 19:00:00, uid = Nmw6lJSS)\n",
" * Record(x = -170, y = 77, datetime = 1900-01-19 09:00:00, uid = x6eLi65N)\n", " * SpaceTimeRecord(x = -23, y = 89, datetime = 1900-01-14 15:00:00, uid = ILOWgH6u)\n",
" * Record(x = -2, y = 7, datetime = 1900-01-12 09:00:00, uid = CjB2Pglt)\n", " * SpaceTimeRecord(x = 145, y = 41, datetime = 1900-01-26 03:00:00, uid = Cii4cemG)\n",
" * Record(x = -175, y = 65, datetime = 1900-01-15 01:00:00, uid = bTB9DkDI)\n", " * SpaceTimeRecord(x = -139, y = -10, datetime = 1900-01-06 21:00:00, uid = 3QFhOlsW)\n",
" * Record(x = 8, y = 83, datetime = 1900-01-04 10:00:00, uid = aYCKIBl9)\n", " * SpaceTimeRecord(x = -146, y = -51, datetime = 1900-01-18 04:00:00, uid = IHZEXm2l)\n",
" * Record(x = 20, y = 60, datetime = 1900-01-24 16:00:00, uid = 8GsD19WF)\n", " * SpaceTimeRecord(x = -96, y = 86, datetime = 1900-01-14 22:00:00, uid = DtwMgmpH)\n",
" * Record(x = 161, y = 40, datetime = 1900-01-24 20:00:00, uid = FIfAABuC)\n", " * SpaceTimeRecord(x = 113, y = 11, datetime = 1900-01-19 00:00:00, uid = OcFalcj8)\n",
" * Record(x = -69, y = -9, datetime = 1900-01-11 11:00:00, uid = uTcS5D4e)\n", " * SpaceTimeRecord(x = -141, y = -26, datetime = 1900-01-08 11:00:00, uid = qoJquq8j)\n",
"- with children:\n", "- with children:\n",
" OctTree:\n", " OctTree:\n",
" - boundary: SpaceTimeRectangle(x = -90.0, y = 45.0, w = 180.0, h = 90.0, t = 1900-01-08 00:00:00, dt = 16 days, 0:00:00)\n", " - boundary: SpaceTimeRectangle(lon=-90.0, lat=45.0, date=datetime.datetime(1900, 1, 8, 0, 0), lon_range=180.0, lat_range=90.0, dt=datetime.timedelta(days=16))\n",
" - capacity: 10\n", " - capacity: 10\n",
" - depth: 1\n", " - depth: 1\n",
" - max_depth: 25\n", " - max_depth: 25\n",
" - contents:\n", " - contents:\n",
" - number of elements: 10\n", " - number of elements: 10\n",
" * Record(x = -156, y = 57, datetime = 1900-01-08 10:00:00, uid = aFheRU2n)\n", " * SpaceTimeRecord(x = -89, y = 62, datetime = 1900-01-02 15:00:00, uid = zo7Thw1L)\n",
" * Record(x = -100, y = 61, datetime = 1900-01-15 09:00:00, uid = Sa1iavle)\n", " * SpaceTimeRecord(x = -34, y = 62, datetime = 1900-01-13 19:00:00, uid = faiiSA9Y)\n",
" * Record(x = -168, y = 88, datetime = 1900-01-03 07:00:00, uid = IlYKGW0N)\n", " * SpaceTimeRecord(x = -61, y = 54, datetime = 1900-01-07 17:00:00, uid = HwzmXILd)\n",
" * Record(x = -80, y = 50, datetime = 1900-01-05 09:00:00, uid = Rg3GHM4d)\n", " * SpaceTimeRecord(x = -154, y = 21, datetime = 1900-01-08 21:00:00, uid = lA7DANvC)\n",
" * Record(x = -92, y = 39, datetime = 1900-01-15 06:00:00, uid = u804YMFB)\n", " * SpaceTimeRecord(x = -59, y = 87, datetime = 1900-01-02 23:00:00, uid = GSzuXTLF)\n",
" * Record(x = -119, y = 60, datetime = 1900-01-12 22:00:00, uid = vdEPjkib)\n", " * SpaceTimeRecord(x = -109, y = 27, datetime = 1900-01-01 21:00:00, uid = qy7npkH7)\n",
" * Record(x = -160, y = 79, datetime = 1900-01-06 08:00:00, uid = QmrPEL6h)\n", " * SpaceTimeRecord(x = -86, y = 46, datetime = 1900-01-10 19:00:00, uid = RJBRR7Rl)\n",
" * Record(x = -95, y = 21, datetime = 1900-01-09 04:00:00, uid = hfjTKSCH)\n", " * SpaceTimeRecord(x = -114, y = 46, datetime = 1900-01-03 05:00:00, uid = Eop56HgI)\n",
" * Record(x = -93, y = 61, datetime = 1900-01-09 20:00:00, uid = SzIrja9S)\n", " * SpaceTimeRecord(x = -148, y = 49, datetime = 1900-01-08 13:00:00, uid = 3bavQs9B)\n",
" * Record(x = -149, y = 34, datetime = 1900-01-05 05:00:00, uid = b02MxQjV)\n", " * SpaceTimeRecord(x = -109, y = 61, datetime = 1900-01-06 06:00:00, uid = zaTZk8xi)\n",
" - with children:\n", " - with children:\n",
" OctTree:\n", " OctTree:\n",
" - boundary: SpaceTimeRectangle(x = -135.0, y = 67.5, w = 90.0, h = 45.0, t = 1900-01-04 00:00:00, dt = 8 days, 0:00:00)\n", " - boundary: SpaceTimeRectangle(lon=-135.0, lat=67.5, date=datetime.datetime(1900, 1, 4, 0, 0), lon_range=90.0, lat_range=45.0, dt=datetime.timedelta(days=8))\n",
" - capacity: 10\n", " - capacity: 10\n",
" - depth: 2\n", " - depth: 2\n",
" - max_depth: 25\n", " - max_depth: 25\n",
" - contents:\n", " - contents:\n",
" - number of elements: 10\n", " - number of elements: 10\n",
" * Record(x = -134, y = 79, datetime = 1900-01-05 14:00:00, uid = 7Q0FKGMk)\n", " * SpaceTimeRecord(x = -93, y = 49, datetime = 1900-01-07 05:00:00, uid = Vq4kU1vN)\n",
" * Record(x = -90, y = 53, datetime = 1900-01-05 03:00:00, uid = LLx7iz2v)\n", " * SpaceTimeRecord(x = -151, y = 84, datetime = 1900-01-03 03:00:00, uid = 1pz9w5ZR)\n",
" * Record(x = -176, y = 50, datetime = 1900-01-06 20:00:00, uid = x6K5DlTl)\n", " * SpaceTimeRecord(x = -154, y = 74, datetime = 1900-01-03 17:00:00, uid = m1jkOrRF)\n",
" * Record(x = -141, y = 52, datetime = 1900-01-02 15:00:00, uid = xTpGPaEy)\n", " * SpaceTimeRecord(x = -121, y = 81, datetime = 1900-01-02 13:00:00, uid = IcrUGFud)\n",
" * Record(x = -116, y = 68, datetime = 1900-01-05 16:00:00, uid = eECSkpdU)\n", " * SpaceTimeRecord(x = -169, y = 52, datetime = 1900-01-03 12:00:00, uid = Rmytp4VV)\n",
" * Record(x = -138, y = 63, datetime = 1900-01-05 02:00:00, uid = Ftf9uhH3)\n", " * SpaceTimeRecord(x = -115, y = 70, datetime = 1900-01-07 17:00:00, uid = XYiNarG0)\n",
" * Record(x = -173, y = 71, datetime = 1900-01-03 03:00:00, uid = mu3vwHM5)\n", " * SpaceTimeRecord(x = -104, y = 51, datetime = 1900-01-02 06:00:00, uid = je7IfMVs)\n",
" * Record(x = -148, y = 49, datetime = 1900-01-05 15:00:00, uid = 8DFDI3CJ)\n", " * SpaceTimeRecord(x = -161, y = 58, datetime = 1900-01-01 19:00:00, uid = szpJUEjI)\n",
" * Record(x = -157, y = 63, datetime = 1900-01-06 19:00:00, uid = mVqLntgh)\n", " * SpaceTimeRecord(x = -133, y = 46, datetime = 1900-01-03 07:00:00, uid = vD2DhqWZ)\n",
" * Record(x = -154, y = 45, datetime = 1900-01-07 11:00:00, uid = 1UoA1NNC)\n", " * SpaceTimeRecord(x = -125, y = 67, datetime = 1900-01-07 12:00:00, uid = rOJwKqKs)\n",
" - with children:\n", " - with children:\n",
" OctTree:\n", " OctTree:\n",
" - boundary: SpaceTimeRectangle(x = -157.5, y = 78.75, w = 45.0, h = 22.5, t = 1900-01-02 00:00:00, dt = 4 days, 0:00:00)\n", " - boundary: SpaceTimeRectangle(lon=-157.5, lat=78.75, date=datetime.datetime(1900, 1, 2, 0, 0), lon_range=45.0, lat_range=22.5, dt=datetime.timedelta(days=4))\n",
" - capacity: 10\n", " - capacity: 10\n",
" - depth: 3\n", " - depth: 3\n",
" - max_depth: 25\n", " - max_depth: 25\n",
" - contents:\n", " - contents:\n",
" - number of elements: 10\n", " - number of elements: 10\n",
" * Record(x = -147, y = 83, datetime = 1900-01-01 18:00:00, uid = WaO5R7fy)\n", " * SpaceTimeRecord(x = -141, y = 82, datetime = 1900-01-01 13:00:00, uid = B8slJoTY)\n",
" * Record(x = -136, y = 72, datetime = 1900-01-02 03:00:00, uid = OWaMqULr)\n", " * SpaceTimeRecord(x = -147, y = 70, datetime = 1900-01-02 00:00:00, uid = bcnBHvsx)\n",
" * Record(x = -176, y = 79, datetime = 1900-01-02 06:00:00, uid = NTjvqz2c)\n", " * SpaceTimeRecord(x = -170, y = 86, datetime = 1900-01-02 17:00:00, uid = jEUfXIsD)\n",
" * Record(x = -152, y = 72, datetime = 1900-01-03 18:00:00, uid = 7rtQIGtn)\n", " * SpaceTimeRecord(x = -180, y = 78, datetime = 1900-01-02 23:00:00, uid = TYebBlyX)\n",
" * Record(x = -162, y = 78, datetime = 1900-01-02 04:00:00, uid = Wi9RsOIX)\n", " * SpaceTimeRecord(x = -135, y = 68, datetime = 1900-01-01 11:00:00, uid = UtzHGKY0)\n",
" * Record(x = -136, y = 79, datetime = 1900-01-01 11:00:00, uid = hSltzeuH)\n", " * SpaceTimeRecord(x = -136, y = 85, datetime = 1900-01-01 02:00:00, uid = HxwAfFf7)\n",
" * Record(x = -176, y = 89, datetime = 1900-01-02 09:00:00, uid = cOLgAely)\n", " * SpaceTimeRecord(x = -169, y = 71, datetime = 1900-01-03 12:00:00, uid = kGOcDjS4)\n",
" * Record(x = -141, y = 75, datetime = 1900-01-03 23:00:00, uid = gH755dC3)\n", " * SpaceTimeRecord(x = -164, y = 79, datetime = 1900-01-03 11:00:00, uid = IFXekuK1)\n",
" * Record(x = -158, y = 72, datetime = 1900-01-02 23:00:00, uid = NUmMfw9K)\n", " * SpaceTimeRecord(x = -138, y = 76, datetime = 1900-01-01 09:00:00, uid = Qs0TDXtf)\n",
" * Record(x = -168, y = 72, datetime = 1900-01-02 01:00:00, uid = ZFcsxYG4)\n", " * SpaceTimeRecord(x = -154, y = 71, datetime = 1900-01-01 16:00:00, uid = E6eJ2eiF)\n",
" - with children:\n", " - with children:\n",
" OctTree:\n", " OctTree:\n",
" - boundary: SpaceTimeRectangle(x = -168.75, y = 84.375, w = 22.5, h = 11.25, t = 1900-01-01 00:00:00, dt = 2 days, 0:00:00)\n", " - boundary: SpaceTimeRectangle(lon=-168.75, lat=84.375, date=datetime.datetime(1900, 1, 1, 0, 0), lon_range=22.5, lat_range=11.25, dt=datetime.timedelta(days=2))\n",
" - capacity: 10\n", " - capacity: 10\n",
" - depth: 4\n", " - depth: 4\n",
" - max_depth: 25\n", " - max_depth: 25\n",
" - contents:\n",
" - number of elements: 6\n",
" * Record(x = -158, y = 86, datetime = 1900-01-01 15:00:00, uid = DOD5jT2l)\n",
" * Record(x = -165, y = 88, datetime = 1900-01-01 13:00:00, uid = kdGlzz41)\n",
" * Record(x = -173, y = 82, datetime = 1900-01-01 04:00:00, uid = aWBwIP4U)\n",
" * Record(x = -180, y = 89, datetime = 1900-01-01 22:00:00, uid = HOxbaCm8)\n",
" * Record(x = -165, y = 81, datetime = 1900-01-01 16:00:00, uid = JtRn9y9e)\n",
" * Record(x = -164, y = 84, datetime = 1900-01-01 03:00:00, uid = vELpx1ij)\n",
" OctTree:\n", " OctTree:\n",
" - boundary: SpaceTimeRectangle(x = -146.25, y = 84.375, w = 22.5, h = 11.25, t = 1900-01-01 00:00:00, dt = 2 days, 0:00:00)\n", " - boundary: SpaceTimeRectangle(lon=-146.25, lat=84.375, date=datetime.datetime(1900, 1, 1, 0, 0), lon_range=22.5, lat_range=11.25, dt=datetime.timedelta(days=2))\n",
" - capacity: 10\n", " - capacity: 10\n",
" - depth: 4\n", " - depth: 4\n",
" - max_depth: 25\n", " - max_depth: 25\n",
" - contents:\n", " - contents:\n",
" - number of elements: 1\n", " - number of elements: 1\n",
" * Record(x = -157, y = 84, datetime = 1900-01-01 17:00:00, uid = 6DlgVOXg)\n", " * SpaceTimeRecord(x = -135, y = 89, datetime = 1900-01-01 04:00:00, uid = MUrEpv1f)\n",
" OctTree:\n", " OctTree:\n",
" - boundary: SpaceTimeRectangle(x = -168.75, y = 73.125, w = 22.5, h = 11.25, t = 1900-01-01 00:00:00, dt = 2 days, 0:00:00)\n", " - boundary: SpaceTimeRectangle(lon=-168.75, lat=73.125, date=datetime.datetime(1900, 1, 1, 0, 0), lon_range=22.5, lat_range=11.25, dt=datetime.timedelta(days=2))\n",
" - capacity: 10\n", " - capacity: 10\n",
" - depth: 4\n", " - depth: 4\n",
" - max_depth: 25\n", " - max_depth: 25\n",
" - contents:\n", " OctTree:\n",
" - number of elements: 2\n" " - boundary: SpaceTimeRectangle(lon=-146.25, lat=73.125, date=datetime.datetime(1900, 1, 1, 0, 0), lon_range=22.5, lat_range=11.25, dt=datetime.timedelta(days=2))\n",
" - capacity: 10\n",
" - depth: 4\n",
" - max_depth: 25\n",
" OctTree:\n",
" - boundary: SpaceTimeRectangle(lon=-168.75, lat=84.375, date=datetime.datetime(1900, 1, 3, 0, 0), lon_range=22.5, lat_range=11.25, dt=datetime.timedelta(days=2))\n",
" - capacity: 10\n",
" - depth: 4\n",
" - max_depth: 25\n"
] ]
} }
], ],
...@@ -301,7 +301,7 @@ ...@@ -301,7 +301,7 @@
"name": "stdout", "name": "stdout",
"output_type": "stream", "output_type": "stream",
"text": [ "text": [
"206 μs ± 3.36 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)\n" "224 μs ± 6.15 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)\n"
] ]
} }
], ],
...@@ -601,7 +601,7 @@ ...@@ -601,7 +601,7 @@
"name": "stdout", "name": "stdout",
"output_type": "stream", "output_type": "stream",
"text": [ "text": [
"5.33 ms ± 20.1 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n" "5.4 ms ± 19.5 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
] ]
} }
], ],
...@@ -631,7 +631,7 @@ ...@@ -631,7 +631,7 @@
"name": "stdout", "name": "stdout",
"output_type": "stream", "output_type": "stream",
"text": [ "text": [
"CPU times: user 2.52 s, sys: 237 ms, total: 2.75 s\n", "CPU times: user 2.5 s, sys: 228 ms, total: 2.73 s\n",
"Wall time: 2.65 s\n" "Wall time: 2.65 s\n"
] ]
} }
...@@ -699,7 +699,7 @@ ...@@ -699,7 +699,7 @@
"name": "python", "name": "python",
"nbconvert_exporter": "python", "nbconvert_exporter": "python",
"pygments_lexer": "ipython3", "pygments_lexer": "ipython3",
"version": "3.12.6" "version": "3.12.7"
} }
}, },
"nbformat": 4, "nbformat": 4,
......
test/test_kdtree.py
View file @ 77fb3001
import unittest import unittest
import random
from numpy import min, argmin from numpy import min, argmin
from GeoSpatialTools import haversine, KDTree, Record from GeoSpatialTools import haversine, KDTree, Record
...@@ -91,12 +93,39 @@ class TestKDTree(unittest.TestCase): ...@@ -91,12 +93,39 @@ class TestKDTree(unittest.TestCase):
def test_wrap(self): def test_wrap(self):
# TEST: Accounts for wrap at -180, 180 # TEST: Accounts for wrap at -180, 180
kt = KDTree(self.records) kt = KDTree(self.records)
kt.insert(Record(-160, -64, uid="G")) bad_rec = Record(-160, -64, uid="G")
kt.insert(bad_rec)
query_rec = Record(-178, -79, uid="E") query_rec = Record(-178, -79, uid="E")
r, _ = kt.query(query_rec) r, _ = kt.query(query_rec)
assert len(r) == 1 assert len(r) == 1
assert r[0].uid == "C" assert r[0].uid == "C"
def test_near_pole_query(self):
test_records = [
Record(-180, 89.5, uid="1"),
Record(-90, 89.9, uid="2"),
Record(0, 89.5, uid="3"),
]
N_others = 50
test_records.extend(
[
Record(
random.choice(range(-180, 180)),
random.choice(range(80, 90)),
)
for _ in range(N_others)
]
)
kt = KDTree(test_records, max_depth=3)
query_rec = Record(90, 89.8, uid="4")
r, d = kt.query(query_rec)
assert len(r) == 1
print(r[0])
print(d)
assert r[0].uid == "2"
if __name__ == "__main__": if __name__ == "__main__":
unittest.main() unittest.main()
test/test_octtree.py
View file @ 77fb3001
import random
import unittest import unittest
from datetime import datetime, timedelta from datetime import datetime, timedelta
...@@ -27,6 +28,50 @@ class TestRect(unittest.TestCase): ...@@ -27,6 +28,50 @@ class TestRect(unittest.TestCase):
res = list(map(rect.contains, points)) res = list(map(rect.contains, points))
assert res == expected assert res == expected
def test_intersection(self):
d = datetime(2009, 1, 1, 0, 0)
dt = timedelta(days=14)
rect = Rectangle(10, 5, d, 20, 10, dt)
test_rects: list[Rectangle] = [
Rectangle(10, 5, d + timedelta(days=2), 18, 8, dt),
Rectangle(25, 5, d, 9, 12, timedelta(hours=7)),
Rectangle(
15, 8, d - timedelta(hours=18), 12, 7, timedelta(hours=4)
),
Rectangle(15, 8, d + timedelta(days=25), 12, 7, dt),
]
expected = [True, False, True, False]
res = list(map(rect.intersects, test_rects))
assert res == expected
def test_wrap(self):
d = datetime(2009, 1, 1, 0, 0)
dt = timedelta(days=14)
rect = Rectangle(170, 45, d, 180, 20, dt)
assert rect.east < 0
assert rect.west > 0
test_points: list[Record] = [
Record(-140, 40, d),
Record(0, 50, d),
Record(100, 45, d - timedelta(hours=2)),
Record(100, 45, d + timedelta(days=12)),
]
expected = [True, False, True, False]
res = list(map(rect.contains, test_points))
assert res == expected
test_rect = Rectangle(
-100, 40, d + timedelta(days=3), 80, 40, timedelta(days=2)
)
assert test_rect.east < rect.west
assert rect.intersects(test_rect)
# TEST: spatially match, time fail
test_rect = Rectangle(
-100, 40, d + timedelta(days=13), 80, 40, timedelta(days=2)
)
assert not rect.intersects(test_rect)
class TestOctTree(unittest.TestCase): class TestOctTree(unittest.TestCase):
def test_divides(self): def test_divides(self):
...@@ -130,6 +175,36 @@ class TestOctTree(unittest.TestCase): ...@@ -130,6 +175,36 @@ class TestOctTree(unittest.TestCase):
assert res == expected assert res == expected
def test_wrap_query(self):
N = 100
d = datetime(2023, 3, 24, 12, 0)
dt = timedelta(days=10)
boundary = Rectangle(0, 0, d, 360, 180, dt)
ot = OctTree(boundary, capacity=3)
quert_rect = Rectangle(
170, 45, d + timedelta(days=4), 60, 10, timedelta(days=8)
)
points_want: list[Record] = [
Record(175, 43, d + timedelta(days=2)),
Record(-172, 49, d + timedelta(days=4)),
]
points: list[Record] = [
Record(
random.choice(range(-150, 130)),
random.choice(range(-90, 91)),
d + timedelta(hours=random.choice(range(-120, 120))),
)
for _ in range(N)
]
points.extend(points_want)
for p in points:
ot.insert(p)
res = ot.query(quert_rect)
assert len(res) == len(points_want)
assert all([p in res for p in points_want])
def test_ellipse_query(self): def test_ellipse_query(self):
d1 = haversine(0, 2.5, 1, 2.5) d1 = haversine(0, 2.5, 1, 2.5)
d2 = haversine(0, 2.5, 0, 3.0) d2 = haversine(0, 2.5, 0, 3.0)
......
test/test_quadtree.py
View file @ 77fb3001
from math import pi import random
import unittest import unittest
from GeoSpatialTools import haversine from GeoSpatialTools import haversine
from GeoSpatialTools.quadtree import QuadTree, Record, Rectangle, Ellipse from GeoSpatialTools.quadtree import QuadTree, Record, Rectangle, Ellipse
...@@ -29,6 +30,23 @@ class TestRect(unittest.TestCase): ...@@ -29,6 +30,23 @@ class TestRect(unittest.TestCase):
res = list(map(rect.intersects, test_rects)) res = list(map(rect.intersects, test_rects))
assert res == expected assert res == expected
def test_wrap(self):
rect = Rectangle(170, 45, 180, 20)
assert rect.east < 0
assert rect.west > 0
test_points: list[Record] = [
Record(-140, 40),
Record(0, 50),
Record(100, 45),
]
expected = [True, False, True]
res = list(map(rect.contains, test_points))
assert res == expected
test_rect = Rectangle(-100, 40, 80, 40)
assert test_rect.east < rect.west
assert rect.intersects(test_rect)
class TestQuadTree(unittest.TestCase): class TestQuadTree(unittest.TestCase):
def test_divides(self): def test_divides(self):
...@@ -98,6 +116,31 @@ class TestQuadTree(unittest.TestCase): ...@@ -98,6 +116,31 @@ class TestQuadTree(unittest.TestCase):
assert res == expected assert res == expected
def test_wrap_query(self):
N = 100
qt_boundary = Rectangle(0, 0, 360, 180)
qt = QuadTree(qt_boundary, capacity=3)
quert_rect = Rectangle(170, 45, 60, 10)
points_want: list[Record] = [
Record(175, 43),
Record(-172, 49),
]
points: list[Record] = [
Record(
random.choice(range(-150, 130)),
random.choice(range(-90, 91)),
)
for _ in range(N)
]
points.extend(points_want)
for p in points:
qt.insert(p)
res = qt.query(quert_rect)
assert len(res) == len(points_want)
assert all([p in res for p in points_want])
def test_ellipse_query(self): def test_ellipse_query(self):
d1 = haversine(0, 2.5, 1, 2.5) d1 = haversine(0, 2.5, 1, 2.5)
d2 = haversine(0, 2.5, 0, 3.0) d2 = haversine(0, 2.5, 0, 3.0)
......
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