---
name: gis-postgis-ops
description: Professional geospatial engineering with PostGIS, GeoPandas, Austrian/global projections (SRID), and interactive Leaflet/MapLibre web maps.
---

# Geospatial Engineering: PostGIS, GeoData & Interactive Web Maps

Comprehensive patterns for storing, querying, transforming, and visualizing spatial data using PostGIS, Python Geo-tooling, and modern web mapping libraries (Leaflet & MapLibre).

---

## 1. PostGIS Database Architecture

PostGIS adds spatial objects and functions to PostgreSQL.

### Enabling & Setup
```sql
CREATE EXTENSION IF NOT EXISTS postgis;

-- Verify version
SELECT PostGIS_Full_Version();
```

### Geometry vs. Geography & Coordinate Reference Systems (CRS)

| Type | Best For | Measurements | Distance Unit |
| :--- | :--- | :--- | :--- |
| **`GEOMETRY`** | Local/national projections, CAD, planar calculations, fast indexing | Planar Cartesian | Units of SRID (e.g., meters) |
| **`GEOGRAPHY`** | Global GPS coordinates across large distances without projection distortion | Great-circle spherical | Always in meters |

#### Essential Coordinate Systems (SRIDs)
* **`EPSG:4326` (WGS 84):** Standard GPS coordinates (`lon, lat`). Use for data exchange, GeoJSON, and API payloads.
* **`EPSG:3857` (Web Mercator):** Standard projected coordinate system for Web maps (Leaflet, OpenStreetMap, Mapbox, Google Maps).
* **`EPSG:31256` / `EPSG:31287` (Austria MGI / Lambert):** High-precision national grids for Austria (e.g. Burgenland VGD, cadastral surveys). Units are meters.

### Table Definition & Spatial Indexing Pattern
Always use **GiST** indexing for spatial columns. Without a GiST index, spatial queries default to full table scans.

```sql
CREATE TABLE spatial_features (
    id UUID PRIMARY KEY DEFAULT gen_random_uuid(),
    title TEXT NOT NULL,
    category TEXT NOT NULL,
    -- Projected geometry in EPSG:4326 (WGS 84)
    geom GEOMETRY(Geometry, 4326) NOT NULL,
    properties JSONB DEFAULT '{}'::jsonb,
    created_at TIMESTAMPTZ DEFAULT NOW()
);

-- Crucial: GiST Index for 2D bounding-box queries
CREATE INDEX idx_spatial_features_geom ON spatial_features USING GIST (geom);

-- Optional: GIN index for attribute filtering alongside spatial search
CREATE INDEX idx_spatial_features_props ON spatial_features USING GIN (properties);
```

### Essential PostGIS Queries

#### 1. Features within Radius (in meters)
Using geography casting for accurate metric distance on EPSG:4326:
```sql
SELECT id, title,
       ST_Distance(geom::geography, ST_SetSRID(ST_MakePoint(16.53, 47.85), 4326)::geography) AS distance_meters
FROM spatial_features
WHERE ST_DWithin(geom::geography, ST_SetSRID(ST_MakePoint(16.53, 47.85), 4326)::geography, 5000)
ORDER BY distance_meters ASC;
```

#### 2. Bounding Box & Polygon Intersection
```sql
-- Fast bounding box intersection using && operator
SELECT id, title
FROM spatial_features
WHERE geom && ST_MakeEnvelope(16.4, 47.7, 16.7, 48.0, 4326)
  AND ST_Intersects(geom, ST_MakeEnvelope(16.4, 47.7, 16.7, 48.0, 4326));
```

#### 3. Export directly as GeoJSON FeatureCollection
```sql
SELECT jsonb_build_object(
    'type', 'FeatureCollection',
    'features', COALESCE(jsonb_agg(
        jsonb_build_object(
            'type', 'Feature',
            'id', id,
            'geometry', ST_AsGeoJSON(geom)::jsonb,
            'properties', jsonb_build_object(
                'title', title,
                'category', category
            ) || properties
        )
    ), '[]'::jsonb)
) AS geojson
FROM spatial_features
WHERE category = 'park';
```

---

## 2. Interactive Web Maps: Leaflet & MapLibre

Leaflet is lightweight, mature, and works seamlessly with standard GeoJSON layers and tile servers.

### Leaflet Modern Setup (Vanilla JS / TypeScript)

```html
<!-- Leaflet Assets -->
<link rel="stylesheet" href="https://unpkg.com/leaflet@1.9.4/dist/leaflet.css" />
<script src="https://unpkg.com/leaflet@1.9.4/dist/leaflet.js"></script>

<div id="map" style="width: 100%; height: 600px; border-radius: 8px;"></div>
```

