GIS Emergency Response Network: Multi-Step Spatial Analytics

NYC dispatch scenario combining st_distance nearest-neighbor ranking, multi-algorithm distance comparison, navigation bearings, and WKT polygon geofencing across 8 hospitals, 12 incidents, and 4 response zones.

Category: spatial

Syntax

-- ============================================================================
-- SETUP, zone, schema, 3 Delta tables (hospitals, incidents, response_zones)
-- ============================================================================
CREATE ZONE IF NOT EXISTS {{zone_name}} TYPE DELTA
    COMMENT 'Delta tables for demo datasets';

CREATE SCHEMA IF NOT EXISTS {{zone_name}}.emergency
    COMMENT 'NYC emergency response, hospitals, incidents, response zones';

CREATE DELTA TABLE IF NOT EXISTS {{zone_name}}.emergency.hospitals (
    hospital_id INT, hospital_name VARCHAR, lat DOUBLE, lng DOUBLE,
    trauma_level INT, bed_capacity INT
) LOCATION '{{data_path}}/emergency_hospitals';

CREATE DELTA TABLE IF NOT EXISTS {{zone_name}}.emergency.incidents (
    incident_id INT, incident_type VARCHAR, lat DOUBLE, lng DOUBLE,
    severity VARCHAR, reported_at VARCHAR
) LOCATION '{{data_path}}/emergency_incidents';

CREATE DELTA TABLE IF NOT EXISTS {{zone_name}}.emergency.response_zones (
    zone_id INT, zone_name VARCHAR, zone_polygon VARCHAR, priority_level INT
) LOCATION '{{data_path}}/emergency_zones';

-- ============================================================================
-- QUERIES, layered GIS analytics
-- ============================================================================
-- Nearest hospital per incident (ROW_NUMBER + st_distance)
WITH ranked AS (
    SELECT i.incident_id, i.incident_type, i.severity, h.hospital_name,
           ROUND(st_distance(i.lat, i.lng, h.lat, h.lng) / 1000.0, 1) AS distance_km,
           ROW_NUMBER() OVER (PARTITION BY i.incident_id
                              ORDER BY st_distance(i.lat, i.lng, h.lat, h.lng)) AS rn
    FROM {{zone_name}}.emergency.incidents i
    CROSS JOIN {{zone_name}}.emergency.hospitals h
)
SELECT incident_id, incident_type, severity,
       hospital_name AS nearest_hospital, distance_km
FROM ranked WHERE rn = 1 ORDER BY incident_id;

-- Distance algorithm comparison (haversine vs sphere vs Vincenty)
SELECT ROUND(st_distance_haversine(h1.lat, h1.lng, h2.lat, h2.lng) / 1000.0, 1) AS haversine_km,
       ROUND(st_distance_sphere  (h1.lat, h1.lng, h2.lat, h2.lng) / 1000.0, 1) AS sphere_km,
       ROUND(st_distance_vincenty(h1.lat, h1.lng, h2.lat, h2.lng) / 1000.0, 1) AS vincenty_km
FROM {{zone_name}}.emergency.hospitals h1 CROSS JOIN {{zone_name}}.emergency.hospitals h2
WHERE h1.hospital_id = 8 AND h2.hospital_id = 4;

-- Dispatch bearings from Bellevue to every critical incident
SELECT i.incident_id, i.incident_type,
       ROUND(st_bearing      (h.lat, h.lng, i.lat, i.lng), 2) AS bearing_deg,
       ROUND(st_azimuth      (h.lat, h.lng, i.lat, i.lng), 4) AS azimuth_rad,
       ROUND(st_final_bearing(h.lat, h.lng, i.lat, i.lng), 2) AS final_bearing_deg
FROM {{zone_name}}.emergency.incidents i
CROSS JOIN {{zone_name}}.emergency.hospitals h
WHERE h.hospital_id = 1 AND i.severity = 'critical'
ORDER BY i.incident_id;

-- Geofencing, st_contains(polygon, lat, lng)
SELECT i.incident_id, i.incident_type, i.severity, z.zone_name, z.priority_level
FROM {{zone_name}}.emergency.incidents i
CROSS JOIN {{zone_name}}.emergency.response_zones z
WHERE st_contains(z.zone_polygon, i.lat, i.lng) = true
ORDER BY i.incident_id;

-- Coordinate round-trip, st_make_point(lng, lat) then st_x / st_y
SELECT hospital_id, hospital_name,
       st_make_point(lng, lat)                    AS wkt_point,
       st_x(st_make_point(lng, lat))              AS extracted_lng,
       st_y(st_make_point(lng, lat))              AS extracted_lat
FROM {{zone_name}}.emergency.hospitals ORDER BY hospital_id;

-- Geodesic zone area, st_area returns m^2
SELECT MAX(CASE WHEN zone_id = 1 THEN ROUND(st_area(zone_polygon)/1e6, 2) END) AS area_km2_downtown,
       MAX(CASE WHEN zone_id = 3 THEN ROUND(st_area(zone_polygon)/1e6, 2) END) AS area_km2_brooklyn
FROM {{zone_name}}.emergency.response_zones;

