GeoSpatialTools

Python module containing useful functions and classes for Spatial Analysis.

Tested on Python versions 3.11, 3.12, and 3.13. Not expected to work with Python 3.9 or below due to usage of type annotations.

Installation

As a dependency with pip

pip install git+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:

def find_nearest(vals: list[Numeric], test: list[Numeric]) -> list[int]:
    ...

Example:

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:

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).

KDTree(
    points: list[Record],
    max_depth: int,
    depth: int = 0  # Internal parameter
)

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.

QuadTree(
    boundary: Rectangle,
    capacity: int,
    max_depth: int,
    depth: int = 0  # Internal parameter
)

Rectangle(
    lon: float,  # Centre longitude
    lat: float,  # Centre latitude
    lon_range: float,  # Full width of rectangle (degrees)
    lat_range: float,  # Full height of rectangle (degrees)
)

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(0, 0, 360, 180)  # 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.

SpaceTimeRecord(
    lon: float,
    lat: float,
    datetime: datetime.datetime,  # datetime no longer optional
    uid: str | None,
    **data
)

SpaceTimeRectangle(
    lon: float,          # Centre longitude
    lat: float,          # Centre latitude
    datetime: datetime,  # Central datetime
    w: float,            # Full width of rectangle (degrees)
    h: float,            # Full height of rectangle (degrees)
    dt: timedelta,       # Time extent of rectangle
)

Example

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(0, 0, datetime(2009, 1, 15, 12), 360, 180, timedelta(days=31))  # 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)