GIS Maritime Shipping: PostGIS-Compatible Geospatial Functions

Exercises all 18 PostGIS-compatible st_* functions against 5 cargo vessels navigating 10 world ports with 40 GPS position reports, distances, bearings, WKT harbor polygons, and point-in-polygon detection.

Category: spatial

Syntax

-- ============================================================================
-- SETUP, zone, schema, 3 Delta tables
-- ============================================================================
CREATE ZONE IF NOT EXISTS {{zone_name}} TYPE EXTERNAL
    COMMENT 'External and Delta tables, demo datasets';

CREATE SCHEMA IF NOT EXISTS {{zone_name}}.maritime
    COMMENT 'Maritime shipping, ports, vessels, and GPS positions';

CREATE DELTA TABLE IF NOT EXISTS {{zone_name}}.maritime.ports (
    port_id INT, port_name VARCHAR, country VARCHAR,
    lat DOUBLE, lng DOUBLE, harbor_wkt VARCHAR
) LOCATION '{{data_path}}/maritime_ports';

CREATE DELTA TABLE IF NOT EXISTS {{zone_name}}.maritime.vessels (
    vessel_id INT, vessel_name VARCHAR, vessel_type VARCHAR,
    flag_country VARCHAR, deadweight_tons INT
) LOCATION '{{data_path}}/maritime_vessels';

CREATE DELTA TABLE IF NOT EXISTS {{zone_name}}.maritime.positions (
    position_id INT, vessel_id INT, lat DOUBLE, lng DOUBLE,
    speed_knots DOUBLE, recorded_at VARCHAR
) LOCATION '{{data_path}}/maritime_positions';

-- ============================================================================
-- QUERIES, one per st_* family
-- ============================================================================
-- Great-circle distances between port pairs (km, haversine)
SELECT p1.port_name || ' -> ' || p2.port_name AS route,
       ROUND(st_distance(p1.lat, p1.lng, p2.lat, p2.lng) / 1000.0, 0) AS distance_km
FROM {{zone_name}}.maritime.ports p1
CROSS JOIN {{zone_name}}.maritime.ports p2
WHERE (p1.port_id, p2.port_id) IN ((1,3),(2,10),(5,1),(7,9),(8,4),(1,4))
ORDER BY distance_km;

-- Algorithm comparison, haversine vs sphere vs Vincenty (Singapore -> Rotterdam)
SELECT ROUND(st_distance_haversine(p1.lat, p1.lng, p2.lat, p2.lng) / 1000.0, 0) AS haversine_km,
       ROUND(st_distance_sphere  (p1.lat, p1.lng, p2.lat, p2.lng) / 1000.0, 0) AS sphere_km,
       ROUND(st_distance_vincenty(p1.lat, p1.lng, p2.lat, p2.lng) / 1000.0, 0) AS vincenty_km
FROM {{zone_name}}.maritime.ports p1 CROSS JOIN {{zone_name}}.maritime.ports p2
WHERE p1.port_id = 1 AND p2.port_id = 2;

-- Navigation bearings
SELECT p1.port_name || ' -> ' || p2.port_name AS route,
       ROUND(st_bearing      (p1.lat, p1.lng, p2.lat, p2.lng), 2) AS bearing_deg,
       ROUND(st_azimuth      (p1.lat, p1.lng, p2.lat, p2.lng), 4) AS azimuth_rad,
       ROUND(st_final_bearing(p1.lat, p1.lng, p2.lat, p2.lng), 2) AS final_bearing_deg
FROM {{zone_name}}.maritime.ports p1 CROSS JOIN {{zone_name}}.maritime.ports p2
WHERE (p1.port_id, p2.port_id) IN ((1,3),(2,10),(5,1),(8,4))
ORDER BY p1.port_id;

-- Point-in-polygon, which vessel positions are inside a harbor?
SELECT pos.position_id, pos.vessel_id, v.vessel_name, p.port_name, pos.lat, pos.lng
FROM {{zone_name}}.maritime.positions pos
JOIN  {{zone_name}}.maritime.vessels  v ON v.vessel_id = pos.vessel_id
CROSS JOIN {{zone_name}}.maritime.ports p
WHERE st_contains(p.harbor_wkt, pos.lat, pos.lng) = true
ORDER BY pos.position_id;

-- Coordinate round-trip
SELECT port_id, port_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}}.maritime.ports ORDER BY port_id;

-- Harbor areas: Singapore (equator, larger) vs Hamburg (53N, smaller)
SELECT MAX(CASE WHEN port_id = 1 THEN ROUND(st_area(harbor_wkt)/1e6, 2) END) AS area_km2_singapore,
       MAX(CASE WHEN port_id = 6 THEN ROUND(st_area(harbor_wkt)/1e6, 2) END) AS area_km2_hamburg
FROM {{zone_name}}.maritime.ports;

