Math and Algorithms Reference for Geospatial Software Engineers¶

This reference collects the computation patterns that appear throughout the book. The goal is to make the math visible enough for software engineers to reason about implementation, testing, and failure modes.

1. Coordinate Tuples and Vectors¶

A 2D point is usually represented as:

p = [x, y]

A 3D point is:

p = [x, y, z]

A spatiotemporal observation is:

o = [x, y, z, t, attributes...]

Use this for points, trajectories, point clouds, 3D scenes, space-time cubes, and moving-object systems.

2. Euclidean Distance¶

For projected planar coordinates:

d(p1, p2) = sqrt((x2 - x1)^2 + (y2 - y1)^2)

For 3D:

d(p1, p2) = sqrt((x2 - x1)^2 + (y2 - y1)^2 + (z2 - z1)^2)

Use only when coordinates are in an appropriate projected coordinate system with linear units.

3. Haversine Great-Circle Distance¶

For approximate distance on a sphere:

a = sin^2((lat2 - lat1) / 2)
    + cos(lat1) * cos(lat2) * sin^2((lon2 - lon1) / 2)

c = 2 * atan2(sqrt(a), sqrt(1 - a))

d = R * c

Where:

lat and lon are in radians.
R is Earth radius.
d is the great-circle distance.

For high-accuracy geodesic work, use a geodesic library instead of implementing this yourself.

4. Affine Transformation¶

Affine transformations support translation, rotation, scaling, and shearing:

[x']   [a b tx] [x]
[y'] = [c d ty] [y]
[1 ]   [0 0 1 ] [1]

Use this for image georeferencing, map rendering, screen-to-map transforms, and simple coordinate conversions inside one reference frame.

5. Raster Geotransform¶

A common raster geotransform maps pixel coordinates to map coordinates:

X_map = GT0 + col * GT1 + row * GT2
Y_map = GT3 + col * GT4 + row * GT5

For north-up rasters, GT2 and GT4 are often 0, GT1 is pixel width, and GT5 is negative pixel height.

6. Polygon Area: Shoelace Formula¶

For a simple planar polygon:

area = 0.5 * abs(sum(i=1..n) (x_i * y_(i+1) - x_(i+1) * y_i))

The ring wraps so point n + 1 is point 1.

Use this only for planar coordinates. For geographic coordinates, use geodesic area or project first.

7. Point-in-Polygon: Ray Crossing Concept¶

A common algorithm casts a ray from a point and counts boundary crossings:

inside = (number_of_ray_crossings % 2) == 1

Boundary cases are tricky. In production, use a tested geometry library such as GEOS, JTS, Shapely, or PostGIS.

8. Line Segment Intersection: Orientation Test¶

For points a, b, and c:

orient(a, b, c) = sign((b_x - a_x) * (c_y - a_y)
                      - (b_y - a_y) * (c_x - a_x))

Two line segments usually intersect when endpoint orientations differ across both segments, with special handling for collinear cases.

9. Bounding Box Filter¶

A bounding box is:

bbox = [minx, miny, maxx, maxy]

Two bounding boxes intersect when:

intersects =
  bbox1.minx <= bbox2.maxx and
  bbox1.maxx >= bbox2.minx and
  bbox1.miny <= bbox2.maxy and
  bbox1.maxy >= bbox2.miny

Spatial indexes often use bounding boxes as fast candidate filters before exact geometry predicates.

10. Buffer Concept¶

A buffer around a geometry is the set of all points within distance r:

buffer(G, r) = { p | distance(p, G) <= r }

Correctness depends on CRS, distance metric, and geometry validity.

11. Graph Shortest Path¶

A weighted graph is:

G = (V, E, w)

Where:

V is vertices or nodes.
E is edges.
w(e) is edge cost.

Shortest path objective:

minimize sum(e in path) w(e)

Use for roads, pipes, utility traces, indoor routing, river networks, and graph-based GeoAI.

12. Raster Map Algebra¶

Raster map algebra computes cell-by-cell:

output[row, col] = f(raster1[row, col], raster2[row, col], ...)

Example normalized difference vegetation index:

NDVI = (NIR - Red) / (NIR + Red)

Always handle nodata, scale factors, band meaning, and division by zero.

13. Bilinear Interpolation¶

For a value inside a raster cell neighborhood:

f(x, y) =
  f00 * (1 - dx) * (1 - dy) +
  f10 * dx       * (1 - dy) +
  f01 * (1 - dx) * dy       +
  f11 * dx       * dy

Use for continuous surfaces such as elevation or temperature. Avoid for categorical classes such as land cover.

14. Inverse Distance Weighting¶

For interpolation from nearby samples:

value(x) = sum(i=1..n) w_i * value_i / sum(i=1..n) w_i

w_i = 1 / distance(x, x_i)^p

The power p controls how quickly influence decays with distance.

15. Spatial Autocorrelation: Moran's I¶

Moran's I measures whether nearby values are more similar or dissimilar than expected:

I = (n / W) *
    (sum_i sum_j w_ij * (x_i - x_bar) * (x_j - x_bar)) /
    (sum_i (x_i - x_bar)^2)

Where:

w_ij is the spatial weight between locations i and j.
W = sum_i sum_j w_ij.
x_bar is the mean.

This is foundational for GeoAI leakage, clustering, spatial statistics, and model validation.

16. Space-Time Cube¶

A space-time cube bins observations by location and time:

cube[cell_id, time_bin] = aggregate(events where
                                   event.location in cell_id and
                                   event.time in time_bin)

Common aggregations:

count
mean(value)
sum(value)
max(value)
trend(value over time)

Conceptually:

X = longitude or projected x
Y = latitude or projected y
T = time
cube = X by Y by T

For 3D or atmospheric systems:

cube = X by Y by Z by T

17. Time Series Lag and Difference¶

Lag:

lag_k(x_t) = x_(t-k)

First difference:

delta_t = x_t - x_(t-1)

Percent change:

pct_change_t = (x_t - x_(t-1)) / x_(t-1)

Use in space-time cubes, change detection, anomaly detection, mobility, weather, sensor streams, and Earth observation time series.

18. Linear Algebra for Embeddings¶

A vector embedding is:

v = [v1, v2, ..., vd]

Dot product:

dot(a, b) = sum(i=1..d) a_i * b_i

Vector norm:

||a|| = sqrt(sum(i=1..d) a_i^2)

Cosine similarity:

cosine_similarity(a, b) = dot(a, b) / (||a|| * ||b||)

This is foundational to vector search, semantic retrieval, GeoRAG, image embeddings, spatial embeddings, and hybrid search.

19. Hybrid Spatial and Semantic Retrieval¶

GeoRAG often combines text relevance, vector similarity, spatial filters, and time filters:

score(document) =
  alpha * semantic_similarity(query_embedding, document_embedding)
  + beta * text_score(query, document)
  + gamma * spatial_relevance(query_geometry, document_geometry)
  + delta * temporal_relevance(query_time, document_time)

Candidate filter:

candidate if
  ST_Intersects(document.geometry, query_area)
  and document.time overlaps query_time_window

This is the mathematical heart of GeoRAG: retrieve by meaning, place, and time.

20. Attention in Transformer Models¶

Scaled dot-product attention:

Attention(Q, K, V) = softmax((Q * K^T) / sqrt(d_k)) * V

Where:

Q is queries.
K is keys.
V is values.
d_k is key dimension.

GeoRAG does not replace this model math. It adds geospatial retrieval and grounding around the model.

21. Spatial Train-Test Split¶

A random split can leak nearby information. A spatial block split assigns geography to folds:

fold_id = spatial_block(location)

train = observations where fold_id != held_out_fold
test  = observations where fold_id == held_out_fold

This better estimates performance in new areas.

22. Error and Residuals¶

Prediction residual:

residual_i = observed_i - predicted_i

Root mean squared error:

RMSE = sqrt(mean(residual_i^2))

Map residuals. A model with acceptable global RMSE can still fail badly in particular regions.

23. Optimization and Suitability¶

Weighted suitability score:

score(location) = sum(i=1..n) weight_i * normalized_criterion_i(location)

Constraint:

eligible(location) =
  not in_exclusion_zone(location)
  and slope(location) <= max_slope
  and distance_to_road(location) <= max_distance

Use for site selection, conservation planning, risk screening, and planning models.

24. Testing Spatial Math in Software¶

Every equation should have known-answer tests:

given input geometry/data
when operation runs
then output equals expected result within tolerance

Example tolerance check:

abs(actual - expected) <= epsilon

For geometry:

area(symmetric_difference(actual, expected)) <= epsilon

25. 3D Gaussian Splatting for Field Inspection¶

3D Gaussian Splatting represents a scene with many anisotropic Gaussian primitives. Each primitive has a position, covariance, opacity, and color or learned radiance parameters:

g_i = {mu_i, Sigma_i, alpha_i, color_i}

mu_i = [x_i, y_i, z_i]
Sigma_i = R_i * S_i * S_i^T * R_i^T

The Gaussian density around a primitive is:

G_i(x) = exp(-0.5 * (x - mu_i)^T * inverse(Sigma_i) * (x - mu_i))

A camera projects a 3D point into image space with camera intrinsics and pose:

u = K * (R * X + t)

Where:

X is a 3D point in world coordinates.
R and t are camera rotation and translation.
K is the camera intrinsic matrix.
u is the projected image coordinate.

Rendered color along a ray uses front-to-back alpha compositing:

C = sum(i=1..n) T_i * alpha_i * color_i
T_i = product(j=1..i-1) (1 - alpha_j)

For inspection, compare a baseline scene and a later scene using geometry, image, or semantic residuals:

inspection_change_score(asset) =
  w_geometry * geometric_residual(asset)
  + w_image * photometric_residual(asset)
  + w_semantic * defect_probability(asset)

Use Gaussian Splatting as visual evidence and inspection context. For authoritative measurement, tie the scene to control points, calibrated camera poses, LiDAR, photogrammetry outputs, or survey-grade measurements.