from datetime import datetime, timedelta
from .distance_metrics import haversine
from .distance_metrics import haversine, destination
from math import degrees, sqrt
class SpaceTimeRecord:
"""
ICOADS Record class.
This is a simple instance of an ICOARDS record, it requires position and
temporal data. It can optionally include a UID and extra data.
The temporal component was designed to use `datetime` values, however all
methods will work with numeric datetime information - for example a pentad,
timestamp, julian day, etc. Note that any uses within an OctTree and
SpaceTimeRectangle must also have timedelta values replaced with numeric
ranges in this case.
Equality is checked only on the required fields + UID if it is specified.
Parameters
----------
lon : float
Horizontal coordinate (longitude).
lat : float
Vertical coordinate (latitude).
datetime : datetime
Datetime of the record. Can also be a numeric value such as pentad.
Comparisons between Records with datetime and Records with numeric
datetime will fail.
uid : str | None
Unique Identifier.
**data
Additional data passed to the SpaceTimeRecord for use by other functions
or classes.
"""
def __init__(
self,
lon: float,
lat: float,
datetime: datetime,
uid: str | None = None,
**data,
) -> None:
self.lon = lon
self.lat = lat
self.datetime = datetime
self.uid = uid
for var, val in data.items():
setattr(self, var, val)
return None
def __str__(self) -> str:
return f"Record(x = {self.lon}, y = {self.lat}, datetime = {self.datetime}, uid = {self.uid})"
def __eq__(self, other: object) -> bool:
return (
isinstance(other, SpaceTimeRecord)
and self.lon == other.lon
and self.lat == other.lat
and self.datetime == other.datetime
and (not (self.uid or other.uid) or self.uid == other.uid)
)
class SpaceTimeRecords(list[SpaceTimeRecord]):
"""List of SpaceTimeRecords"""
class SpaceTimeRectangle:
"""
A simple Space Time SpaceTimeRectangle class.
This constructs a simple Rectangle object.
The defining coordinates are the centres of the box, and the extents
are the full width, height, and time extent.
Whilst the rectangle is assumed to lie on the surface of Earth, this is
a projection as the rectangle is defined by a longitude/latitude range.
The temporal components are defined in the same way as the spatial
components, that is that the `datetime` component (t) is the "centre", and
the time extent (dt) is the full time range of the box.
Parameters
----------
lon : float
Horizontal centre of the rectangle (longitude).
lat : float
Vertical centre of the rectangle (latitude).
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
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
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})"
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
)
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
)
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
)
def nearby(
self,
point: SpaceTimeRecord,
dist: float,
t_dist: timedelta,
) -> bool:
"""
Check if point is nearby the Rectangle
Determines if a SpaceTimeRecord that falls on the surface of Earth is
nearby to the rectangle in space and time. This calculation uses the
Haversine distance metric.
Distance from rectangle to point is challenging on the surface of a
sphere, this calculation will return false positives as a check based
on the distance from the centre of the rectangle to the corners, or
to its Eastern edge (if the rectangle crosses the equator) is used in
combination with the input distance.
The primary use-case of this method is for querying an OctTree for
nearby Records.
Parameters
----------
point : SpaceTimeRecord
dist : float,
t_dist : 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
):
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
)
class SpaceTimeEllipse:
"""
A simple Ellipse Class for an ellipse on the surface of a sphere.
Parameters
----------
lon : float
Horizontal centre of the ellipse
lat : float
Vertical centre of the ellipse
datetime : datetime
Datetime centre of the ellipse.
a : float
Length of the semi-major axis
b : float
Length of the semi-minor axis
theta : float
Angle of the semi-major axis from horizontal anti-clockwise in radians
dt : timedelta
(full) time extent of the ellipse.
"""
def __init__(
self,
lon: float,
lat: float,
datetime: datetime,
a: float,
b: float,
theta: float,
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
)
class OctTree:
"""
A Simple OctTree class for PyCOADS.
Acts as a space-time OctTree on the surface of Earth, allowing for querying
nearby points faster than searching a full DataFrame. As SpaceTimeRecords
are added to the OctTree, the OctTree divides into 8 children as the
capacity is reached. Additional SpaceTimeRecords are then added to the
children where they fall within the child OctTree's boundary.
SpaceTimeRecords already part of the OctTree before divided are not
distributed to the children OctTrees.
Whilst the OctTree has a temporal component, and was designed to utilise
datetime / timedelta objects, numeric values and ranges can be used. This
usage must be consistent for the boundary and all SpaceTimeRecords that
are part of the OctTree. This allows for usage of pentad, timestamp,
Julian day, etc. as datetime values.
Parameters
----------
boundary : SpaceTimeRectangle
The bounding SpaceTimeRectangle of the QuadTree
capacity : int
The capacity of each cell, if max_depth is set then a cell at the
maximum depth may contain more points than the capacity.
depth : int
The current depth of the cell. Initialises to zero if unset.
max_depth : int | None
The maximum depth of the QuadTree. If set, this can override the
capacity for cells at the maximum depth.
"""
def __init__(
self,
boundary: SpaceTimeRectangle,
capacity: int = 5,
depth: int = 0,
max_depth: int | None = None,
) -> None:
self.boundary = boundary
self.capacity = capacity
self.depth = depth
self.max_depth = max_depth
self.points = SpaceTimeRecords()
self.divided: bool = False
return None
def __str__(self) -> str:
indent = " " * self.depth
out = f"{indent}OctTree:\n"
out += f"{indent}- boundary: {self.boundary}\n"
out += f"{indent}- capacity: {self.capacity}\n"
out += f"{indent}- depth: {self.depth}\n"
if self.max_depth:
out += f"{indent}- max_depth: {self.max_depth}\n"
if self.points:
out += f"{indent}- contents:\n"
out += f"{indent}- number of elements: {len(self.points)}\n"
for p in self.points:
out += f"{indent} * {p}\n"
if self.divided:
out += f"{indent}- with children:\n"
out += f"{self.northwestback}"
out += f"{self.northeastback}"
out += f"{self.southwestback}"
out += f"{self.southeastback}"
out += f"{self.northwestfwd}"
out += f"{self.northeastfwd}"
out += f"{self.southwestfwd}"
out += f"{self.southeastfwd}"
return out
def divide(self):
"""Divide the QuadTree"""
self.northwestfwd = 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.lon_range / 2,
self.boundary.lat_range / 2,
self.boundary.dt / 2,
),
capacity=self.capacity,
depth=self.depth + 1,
max_depth=self.max_depth,
)
self.northeastfwd = 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.lon_range / 2,
self.boundary.lat_range / 2,
self.boundary.dt / 2,
),
capacity=self.capacity,
depth=self.depth + 1,
max_depth=self.max_depth,
)
self.southwestfwd = 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.lon_range / 2,
self.boundary.lat_range / 2,
self.boundary.dt / 2,
),
capacity=self.capacity,
depth=self.depth + 1,
max_depth=self.max_depth,
)
self.southeastfwd = 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.lon_range / 2,
self.boundary.lat_range / 2,
self.boundary.dt / 2,
),
capacity=self.capacity,
depth=self.depth + 1,
max_depth=self.max_depth,
)
self.northwestback = 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.lon_range / 2,
self.boundary.lat_range / 2,
self.boundary.dt / 2,
),
capacity=self.capacity,
depth=self.depth + 1,
max_depth=self.max_depth,
)
self.northeastback = 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.lon_range / 2,
self.boundary.lat_range / 2,
self.boundary.dt / 2,
),
capacity=self.capacity,
depth=self.depth + 1,
max_depth=self.max_depth,
)
self.southwestback = 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.lon_range / 2,
self.boundary.lat_range / 2,
self.boundary.dt / 2,
),
capacity=self.capacity,
depth=self.depth + 1,
max_depth=self.max_depth,
)
self.southeastback = 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.lon_range / 2,
self.boundary.lat_range / 2,
self.boundary.dt / 2,
),
capacity=self.capacity,
depth=self.depth + 1,
max_depth=self.max_depth,
)
self.divided = True
def _datetime_is_numeric(self) -> bool:
return not isinstance(self.boundary.datetime, datetime)
def insert(self, point: SpaceTimeRecord) -> bool:
"""
Insert a SpaceTimeRecord into the QuadTree.
Note that the SpaceTimeRecord can have numeric datetime values if that
is consistent with the OctTree.
"""
if not self.boundary.contains(point):
return False
elif self.max_depth and self.depth == self.max_depth:
self.points.append(point)
return True
elif len(self.points) < self.capacity:
self.points.append(point)
return True
else:
if not self.divided:
self.divide()
if self.northwestback.insert(point):
return True
elif self.northeastback.insert(point):
return True
elif self.southwestback.insert(point):
return True
elif self.southeastback.insert(point):
return True
elif self.northwestfwd.insert(point):
return True
elif self.northeastfwd.insert(point):
return True
elif self.southwestfwd.insert(point):
return True
elif self.southeastfwd.insert(point):
return True
return False
def query(
self,
rect: SpaceTimeRectangle,
points: SpaceTimeRecords | None = None,
) -> SpaceTimeRecords:
"""Get points that fall in a SpaceTimeRectangle"""
if not points:
points = SpaceTimeRecords()
if not self.boundary.intersects(rect):
return points
for point in self.points:
if rect.contains(point):
points.append(point)
if self.divided:
points = self.northwestfwd.query(rect, points)
points = self.northeastfwd.query(rect, points)
points = self.southwestfwd.query(rect, points)
points = self.southeastfwd.query(rect, points)
points = self.northwestback.query(rect, points)
points = self.northeastback.query(rect, points)
points = self.southwestback.query(rect, points)
points = self.southeastback.query(rect, points)
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,
dist: float,
t_dist: timedelta,
points: SpaceTimeRecords | None = None,
) -> SpaceTimeRecords:
"""
Get all points that are nearby another point.
Query the OctTree to find all SpaceTimeRecords within the OctTree that
are nearby to the query SpaceTimeRecord. This search should be faster
than searching through all records, since only OctTree children whose
boundaries are close to the query SpaceTimeRecord are evaluated.
Parameters
----------
point : SpaceTimeRecord
The query point.
dist : float
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
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
values and ranges.
points : SpaceTimeRecords | None
List of SpaceTimeRecords already found. Most use cases will be to
not set this value, since it's main use is for passing onto the
children OctTrees.
Returns
-------
SpaceTimeRecords : A list of SpaceTimeRecords whose distance to the
query SpaceTimeRecord is <= dist, and the datetimes of the
SpaceTimeRecords fall within the datetime range of the query
SpaceTimeRecord.
"""
if not points:
points = SpaceTimeRecords()
if not self.boundary.nearby(point, dist, t_dist):
return points
for test_point in self.points:
if (
haversine(point.lon, point.lat, test_point.lon, test_point.lat)
<= dist
and test_point.datetime <= point.datetime + t_dist
and test_point.datetime >= point.datetime - t_dist
):
points.append(test_point)
if self.divided:
points = self.northwestback.nearby_points(
point, dist, t_dist, points
)
points = self.northeastback.nearby_points(
point, dist, t_dist, points
)
points = self.southwestback.nearby_points(
point, dist, t_dist, points
)
points = self.southeastback.nearby_points(
point, dist, t_dist, points
)
points = self.northwestfwd.nearby_points(
point, dist, t_dist, points
)
points = self.northeastfwd.nearby_points(
point, dist, t_dist, points
)
points = self.southwestfwd.nearby_points(
point, dist, t_dist, points
)
points = self.southeastfwd.nearby_points(
point, dist, t_dist, points
)
return points