# GeoSpatialTools Python module containing useful functions and classes for Spatial Analysis. Tested on Python versions 3.9 to 3.13. ## Installation As a dependency with `pip` ```bash pip install git+ssh://git@git.noc.ac.uk/nocsurfaceprocesses/geospatialtools.git ``` ## Neighbours ### 1d neighbours For example, finding the closest time value in a list of _sorted_ time values: ```python def find_nearest(vals: list[Numeric], test: list[Numeric]) -> list[int]: ... ``` Example: ```python from GeoSpatialTools import find_nearest import numpy as np # Generate random neighbours and test points potential_neighbours = list(np.random.randn(100)) potential_neighbours.sort() test_values = list(np.random.randn(3)) neighbour_i = find_nearest(potential_neighbours, test_values) ``` Returns the index of the value in `potential_neighbours` that is closest to each value in `test_values`. This function makes use of `bisect.bisect` from the python standard library and requires that the list `potential_neighbours` is sorted. No check on sorted-ness is performed, if the list is unsorted then the incorrect answer will be returned. This calculation is performed in `O(log(n))` time. ### 2d neighbours An implementation of `KDTree` (specifically a `2DTree`) using Haversine distance as the query distance and accounting for longitudinal wrapping at -180, 180. Uses a `Record` class: ```python Record( lon: float, lat: float, datetime: datetime.datetime | None, uid: str | None, **data ) ``` Fast look-up of nearest neighbours in 2 spatial dimensions (longitude and latitude). ```python KDTree( points: list[Record], max_depth: int, depth: int = 0 # Internal parameter ) ``` ```python from GeoSpatialTools import KDTree, Record from random import choice lon_range = list(range(-180, 180)) lat_range = list(range(-90, 90)) N_samples = 1000 records: list[Record] = [Record(choice(lon_range), choice(lat_range)) for _ in range(N_samples)] # Construct Tree kt = KDTree(records) test_value: Record = Record(lon=47.6, lat=-31.1) neighbours: list[Record] = [] neighbours, dist = kt.query(test_value) ``` ### Points within distance (2d \& 3d) Also included is an implementation of `QuadTree` and `OctTree` using Haversine distance for querying with a test `Record` and distance. ```python QuadTree( boundary: Rectangle, capacity: int, max_depth: int, depth: int = 0 # Internal parameter ) ``` The boundary is defined by a `Rectangle` class, which is defined by the bounding box. ```python Rectangle( west: float, # Western edge east: float, # Eastern edge south: float, # Southern edge north: float, # Northern edge ) ``` The `Rectangle` class will raise an error if the northern or southern boundary go beyond the north or south pole. ```python from GeoSpatialTools import QuadTree, Record, Rectangle from random import choice lon_range = list(range(-180, 180)) lat_range = list(range(-90, 90)) N_samples = 1000 # Construct Tree boundary = Rectangle(-180, 180, -90, 90) # Full domain qt = QuadTree(boundary) records: list[Record] = [Record(choice(lon_range), choice(lat_range)) for _ in range(N_samples)] for record in records: qt.insert(record) test_value: Record = Record(lon=47.6, lat=-31.1) dist: float = 340 # km neighbours: list[Record] = qt.nearby_points(test_value, dist) ``` #### OctTree - 3d QuadTree Adds `SpaceTimeRecord`, `SpaceTimeRectangle` and `OctTree` classes. This allows for querying in a third dimension, specifically this adds a time dimension. Typically the time dimension is assumed to be of `datetime.datetime` type, however this is expected to work with numeric values, for example pentad, day of year. However, this non-datetime behaviour is not intended. ```python SpaceTimeRecord( lon: float, lat: float, datetime: datetime.datetime, # datetime no longer optional uid: str | None, **data ) ``` As with the `Rectangle` class for the `QuadTree`, the `SpaceTimeRectangle` defines the boundary of an `OctTree` class, and is defined by the space-time bounding box. ```python SpaceTimeRectangle( west: float, # Western edge east: float, # Eastern edge south: float, # Southern edge north: float, # Northern edge start: datetime.datetime, # Start datetime end: datetime.datetime, # End datetime ) ``` Example ```python from GeoSpatialTools.octtree import OctTree, SpaceTimeRecord, SpaceTimeRectangle from datetime import datetime, timedelta from random import choice from pandas import date_range lon_range = list(range(-180, 180)) lat_range = list(range(-90, 90)) dates = date_range( start=datetime(2009, 1, 1, 0, 0), end=datetime(2009, 2, 1, 0, 0), interval=timedelta(hours=1), inclusive="left", ) N_samples = 1000 # Construct Tree boundary = SpaceTimeRectangle(-180, 180, -90, 90, datetime(2009, 1, 1, 0), datetime(2009, 1, 2, 23)) # Full domain ot = OctTree(boundary) records: list[SpaceTimeRecord] = [ SpaceTimeRecord(choice(lon_range), choice(lat_range), choice(dates)) for _ in range(N_samples)] for record in records: ot.insert(record) test_value: SpaceTimeRecord = SpaceTimeRecord(lon=47.6, lat=-31.1, datetime=datetime(2009, 1, 23, 17, 41)) dist: float = 340 # km t_dist = timedelta(hours=4) neighbours: list[Record] = ot.nearby_points(test_value, dist, t_dist) ```