-- Nearest Level-1 trauma center per critical incident (distance + bearing)
WITH trauma1_ranked AS (
    SELECT i.incident_id, i.incident_type, h.hospital_name,
           ROUND(st_distance(i.lat, i.lng, h.lat, h.lng) / 1000.0, 1) AS distance_km,
           ROUND(st_bearing (h.lat, h.lng, i.lat, i.lng), 2)         AS dispatch_bearing,
           ROW_NUMBER() OVER (PARTITION BY i.incident_id
                              ORDER BY st_distance(i.lat, i.lng, h.lat, h.lng)) AS rn
    FROM {{zone_name}}.emergency.incidents i
    CROSS JOIN {{zone_name}}.emergency.hospitals h
    WHERE i.severity = 'critical' AND h.trauma_level = 1
)
SELECT incident_id, incident_type, hospital_name AS nearest_trauma1,
       distance_km, dispatch_bearing
FROM trauma1_ranked WHERE rn = 1 ORDER BY incident_id;

-- Coverage, incidents within 3 km of each hospital
SELECT h.hospital_name, h.trauma_level, COUNT(*) AS incidents_in_range
FROM {{zone_name}}.emergency.hospitals h
CROSS JOIN {{zone_name}}.emergency.incidents i
WHERE st_distance(h.lat, h.lng, i.lat, i.lng) <= 3000
GROUP BY h.hospital_name, h.trauma_level
HAVING COUNT(*) > 0
ORDER BY incidents_in_range DESC, h.hospital_name;

-- ============================================================================
-- CLEANUP
-- ============================================================================
DROP DELTA TABLE IF EXISTS {{zone_name}}.emergency.response_zones WITH FILES;
DROP DELTA TABLE IF EXISTS {{zone_name}}.emergency.incidents     WITH FILES;
DROP DELTA TABLE IF EXISTS {{zone_name}}.emergency.hospitals     WITH FILES;
DROP SCHEMA IF EXISTS {{zone_name}}.emergency;

Description

## When to Use Reach for this demo when you need to compose several GIS primitives into a single analytical pipeline: nearest-facility lookup, dispatch bearing, and polygon geofencing on the same dataset. It is the counterpart to the maritime demo, which exercises each st_* function in isolation, here the emphasis is on layered queries (window functions + cross joins + containment) that mirror how a real dispatch service queries Delta tables. ## What You Will Learn 1. How to rank nearest facilities per incident using `ROW_NUMBER() OVER (PARTITION BY ... ORDER BY st_distance(...))` on a CROSS JOIN 2. How `st_distance` (haversine), `st_distance_sphere`, and `st_distance_vincenty` differ in accuracy and why Vincenty diverges on longer routes 3. How to compute dispatch direction with `st_bearing`, `st_azimuth` (radians), and `st_final_bearing` 4. How `st_contains(polygon, lat, lng)` and the inverted `st_within(lat, lng, polygon)` implement geofencing against axis-aligned WKT POLYGON strings 5. How to round-trip coordinates via `st_make_point(lng, lat)` → `st_x` / `st_y` and compute geodesic area with `st_area` ## Prerequisites - DeltaForge GUI running with a workspace open - Write access to a DELTA zone and permission to CREATE SCHEMA - Variables: `{{zone_name}}` (defaults to `external`), `{{data_path}}` pointing at writable storage for three Delta tables (`hospitals`, `incidents`, `response_zones`)

Pitfalls

st_distance returns meters (great-circle haversine), not kilometers, divide by 1000.0 before comparing to human-scale distances or your ASSERT VALUE checks will be off by three orders of magnitude.
Argument order diverges between st_contains(polygon, lat, lng) and st_within(lat, lng, polygon); swapping them silently returns false for every row because the WKT parser will treat a lat/lng pair as invalid geometry and most rows fall outside the implied envelope.
st_make_point takes arguments in (longitude, latitude) order to match WKT POINT convention, while every other function in this demo takes (lat, lng). Passing (lat, lng) produces a valid-looking POINT with the coordinates flipped, and st_x / st_y still return them, so the round-trip assertion passes while the geometry is actually on the wrong side of the planet.
Haversine and spherical-law-of-cosines agree to ~0.3% but Vincenty uses the WGS84 ellipsoid and diverges by several kilometers on long routes; do not assert equality across the three algorithms, assert each to its own independently computed reference value (this demo hard-codes haversine_km = 31.4 and leaves Vincenty unasserted for exactly this reason).
WKT polygons in this demo are axis-aligned rectangles in (lng lat) order and MUST be closed (first vertex = last vertex). An unclosed ring is rejected by the WKT parser with a generic 'invalid geometry' error rather than silently being auto-closed.
st_area returns m^2 using a geodesic (not planar) calculation, so the number depends on latitude, two polygons with identical 0.10° × 0.10° spans have different areas at NYC vs the equator. Do not hard-code an area from a flat-earth calculation.
st_bearing / st_final_bearing differ on long great-circle routes because the great-circle heading rotates along the path; short NYC dispatches make them look identical, do not assume they are interchangeable on intercontinental routes.
CROSS JOIN + st_distance is an O(N × M) scan; the 8 × 12 table here is fine but scaling the pattern to production-sized hospital and incident tables requires pre-filtering with a bounding-box WHERE clause or an H3-indexed join before the distance function is called.
DROP DELTA TABLE ... WITH FILES is required to remove the physical Delta Lake directory, a plain DROP TABLE targets DataFusion in-memory tables only and leaves the storage path intact for the next demo run, which then fails with a 'table already exists' error against the on-disk Delta log.

GIS Emergency Response Network: Multi-Step Spatial Analytics

Syntax

Description

Pitfalls

See Also