-- Nearest port to each vessel's last at-sea position (st_distance + ROW_NUMBER)
WITH last_sea AS (
    SELECT pos.vessel_id, pos.lat, pos.lng
    FROM {{zone_name}}.maritime.positions pos
    WHERE pos.speed_knots > 0
      AND pos.position_id = (SELECT MAX(p2.position_id)
                             FROM {{zone_name}}.maritime.positions p2
                             WHERE p2.vessel_id = pos.vessel_id AND p2.speed_knots > 0)
),
distances AS (
    SELECT v.vessel_name, ls.lat AS vessel_lat, ls.lng AS vessel_lng,
           p.port_name,
           ROUND(st_distance(ls.lat, ls.lng, p.lat, p.lng) / 1000.0, 0) AS dist_km,
           ROW_NUMBER() OVER (PARTITION BY v.vessel_id
                              ORDER BY st_distance(ls.lat, ls.lng, p.lat, p.lng)) AS rn
    FROM {{zone_name}}.maritime.vessels v
    JOIN last_sea ls ON ls.vessel_id = v.vessel_id
    CROSS JOIN {{zone_name}}.maritime.ports p
)
SELECT vessel_name, vessel_lat, vessel_lng,
       port_name AS nearest_port, dist_km AS nearest_dist_km
FROM distances WHERE rn = 1 ORDER BY vessel_name;

-- ============================================================================
-- CLEANUP
-- ============================================================================
DROP DELTA TABLE IF EXISTS {{zone_name}}.maritime.positions WITH FILES;
DROP DELTA TABLE IF EXISTS {{zone_name}}.maritime.vessels   WITH FILES;
DROP DELTA TABLE IF EXISTS {{zone_name}}.maritime.ports     WITH FILES;
DROP SCHEMA IF EXISTS {{zone_name}}.maritime;

Description

## When to Use Reach for this demo when onboarding to DeltaForge's PostGIS-compatible geometry layer, it is the systematic reference that exercises every st_* function in isolation against independently computed geodesic ground-truth values. Use it to verify that a new storage backend, a new engine build, or a custom function pack returns the expected meters, bearings, and containment results before trusting layered analytics. ## What You Will Learn 1. How to compute intercontinental great-circle distances with `st_distance` (haversine) and compare against `st_distance_sphere` and `st_distance_vincenty` (WGS84 ellipsoid) 2. How to read initial and final compass bearings with `st_bearing`, `st_azimuth` (radians), and `st_final_bearing` 3. How to detect vessels inside a harbor using `st_contains(harbor_wkt, lat, lng)` and its inverse `st_within(lat, lng, harbor_wkt)` 4. How `st_make_point(lng, lat)` round-trips through `st_x` and `st_y` for coordinate construction/extraction 5. How `st_area` on a geodesic polygon shrinks with increasing latitude even when the lat/lng span is fixed (Singapore vs Hamburg at the same 0.02° × 0.02° box) 6. How to combine `st_distance` with a correlated subquery and `ROW_NUMBER()` to find each vessel's nearest port from its last at-sea position ## Prerequisites - DeltaForge GUI running with a workspace open - Write access to an EXTERNAL zone and permission to CREATE SCHEMA - Variables: `{{zone_name}}` (defaults to `external`), `{{data_path}}` for three Delta tables (`ports`, `vessels`, `positions`)

Pitfalls

st_distance_vincenty diverges from st_distance_haversine by ~3 km on a 10,500 km route (Singapore → Rotterdam) because Vincenty models the WGS84 ellipsoid while haversine assumes a sphere; asserting haversine = vincenty breaks on long routes, this demo asserts haversine_km = 10536 and vincenty_km = 10539 independently.
The harbor WKT rings in this demo are ±0.01° boxes (~1–2 km), so a vessel GPS ping at (lat ± 0.005, lng ± 0.005) reads as 'in harbor' even while well outside a real mooring, the lesson is that st_contains is only as precise as the polygon you feed it, not a built-in definition of 'in port'.
WKT POLYGON strings use (lng lat) order, but st_contains and st_within take their scalar arguments as (lat, lng); the two conventions flipping around the same query is the #1 cause of 'every row is false' bugs when copying between formats.
st_make_point(lng, lat) is the PostGIS convention and conflicts with the (lat, lng) argument order used by every other function in this demo, flipping the pair still produces a valid WKT POINT and a successful st_x/st_y round-trip, so the coordinate-swap error silently passes assertions while the point sits in the wrong hemisphere.
Point-in-polygon at a pole or across the ±180° antimeridian breaks with the axis-aligned rectangular polygons used here; none of this demo's ports cross the dateline, but a harbor drawn from (179.9, lat) to (-179.9, lat) is interpreted as a near-global ring and returns true for almost every vessel.
st_area returns square meters with a geodesic (not planar) calculation, so the same ±0.01° box at Singapore (1.26°N) is ~2.47 km² while at Hamburg (53.54°N) it is ~1.47 km²; hard-coded equality assertions across latitudes fail, this demo uses ASSERT EXPRESSION area_singapore > area_hamburg instead.
Vessels at berth (speed_knots = 0) appear 11 times in the positions table but map to only 7 unique ports in the query results, there is no automatic dedup, so 'positions in harbor' counts ports-times-positions if you do not GROUP BY vessel_id or filter by position_id.
The last-at-sea CTE uses a correlated scalar subquery (SELECT MAX(position_id) WHERE vessel_id = pos.vessel_id); DeltaForge supports this pattern only inside WHERE-clause equality, and NOT inside the SELECT projection, see project_datafusion_limitations for the workaround (CTE + ROW_NUMBER + JOIN).
Always use DROP DELTA TABLE ... WITH FILES to delete the physical Delta directory between demo runs, the plain DROP TABLE variant targets DataFusion in-memory tables and leaves the Delta log on disk, which causes the next CREATE to fail with 'table already exists at location'.

GIS Maritime Shipping: PostGIS-Compatible Geospatial Functions

Syntax

Description

Pitfalls

See Also