```javascript
// Initialize Map centered on Eastern Austria / Neusiedlersee (e.g. 47.85, 16.65)
const map = L.map('map', {
    zoomControl: true,
    scrollWheelZoom: true
}).setView([47.85, 16.65], 11);

// Layer 1: OpenStreetMap Standard
const osmLayer = L.tileLayer('https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
    maxZoom: 19,
    attribution: '&copy; OpenStreetMap contributors'
});

// Layer 2: Basemap.at (High precision Austrian official map)
const basemapAt = L.tileLayer('https://mapsneu.wien.gv.at/basemap/geolandbasemap/normal/google3857/{z}/{y}/{x}.png', {
    maxZoom: 19,
    attribution: 'Datenquelle: <a href="https://basemap.at">basemap.at</a>'
});

// Layer 3: CartoDB Dark Matter (High-tech Dark Theme)
const darkLayer = L.tileLayer('https://{s}.basemaps.cartocdn.com/dark_all/{z}/{x}/{y}{r}.png', {
    maxZoom: 19,
    attribution: '&copy; CartoDB &copy; OpenStreetMap'
});

// Add default basemap
darkLayer.addTo(map);

// Add Layer Control
L.control.layers({
    "Dark Matter": darkLayer,
    "basemap.at (Österreich)": basemapAt,
    "OpenStreetMap": osmLayer
}).addTo(map);

// Dynamic GeoJSON Loading with Popups & Theming
async function loadGeoJsonLayer(url) {
    const res = await fetch(url);
    const geojsonData = await res.json();

    const geoJsonLayer = L.geoJSON(geojsonData, {
        style: (feature) => ({
            color: feature.properties.color || '#3b82f6',
            weight: 2,
            opacity: 0.8,
            fillColor: feature.properties.color || '#60a5fa',
            fillOpacity: 0.35
        }),
        pointToLayer: (feature, latlng) => {
            return L.circleMarker(latlng, {
                radius: 6,
                fillColor: '#38bdf8',
                color: '#ffffff',
                weight: 1.5,
                fillOpacity: 0.9
            });
        },
        onEachFeature: (feature, layer) => {
            const props = feature.properties || {};
            layer.bindPopup(`
                <div style="font-family: sans-serif; font-size: 13px;">
                    <h4 style="margin: 0 0 6px 0; color: #0f172a;">${props.title || 'Feature'}</h4>
                    <p style="margin: 0; color: #475569;">${props.description || ''}</p>
                </div>
            `);
        }
    }).addTo(map);

    // Zoom map to fit the loaded data
    if (geoJsonLayer.getLayers().length > 0) {
        map.fitBounds(geoJsonLayer.getBounds(), { padding: [30, 30] });
    }
}
```

---

## 3. Python Geospatial Stack: GeoPandas & Shapely

### Reading, Transforming, and Writing Spatial Data

```python
import geopandas as gpd
from shapely.geometry import Point
from sqlalchemy import create_engine

# 1. Read GPX, Shapefile, or GeoPackage
# GeoPandas reads directly with fiona/pyogrio:
gdf = gpd.read_file("tracking_data.gpx", layer="tracks")

# 2. Check and Reproject Coordinate System
if gdf.crs != "EPSG:4326":
    gdf = gdf.to_crs(epsg=4326)

# 3. Spatial Join / Point-in-Polygon Check
points_df = gpd.GeoDataFrame(
    [{"name": "Station A", "geometry": Point(16.53, 47.85)}],
    crs="EPSG:4326"
)

# 4. Stream to PostGIS Database
engine = create_engine("postgresql+psycopg2://user:pass@127.0.0.1:5432/gis_db")
gdf.to_postgis("survey_tracks", engine, if_exists="replace", index=False)
```

---

## 4. SurrealDB (Open Notebook Reference)

Open Notebook on your `ki-server` uses **SurrealDB** for document-graph-vector storage.

### Operational Commands & SurrealQL
```bash
# Connect to SurrealDB inside Docker container
docker exec -it surrealdb /surreal sql --endpoint http://localhost:8000 --username root --password root --namespace open_notebook --database open_notebook

# Export full database dump
docker exec surrealdb /surreal export --endpoint http://localhost:8000 --user root --pass root --ns open_notebook --db open_notebook /mydata/backup.surql

# Key SurrealQL queries:
# INFO FOR DB;
# SELECT * FROM source LIMIT 5;
# SELECT * FROM notebook;
```