Geospatial Data Pipeline

Geospatial Data Pipeline

Expert in processing, optimizing, and visualizing geospatial data at scale.

When to Use

✅ Use for:

• Drone imagery processing and annotation
• GPS track analysis and visualization
• Location-based search (find nearby X)
• Map tile generation for web/mobile
• Coordinate system transformations
• Geofencing and spatial queries
• GeoJSON optimization for web

❌ NOT for:

• Simple address validation (use address APIs)
• Basic distance calculations (use Haversine formula)
• Static map embeds (use Mapbox Static API)
• Geocoding (use Nominatim or Google Geocoding API)

---

Technology Selection

Database: PostGIS vs MongoDB Geospatial

Timeline:

• 2005: PostGIS 1.0 released
• 2012: MongoDB adds geospatial indexes
• 2020: PostGIS 3.0 with improved performance
• 2024: PostGIS remains gold standard for GIS workloads

---

Common Anti-Patterns

Anti-Pattern 1: Storing Coordinates as Strings

Novice thinking: "I'll just store lat/lon as text, it's simple"

Problem: Can't use spatial indexes, queries are slow, no validation.

Wrong approach:

// ❌ String storage, no spatial features
interface Location {
  id: string;
  name: string;
  latitude: string;   // "37.7749"
  longitude: string;  // "-122.4194"
}

// Linear scan for "nearby" queries
async function findNearby(lat: string, lon: string): Promise<Location[]> {
  const all = await db.locations.findAll();

  return all.filter(loc => {
    const distance = calculateDistance(
      parseFloat(lat),
      parseFloat(lon),
      parseFloat(loc.latitude),
      parseFloat(loc.longitude)
    );
    return distance < 5000; // 5km
  });
}

Why wrong: O(N) linear scan, no spatial index, string parsing overhead.

Correct approach:

// ✅ PostGIS GEOGRAPHY type with spatial index
CREATE TABLE locations (
  id SERIAL PRIMARY KEY,
  name VARCHAR(255),
  location GEOGRAPHY(POINT, 4326)  -- WGS84 coordinates
);

-- Spatial index (GiST)
CREATE INDEX idx_locations_geography ON locations USING GIST(location);

-- TypeScript query
async function findNearby(lat: number, lon: number, radiusMeters: number): Promise<Location[]> {
  const query = `
    SELECT id, name, ST_AsGeoJSON(location) as geojson
    FROM locations
    WHERE ST_DWithin(
      location,
      ST_SetSRID(ST_MakePoint($1, $2), 4326)::geography,
      $3
    )
    ORDER BY location <-> ST_SetSRID(ST_MakePoint($1, $2), 4326)::geography
    LIMIT 100
  `;

  return db.query(query, [lon, lat, radiusMeters]);  // &lt;10ms with index
}

Timeline context:

• 2000s: Stored lat/lon as FLOAT columns, did math in app code
• 2010s: PostGIS adoption, spatial indexes
• 2024: GEOGRAPHY type handles Earth curvature automatically

---

Anti-Pattern 2: Not Using Spatial Indexes

Problem: Proximity queries do full table scans.

Wrong approach:

-- ❌ No index, sequential scan
CREATE TABLE drone_images (
  id SERIAL PRIMARY KEY,
  image_url VARCHAR(255),
  location GEOGRAPHY(POINT, 4326)
);

-- This query scans ALL rows
SELECT * FROM drone_images
WHERE ST_DWithin(
  location,
  ST_SetSRID(ST_MakePoint(-122.4194, 37.7749), 4326)::geography,
  1000  -- 1km
);

EXPLAIN output: Seq Scan on drone_images (cost=0.00..1234.56 rows=1 width=123)

Correct approach:

-- ✅ GiST index for spatial queries
CREATE INDEX idx_drone_images_location ON drone_images USING GIST(location);

-- Same query, now uses index
SELECT * FROM drone_images
WHERE ST_DWithin(
  location,
  ST_SetSRID(ST_MakePoint(-122.4194, 37.7749), 4326)::geography,
  1000
);

EXPLAIN output: Bitmap Index Scan on idx_drone_images_location (cost=4.30..78.30 rows=50 width=123)

Performance impact: 10M points, 5km radius query

• Without index: 3.2 seconds (full scan)
• With GiST index: 12ms (99.6% faster)

---

Anti-Pattern 3: Mixing Coordinate Systems

Novice thinking: "Coordinates are just numbers, I can mix them"

Problem: Incorrect distances, misaligned map features.

Wrong approach:

// ❌ Mixing EPSG:4326 (WGS84) and EPSG:3857 (Web Mercator)
const userLocation = {
  lat: 37.7749,   // WGS84
  lon: -122.4194
};

const droneImage = {
  x: -13634876,  // Web Mercator (EPSG:3857)
  y: 4545684
};

// Comparing apples to oranges!
const distance = Math.sqrt(
  Math.pow(userLocation.lon - droneImage.x, 2) +
  Math.pow(userLocation.lat - droneImage.y, 2)
);

Result: Wildly incorrect distance (millions of "units").

Correct approach:

-- ✅ Transform to common coordinate system
SELECT ST_Distance(
  ST_Transform(
    ST_SetSRID(ST_MakePoint(-122.4194, 37.7749), 4326),  -- WGS84
    3857  -- Transform to Web Mercator
  ),
  ST_SetSRID(ST_MakePoint(-13634876, 4545684), 3857)  -- Already Web Mercator
) AS distance_meters;

Or better: Always store in one system (WGS84), transform on display only.

Timeline:

• 2005: Web Mercator (EPSG:3857) introduced by Google Maps
• 2010: Confusion peaks as apps mix WGS84 data with Web Mercator tiles
• 2024: Best practice: Store WGS84, transform to 3857 only for tile rendering

---

Anti-Pattern 4: Loading Huge GeoJSON Files

Problem: 50MB GeoJSON file crashes browser.

Wrong approach:

// ❌ Load entire file into memory
const geoJson = await fetch('/drone-survey-data.geojson').then(r => r.json());

// 50MB of GeoJSON = browser freeze
map.addSource('drone-data', {
  type: 'geojson',
  data: geoJson  // All 10,000 polygons loaded at once
});

Correct approach 1: Vector tiles (pre-chunked)

// ✅ Serve as vector tiles (MBTiles or PMTiles)
map.addSource('drone-data', {
  type: 'vector',
  tiles: ['https://api.example.com/tiles/{z}/{x}/{y}.pbf'],
  minzoom: 10,
  maxzoom: 18
});

// Browser only loads visible tiles

Correct approach 2: GeoJSON simplification + chunking

# Simplify geometry (reduce points)
npm install -g @mapbox/geojson-precision
geojson-precision -p 5 input.geojson output.geojson

# Split into tiles
npm install -g geojson-vt
# Generate tiles programmatically (see scripts/tile_generator.ts)

Correct approach 3: Server-side filtering

// ✅ Only fetch visible bounds
async function fetchVisibleFeatures(bounds: Bounds): Promise<GeoJSON> {
  const response = await fetch(
    `/api/features?bbox=${bounds.west},${bounds.south},${bounds.east},${bounds.north}`
  );
  return response.json();
}

map.on('moveend', async () => {
  const bounds = map.getBounds();
  const geojson = await fetchVisibleFeatures(bounds);
  map.getSource('dynamic-data').setData(geojson);
});

---

Anti-Pattern 5: Euclidean Distance on Spherical Earth

Novice thinking: "Distance is just Pythagorean theorem"

Problem: Incorrect at scale, worse near poles.

Wrong approach:

// ❌ Flat Earth distance (wrong!)
function distanceKm(lat1: number, lon1: number, lat2: number, lon2: number): number {
  const dx = lon2 - lon1;
  const dy = lat2 - lat1;

  return Math.sqrt(dx * dx + dy * dy) * 111.32;  // 111.32 km/degree (WRONG)
}

// Example: San Francisco to New York
const distance = distanceKm(37.7749, -122.4194, 40.7128, -74.0060);
// Returns: ~55 km (WRONG! Actual: ~4,130 km)

Why wrong: Earth is a sphere, not a flat plane.

Correct approach 1: Haversine formula (great circle distance)

// ✅ Haversine formula (spherical Earth)
function haversineKm(lat1: number, lon1: number, lat2: number, lon2: number): number {
  const R = 6371; // Earth radius in km

  const dLat = toRadians(lat2 - lat1);
  const dLon = toRadians(lon2 - lon1);

  const a =
    Math.sin(dLat / 2) * Math.sin(dLat / 2) +
    Math.cos(toRadians(lat1)) * Math.cos(toRadians(lat2)) *
    Math.sin(dLon / 2) * Math.sin(dLon / 2);

  const c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));

  return R * c;
}

// San Francisco to New York
const distance = haversineKm(37.7749, -122.4194, 40.7128, -74.0060);
// Returns: ~4,130 km ✅

Correct approach 2: PostGIS (handles curvature automatically)

-- ✅ PostGIS ST_Distance with GEOGRAPHY
SELECT ST_Distance(
  ST_SetSRID(ST_MakePoint(-122.4194, 37.7749), 4326)::geography,
  ST_SetSRID(ST_MakePoint(-74.0060, 40.7128), 4326)::geography
) / 1000 AS distance_km;
-- Returns: 4130.137 km ✅

Accuracy comparison:

---

Production Checklist

□ PostGIS extension installed and spatial indexes created
□ All coordinates stored in consistent SRID (recommend: 4326)
□ GeoJSON files optimized (<1MB) or served as vector tiles
□ Coordinate transformations use ST_Transform, not manual math
□ Distance calculations use ST_Distance with GEOGRAPHY type
□ Bounding box queries use ST_MakeEnvelope + ST_Intersects
□ Large geometries chunked (not >100KB per feature)
□ Map tiles pre-generated for common zoom levels
□ CORS configured for tile servers
□ Rate limiting on geocoding/reverse geocoding endpoints

---

When to Use vs Avoid

---

References

• /references/coordinate-systems.md - EPSG codes, transformations, Web Mercator vs WGS84
• /references/postgis-guide.md - PostGIS setup, spatial indexes, common queries
• /references/geojson-optimization.md - Simplification, chunking, vector tiles

Scripts

• scripts/geospatial_processor.ts - Process drone imagery, GPS tracks, GeoJSON validation
• scripts/tile_generator.ts - Generate vector tiles (MBTiles/PMTiles) from GeoJSON

---

This skill guides: Geospatial data | PostGIS | GeoJSON | Map tiles | Coordinate systems | Drone data processing | Spatial queries

Coordinate Systems & Transformations

Complete guide to EPSG codes, coordinate systems, and transformations for geospatial data.

The Problem

Different coordinate systems measure the Earth differently.

• WGS84 (EPSG:4326): Latitude/longitude on spherical Earth
• Web Mercator (EPSG:3857): Projected coordinates for web maps
• UTM zones: Localized projections with minimal distortion

Mixing them causes:

• Incorrect distances
• Misaligned map features
• Wrong spatial queries

---

Common Coordinate Systems

EPSG:4326 (WGS84)

What it is: Geographic coordinates (lat/lon) on an ellipsoid

Units: Degrees

• Longitude: -180° to 180°
• Latitude: -90° to 90°

Used for:

• GPS coordinates
• GeoJSON standard
• Database storage

Characteristics:

• Spherical/ellipsoidal model of Earth
• Not equal-area (distorts at poles)
• Not conformal (distorts shapes)
• Global coverage

Example:

{
  "type": "Point",
  "coordinates": [-122.4194, 37.7749]  // [lon, lat] in degrees
}

Distance calculation: Use great circle distance (Haversine or Vincenty)

---

EPSG:3857 (Web Mercator / Pseudo-Mercator)

What it is: Projected coordinates for web mapping

Units: Meters (but not true meters - see note)

• X (easting): -20,037,508 to 20,037,508
• Y (northing): -20,037,508 to 20,037,508

Used for:

• Google Maps, Mapbox, OpenStreetMap
• Map tiles (z/x/y scheme)
• Web map rendering

Characteristics:

• Cylindrical projection
• Conformal (preserves angles/shapes)
• NOT equal-area (distorts at poles)
• Cuts off at ~85°N/S

Example:

{
  "type": "Point",
  "coordinates": [-13634876, 4545684]  // [x, y] in meters (projected)
}

Important: Distances in Web Mercator are NOT accurate. Always transform to WGS84 for distance calculations.

---

UTM Zones (EPSG:326xx / 327xx)

What it is: Universal Transverse Mercator - localized projections

Coverage: 60 zones, each 6° wide

• Northern Hemisphere: EPSG:32601-32660
• Southern Hemisphere: EPSG:32701-32760

Units: Meters

Used for:

• Engineering projects
• Survey data
• High-accuracy local mapping

Characteristics:

• Conformal (preserves shapes)
• Minimal distortion within zone
• Accurate distance calculations
• Limited coverage (one zone at a time)

Example (San Francisco is in UTM Zone 10N):

-- Transform WGS84 to UTM 10N
SELECT ST_Transform(
  ST_SetSRID(ST_MakePoint(-122.4194, 37.7749), 4326),
  32610  -- UTM 10N
);
-- Result: POINT(551420 4182976)

---

Coordinate Order

CRITICAL: Different systems use different order!

Rule of thumb: GeoJSON and PostGIS use (X, Y) = (lon, lat)

Common mistake:

// ❌ WRONG - swapped coordinates
const point = {
  type: "Point",
  coordinates: [37.7749, -122.4194]  // lat, lon (WRONG!)
};

// ✅ CORRECT - GeoJSON is [lon, lat]
const point = {
  type: "Point",
  coordinates: [-122.4194, 37.7749]  // lon, lat ✅
};

---

Transformations

PostGIS Transformations

Convert between coordinate systems:

-- WGS84 → Web Mercator
SELECT ST_Transform(
  ST_SetSRID(ST_MakePoint(-122.4194, 37.7749), 4326),
  3857
);
-- Result: POINT(-13634876 4545684)

-- Web Mercator → WGS84
SELECT ST_Transform(
  ST_SetSRID(ST_MakePoint(-13634876, 4545684), 3857),
  4326
);
-- Result: POINT(-122.4194 37.7749)

-- WGS84 → UTM 10N
SELECT ST_Transform(
  ST_SetSRID(ST_MakePoint(-122.4194, 37.7749), 4326),
  32610
);
-- Result: POINT(551420 4182976)

Always specify SRID:

-- ❌ WRONG - no SRID, PostGIS doesn't know what it is
SELECT ST_Transform(ST_MakePoint(-122.4194, 37.7749), 3857);

-- ✅ CORRECT - explicit SRID
SELECT ST_Transform(
  ST_SetSRID(ST_MakePoint(-122.4194, 37.7749), 4326),
  3857
);

---

JavaScript Transformations (proj4js)

Install:

npm install proj4

Usage:

import proj4 from 'proj4';

// Define projections (built-in for common ones)
const wgs84 = 'EPSG:4326';
const webMercator = 'EPSG:3857';

// WGS84 → Web Mercator
const [x, y] = proj4(wgs84, webMercator, [-122.4194, 37.7749]);
console.log(x, y);  // -13634876, 4545684

// Web Mercator → WGS84
const [lon, lat] = proj4(webMercator, wgs84, [-13634876, 4545684]);
console.log(lon, lat);  // -122.4194, 37.7749

Custom projections:

// Define UTM 10N
proj4.defs([
  [
    'EPSG:32610',
    '+proj=utm +zone=10 +datum=WGS84 +units=m +no_defs'
  ]
]);

// Transform
const [easting, northing] = proj4('EPSG:4326', 'EPSG:32610', [-122.4194, 37.7749]);

---

Finding the Right UTM Zone

Formula: Zone = floor((longitude + 180) / 6) + 1

function getUTMZone(lon: number, lat: number): number {
  const zone = Math.floor((lon + 180) / 6) + 1;

  // Northern hemisphere: 32600 + zone
  // Southern hemisphere: 32700 + zone
  const epsg = lat >= 0 ? 32600 + zone : 32700 + zone;

  return epsg;
}

// San Francisco (-122.4194, 37.7749)
const zone = getUTMZone(-122.4194, 37.7749);
// Returns: 32610 (UTM 10N)

---

Distance Calculations

Wrong: Euclidean Distance in WGS84

// ❌ WRONG - treats degrees as Cartesian units
function badDistance(lon1, lat1, lon2, lat2) {
  const dx = lon2 - lon1;
  const dy = lat2 - lat1;
  return Math.sqrt(dx * dx + dy * dy) * 111320;  // WRONG!
}

Correct: Haversine Formula

// ✅ CORRECT - great circle distance
function haversineDistance(lon1, lat1, lon2, lat2) {
  const R = 6371000; // Earth radius in meters

  const dLat = toRadians(lat2 - lat1);
  const dLon = toRadians(lon2 - lon1);

  const a =
    Math.sin(dLat / 2) * Math.sin(dLat / 2) +
    Math.cos(toRadians(lat1)) * Math.cos(toRadians(lat2)) *
    Math.sin(dLon / 2) * Math.sin(dLon / 2);

  const c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));

  return R * c;
}

Best: PostGIS with GEOGRAPHY

-- PostGIS handles everything
SELECT ST_Distance(
  ST_SetSRID(ST_MakePoint(-122.4194, 37.7749), 4326)::geography,
  ST_SetSRID(ST_MakePoint(-74.0060, 40.7128), 4326)::geography
) AS distance_meters;

Note: ::geography tells PostGIS to use spheroid model

---

Common Scenarios

Scenario 1: Storing User Locations

Best practice: Store in WGS84 (EPSG:4326)

CREATE TABLE user_locations (
  id SERIAL PRIMARY KEY,
  user_id INTEGER,
  location GEOGRAPHY(POINT, 4326),  -- WGS84
  created_at TIMESTAMP
);

-- Index for spatial queries
CREATE INDEX idx_user_locations ON user_locations USING GIST(location);

-- Insert (from GPS)
INSERT INTO user_locations (user_id, location)
VALUES (123, ST_SetSRID(ST_MakePoint(-122.4194, 37.7749), 4326));

Why WGS84:

• GPS coordinates are in WGS84
• Database can handle distance calculations
• Global coverage

---

Scenario 2: Rendering Map Tiles

Convert to Web Mercator for rendering:

// Client-side: Transform WGS84 → Web Mercator for Mapbox
const features = await fetch('/api/features').then(r => r.json());

// Mapbox expects WGS84, but tiles are in Web Mercator
map.addSource('features', {
  type: 'geojson',
  data: features  // GeoJSON in WGS84
});

// Mapbox handles transformation internally

Server-side: Generate tiles in Web Mercator

-- Tile query (Web Mercator)
SELECT ST_AsMVT(q, 'layer', 4096, 'geom') AS mvt
FROM (
  SELECT
    id,
    name,
    ST_AsMVTGeom(
      ST_Transform(location, 3857),  -- Transform to Web Mercator
      ST_TileEnvelope(:z, :x, :y),
      4096,
      256
    ) AS geom
  FROM locations
  WHERE ST_Intersects(
    location,
    ST_Transform(ST_TileEnvelope(:z, :x, :y), 4326)
  )
) AS q;

---

Scenario 3: Engineering/Survey Data

Use local UTM zone for accuracy:

-- Survey points in UTM 10N (San Francisco area)
CREATE TABLE survey_points (
  id SERIAL PRIMARY KEY,
  point_id VARCHAR(50),
  location GEOMETRY(POINT, 32610),  -- UTM 10N
  elevation FLOAT,
  surveyed_at TIMESTAMP
);

-- Accurate distance calculation (within UTM zone)
SELECT ST_Distance(
  (SELECT location FROM survey_points WHERE point_id = 'A'),
  (SELECT location FROM survey_points WHERE point_id = 'B')
) AS distance_meters;

---

Debugging Coordinate Issues

Check SRID

SELECT ST_SRID(location) FROM locations LIMIT 1;
-- Should return 4326 for WGS84

Visualize Coordinates

// If coordinates look huge, probably Web Mercator
const coords = [-13634876, 4545684];  // Web Mercator
console.log('Magnitude:', Math.abs(coords[0]));  // ~13 million = Web Mercator

// If coordinates are small, probably WGS84
const coords2 = [-122.4194, 37.7749];  // WGS84
console.log('Magnitude:', Math.abs(coords2[0]));  // ~122 = WGS84

Validate Bounds

// WGS84 bounds
const isValidWGS84 = (lon, lat) => {
  return lon >= -180 && lon <= 180 && lat >= -90 && lat <= 90;
};

// Web Mercator bounds
const isValidWebMercator = (x, y) => {
  const limit = 20037508.34;
  return Math.abs(x) <= limit && Math.abs(y) <= limit;
};

---

Quick Reference

---

Resources

• EPSG.io - Search coordinate systems
• Proj4 - Transformation library
• PostGIS Transformations
• Understanding Web Mercator

GeoJSON Optimization

Techniques for optimizing GeoJSON files for web performance: simplification, chunking, precision, and vector tiles.

The Problem

Large GeoJSON files crash browsers:

• 50MB GeoJSON = browser freeze
• 10,000 polygons = slow rendering
• High precision coordinates = wasted bytes

Symptoms:

• Map takes 10+ seconds to load
• Browser tab crashes
• Laggy panning/zooming

---

Optimization Strategies

1. Coordinate Precision

Problem: Default precision is often excessive.

// ❌ Excessive precision (8 decimals = ~1mm accuracy)
{
  "type": "Point",
  "coordinates": [-122.41941234, 37.77492345]
}

// ✅ Appropriate precision (5 decimals = ~1m accuracy)
{
  "type": "Point",
  "coordinates": [-122.4194, 37.7749]
}

Precision vs Accuracy:

Tool: geojson-precision

npm install -g @mapbox/geojson-precision

# Reduce to 5 decimals
geojson-precision -p 5 input.geojson output.geojson

Savings: 8 decimals → 5 decimals = ~40% file size reduction

---

2. Geometry Simplification

Problem: Polygons with too many vertices.

Douglas-Peucker Algorithm: Remove points that don't significantly change shape.

npm install -g @mapbox/simplify-geojson

# Tolerance: higher = more aggressive
simplify-geojson -t 0.001 input.geojson output.geojson

Example:

• Before: 10,000 points
• After (tolerance 0.001): 2,500 points
• Reduction: 75%

Visual comparison:

Before: ⚫⚫⚫⚫⚫⚫⚫⚫⚫⚫ (10,000 points)
After:  ⚫-----⚫-----⚫ (2,500 points)

PostGIS:

-- Simplify geometry
SELECT ST_Simplify(geometry, 0.001) FROM polygons;

-- Simplify preserving topology (prevents self-intersections)
SELECT ST_SimplifyPreserveTopology(geometry, 0.001) FROM polygons;

Tuning tolerance:

• 0.0001: Subtle simplification (~10% reduction)
• 0.001: Moderate simplification (~40% reduction)
• 0.01: Aggressive simplification (~70% reduction)

---

3. Remove Unnecessary Properties

Problem: Feature properties bloat file size.

// ❌ Before: 500 bytes per feature
{
  "type": "Feature",
  "properties": {
    "id": 123,
    "name": "Feature 123",
    "description": "Long description...",
    "created_at": "2024-01-15T10:30:00Z",
    "updated_at": "2024-01-16T14:20:00Z",
    "metadata": { /* ... */ },
    "tags": ["tag1", "tag2", "tag3"]
  },
  "geometry": { /* ... */ }
}

// ✅ After: 100 bytes per feature
{
  "type": "Feature",
  "properties": {
    "id": 123,
    "name": "Feature 123"
  },
  "geometry": { /* ... */ }
}

Savings: 500 bytes → 100 bytes = 80% reduction per feature

Script:

const geojson = require('./input.json');

geojson.features = geojson.features.map(f => ({
  type: 'Feature',
  properties: {
    id: f.properties.id,
    name: f.properties.name
  },
  geometry: f.geometry
}));

fs.writeFileSync('output.json', JSON.stringify(geojson));

---

4. Vector Tiles

Problem: Loading entire dataset at once.

Solution: Split into zoom-dependent tiles.

Benefits:

• Only load visible area
• Progressive loading
• Smaller payloads

Tile Pyramid:

z0: 1 tile (whole world)
z5: 1,024 tiles
z10: 1,048,576 tiles
z14: 268,435,456 tiles

Generate with Tippecanoe:

npm install -g @mapbox/tippecanoe

# Generate tiles
tippecanoe -o output.mbtiles -zg --drop-densest-as-needed input.geojson

# Options:
# -zg: Auto-detect max zoom
# --drop-densest-as-needed: Thin out features at low zoom
# -Z<minzoom> -z<maxzoom>: Set zoom range

Serve tiles:

// Express endpoint
app.get('/tiles/:z/:x/:y.pbf', async (req, res) => {
  const { z, x, y } = req.params;

  // Query PostGIS for tile
  const tile = await db.query(
    'SELECT ST_AsMVT(q, $1) FROM (...) q',
    ['layer']
  );

  res.set('Content-Type', 'application/x-protobuf');
  res.send(tile.rows[0].st_asmvt);
});

Client-side:

// Mapbox GL JS
map.addSource('tiles', {
  type: 'vector',
  tiles: ['https://api.example.com/tiles/{z}/{x}/{y}.pbf'],
  minzoom: 0,
  maxzoom: 14
});

---

5. Feature Chunking

Problem: Single GeoJSON file with 10,000 features.

Solution: Split into multiple files.

Strategy 1: Spatial chunking (by bounds)

// Split by bounding box grid
const chunks = [
  { name: 'nw', bounds: [-180, 0, -90, 90] },
  { name: 'ne', bounds: [-90, 0, 0, 90] },
  { name: 'sw', bounds: [-180, -90, -90, 0] },
  { name: 'se', bounds: [-90, -90, 0, 0] }
];

chunks.forEach(chunk => {
  const features = allFeatures.filter(f =>
    featureIntersectsBounds(f, chunk.bounds)
  );

  fs.writeFileSync(
    `${chunk.name}.geojson`,
    JSON.stringify({ type: 'FeatureCollection', features })
  );
});

Strategy 2: Zoom-based chunking

// High-detail features for high zoom only
const detailedFeatures = features.filter(f => getComplexity(f) > 100);
const simpleFeatures = features.filter(f => getComplexity(f) <= 100);

// Save separately
fs.writeFileSync('simple.geojson', JSON.stringify({ type: 'FeatureCollection', features: simpleFeatures }));
fs.writeFileSync('detailed.geojson', JSON.stringify({ type: 'FeatureCollection', features: detailedFeatures }));

// Load conditionally
if (map.getZoom() > 12) {
  loadDetailed();
} else {
  loadSimple();
}

---

6. Compression

Gzip Compression:

# Compress
gzip -k input.geojson  # Creates input.geojson.gz

# Server: nginx auto-gzip
gzip on;
gzip_types application/json;

Brotli Compression (better than gzip):

# Install
npm install -g brotli-cli

# Compress
brotli input.geojson

# Server: nginx brotli module
brotli on;
brotli_types application/json;

Comparison:

---

7. TopoJSON

Problem: Adjacent polygons share borders = duplicate coordinates.

Solution: TopoJSON encodes topology once.

npm install -g topojson

# Convert GeoJSON → TopoJSON
geo2topo input.geojson > output.topojson

# Simplify topology
toposimplify -s 0.001 output.topojson > simplified.topojson

Savings: 30-80% file size reduction (depends on shared borders)

Use when:

• Choropleth maps (countries, states, counties)
• Adjacent polygons with shared boundaries
• NOT useful for scattered points

Client-side:

import { feature } from 'topojson-client';

// Load TopoJSON
const topology = await fetch('/data.topojson').then(r => r.json());

// Convert to GeoJSON
const geojson = feature(topology, topology.objects.layer);

// Add to map
map.addSource('data', { type: 'geojson', data: geojson });

---

Optimization Workflow

Step 1: Analyze

# Check file size
ls -lh input.geojson

# Count features
jq '.features | length' input.geojson

# Check coordinate precision
jq '.features[0].geometry.coordinates' input.geojson

Step 2: Reduce Precision

geojson-precision -p 5 input.geojson temp1.geojson

Step 3: Simplify Geometry

npx @mapbox/geojson-precision -t 0.001 temp1.geojson temp2.geojson

Step 4: Remove Properties

const data = require('./temp2.geojson');
data.features = data.features.map(f => ({
  type: 'Feature',
  properties: { id: f.properties.id, name: f.properties.name },
  geometry: f.geometry
}));
fs.writeFileSync('temp3.geojson', JSON.stringify(data));

Step 5: Compress

gzip -k temp3.geojson

Step 6: Measure

ls -lh input.geojson input.geojson.gz temp3.geojson temp3.geojson.gz

---

Performance Benchmarks

Test: 50MB GeoJSON with 10,000 polygons

Winner: Vector tiles (100x improvement)

---

When to Use Each Strategy

---

Checklist

□ Reduce coordinate precision to 5 decimals
□ Simplify geometry for features >100 points
□ Remove unused properties from features
□ Enable gzip/brotli compression on server
□ Consider vector tiles for files >5MB
□ Use TopoJSON for adjacent polygons
□ Split large files into zoom-based chunks
□ Test load time on slow connection (3G)

---

Resources

• geojson-precision
• Tippecanoe
• TopoJSON
• Simplify.js

PostGIS Guide

Complete guide to setting up PostGIS, creating spatial indexes, and writing efficient spatial queries.

Installation

macOS (Homebrew)

# Install PostgreSQL with PostGIS
brew install postgresql postgis

# Start PostgreSQL
brew services start postgresql

# Create database
createdb myapp_dev

# Enable PostGIS extension
psql myapp_dev -c "CREATE EXTENSION postgis;"
psql myapp_dev -c "CREATE EXTENSION postgis_topology;"  # Optional

# Verify
psql myapp_dev -c "SELECT PostGIS_version();"

Ubuntu/Debian

# Install
sudo apt-get update
sudo apt-get install postgresql postgresql-contrib postgis

# Enable extension
sudo -u postgres psql -d myapp_dev -c "CREATE EXTENSION postgis;"

Docker

# Official PostGIS image
docker run --name postgis -e POSTGRES_PASSWORD=password -p 5432:5432 -d postgis/postgis

# Connect
psql -h localhost -U postgres

---

Data Types

GEOMETRY vs GEOGRAPHY

GEOMETRY (projected, flat Earth):

CREATE TABLE locations_geometry (
  id SERIAL PRIMARY KEY,
  name VARCHAR(255),
  location GEOMETRY(POINT, 4326)  -- Stores as degrees, treats as flat
);

-- Distance in "degree units" (NOT meters!)
SELECT ST_Distance(
  ST_MakePoint(-122.4194, 37.7749),
  ST_MakePoint(-122.4294, 37.7849)
);
-- Returns: 0.0141 (degrees, meaningless for distance)

GEOGRAPHY (spherical, curved Earth):

CREATE TABLE locations_geography (
  id SERIAL PRIMARY KEY,
  name VARCHAR(255),
  location GEOGRAPHY(POINT, 4326)  -- Treats Earth as sphere
);

-- Distance in meters (accurate!)
SELECT ST_Distance(
  ST_SetSRID(ST_MakePoint(-122.4194, 37.7749), 4326)::geography,
  ST_SetSRID(ST_MakePoint(-122.4294, 37.7849), 4326)::geography
);
-- Returns: 1568.5 (meters) ✅

Rule of thumb:

• Use GEOGRAPHY for WGS84 data (GPS coordinates)
• Use GEOMETRY for projected data (UTM, Web Mercator)

---

Geometry Types

-- Point
GEOMETRY(POINT, 4326)
-- Example: ST_MakePoint(-122.4194, 37.7749)

-- LineString
GEOMETRY(LINESTRING, 4326)
-- Example: ST_MakeLine(ARRAY[point1, point2, point3])

-- Polygon
GEOMETRY(POLYGON, 4326)
-- Example: ST_MakePolygon(ST_MakeLine(ARRAY[p1, p2, p3, p4, p1]))

-- MultiPoint
GEOMETRY(MULTIPOINT, 4326)

-- MultiLineString (e.g., roads network)
GEOMETRY(MULTILINESTRING, 4326)

-- MultiPolygon (e.g., islands)
GEOMETRY(MULTIPOLYGON, 4326)

-- GeometryCollection (mixed types)
GEOMETRY(GEOMETRYCOLLECTION, 4326)

---

Creating Spatial Data

From Lat/Lon

-- Single point
INSERT INTO locations (name, location)
VALUES ('San Francisco', ST_SetSRID(ST_MakePoint(-122.4194, 37.7749), 4326));

-- From coordinates table
UPDATE locations
SET location = ST_SetSRID(ST_MakePoint(longitude, latitude), 4326)
WHERE location IS NULL;

From GeoJSON

-- Insert GeoJSON feature
INSERT INTO locations (name, location)
VALUES (
  'Golden Gate Park',
  ST_SetSRID(
    ST_GeomFromGeoJSON('{"type":"Point","coordinates":[-122.4862,37.7694]}'),
    4326
  )
);

-- Bulk import from GeoJSON file
\copy locations(name, location)
FROM PROGRAM 'cat features.geojson | jq -r ''.features[] | [.properties.name, .geometry | tostring] | @csv'''
WITH CSV DELIMITER ',';

From WKT (Well-Known Text)

INSERT INTO locations (name, location)
VALUES (
  'Line',
  ST_SetSRID(
    ST_GeomFromText('LINESTRING(-122.4194 37.7749, -122.4294 37.7849)'),
    4326
  )
);

---

Spatial Indexes

CRITICAL: Without indexes, spatial queries are O(N) sequential scans.

GiST Index (Most Common)

-- Create table
CREATE TABLE drone_images (
  id SERIAL PRIMARY KEY,
  image_url VARCHAR(255),
  location GEOGRAPHY(POINT, 4326),
  captured_at TIMESTAMP
);

-- GiST index for spatial queries
CREATE INDEX idx_drone_images_location ON drone_images USING GIST(location);

-- Analyze table (update statistics)
ANALYZE drone_images;

Performance impact:

• 10M points, 5km radius query
• Without index: 3.2 seconds (seq scan)
• With GiST index: 12ms (99.6% faster)

SP-GiST Index (Space-Partitioning)

-- Better for point data with high cardinality
CREATE INDEX idx_locations_spgist ON locations USING SPGIST(location);

When to use:

• GiST: General-purpose, works for all geometry types
• SP-GiST: Optimized for points, faster for kNN queries

---

Common Spatial Queries

1. Find Nearby (Within Radius)

-- Find all locations within 5km of point
SELECT id, name, ST_AsGeoJSON(location) as geojson
FROM locations
WHERE ST_DWithin(
  location,
  ST_SetSRID(ST_MakePoint(-122.4194, 37.7749), 4326)::geography,
  5000  -- 5000 meters = 5km
)
ORDER BY location <-> ST_SetSRID(ST_MakePoint(-122.4194, 37.7749), 4326)::geography
LIMIT 10;

Key functions:

• ST_DWithin: Returns true if distance <= threshold
• <->: Distance operator (for ORDER BY)

---

2. Calculate Distance

-- Distance between two points (meters)
SELECT ST_Distance(
  (SELECT location FROM locations WHERE id = 1)::geography,
  (SELECT location FROM locations WHERE id = 2)::geography
) AS distance_meters;

-- Distance from all locations to a point
SELECT
  id,
  name,
  ST_Distance(
    location,
    ST_SetSRID(ST_MakePoint(-122.4194, 37.7749), 4326)::geography
  ) / 1000 AS distance_km
FROM locations
ORDER BY distance_km
LIMIT 10;

---

3. Point in Polygon

-- Check if point is inside polygon
SELECT ST_Within(
  ST_SetSRID(ST_MakePoint(-122.4194, 37.7749), 4326),
  (SELECT boundary FROM jurisdictions WHERE name = 'San Francisco')
) AS is_inside;

-- Find all points inside a polygon
SELECT p.id, p.name
FROM points p
WHERE ST_Within(
  p.location,
  (SELECT boundary FROM jurisdictions WHERE name = 'California')
);

---

4. Bounding Box Query

-- Find all points in bounding box
SELECT id, name, ST_AsGeoJSON(location) as geojson
FROM locations
WHERE ST_Intersects(
  location,
  ST_MakeEnvelope(
    -122.5194, 37.7049,  -- west, south
    -122.3194, 37.8449,  -- east, north
    4326
  )
);

Optimization: Bounding box queries use spatial index efficiently.

---

5. Nearest Neighbor (kNN)

-- Find 10 nearest locations
SELECT id, name, ST_AsGeoJSON(location) as geojson
FROM locations
ORDER BY location <-> ST_SetSRID(ST_MakePoint(-122.4194, 37.7749), 4326)::geography
LIMIT 10;

Note: <-> operator uses index for fast kNN search.

---

6. Clustering (Group by Proximity)

-- Group locations within 1km of each other
SELECT
  ST_ClusterDBSCAN(location, eps := 1000, minpoints := 2) OVER () AS cluster_id,
  id,
  name,
  ST_AsGeoJSON(location) as geojson
FROM locations;

---

Output Formats

GeoJSON

-- Single feature
SELECT ST_AsGeoJSON(location) FROM locations WHERE id = 1;

-- Feature with properties
SELECT json_build_object(
  'type', 'Feature',
  'geometry', ST_AsGeoJSON(location)::json,
  'properties', json_build_object('name', name, 'id', id)
) FROM locations WHERE id = 1;

-- FeatureCollection
SELECT json_build_object(
  'type', 'FeatureCollection',
  'features', json_agg(
    json_build_object(
      'type', 'Feature',
      'geometry', ST_AsGeoJSON(location)::json,
      'properties', json_build_object('name', name, 'id', id)
    )
  )
) FROM locations;

WKT (Well-Known Text)

SELECT ST_AsText(location) FROM locations WHERE id = 1;
-- Returns: POINT(-122.4194 37.7749)

Coordinates

-- Longitude
SELECT ST_X(location) FROM locations WHERE id = 1;

-- Latitude
SELECT ST_Y(location) FROM locations WHERE id = 1;

---

Performance Optimization

Use Bounding Box Before Expensive Operations

-- ❌ Slow: ST_Distance on all rows
SELECT id, name
FROM locations
WHERE ST_Distance(
  location,
  ST_SetSRID(ST_MakePoint(-122.4194, 37.7749), 4326)::geography
) < 5000;

-- ✅ Fast: Bounding box pre-filter + ST_Distance
SELECT id, name
FROM locations
WHERE ST_DWithin(
  location,
  ST_SetSRID(ST_MakePoint(-122.4194, 37.7749), 4326)::geography,
  5000
)
AND ST_Distance(
  location,
  ST_SetSRID(ST_MakePoint(-122.4194, 37.7749), 4326)::geography
) < 5000;

Why: ST_DWithin uses bounding box from index, then ST_Distance refines.

---

Simplify Geometries

-- Simplify polygon (reduce points)
SELECT ST_Simplify(geometry, 0.001) FROM complex_polygons;

-- Simplify preserving topology
SELECT ST_SimplifyPreserveTopology(geometry, 0.001) FROM complex_polygons;

Use when: Rendering at low zoom levels (don't need full detail).

---

Avoid GEOGRAPHY for Local Queries

-- ❌ Slow: GEOGRAPHY on local data
CREATE TABLE local_points (
  location GEOGRAPHY(POINT, 4326)
);

-- ✅ Fast: GEOMETRY in local UTM projection
CREATE TABLE local_points (
  location GEOMETRY(POINT, 32610)  -- UTM 10N
);

Why: GEOGRAPHY calculations are slower (spheroid math). For local areas, projected coordinates are faster and accurate.

---

Vector Tile Generation

Mapbox Vector Tiles (MVT)

-- Generate tile at z/x/y
SELECT ST_AsMVT(q, 'layer', 4096, 'geom') AS mvt
FROM (
  SELECT
    id,
    name,
    ST_AsMVTGeom(
      ST_Transform(location, 3857),  -- Transform to Web Mercator
      ST_TileEnvelope(:z, :x, :y),
      4096,  -- Tile extent
      256    -- Buffer
    ) AS geom
  FROM locations
  WHERE ST_Intersects(
    location,
    ST_Transform(ST_TileEnvelope(:z, :x, :y), 4326)
  )
) AS q;

Use in Express:

app.get('/tiles/:z/:x/:y.mvt', async (req, res) => {
  const { z, x, y } = req.params;

  const result = await db.query(
    'SELECT ST_AsMVT(q, $1, 4096, $2) AS mvt FROM ...',
    ['layer', 'geom']
  );

  res.set('Content-Type', 'application/x-protobuf');
  res.send(result.rows[0].mvt);
});

---

Common Pitfalls

1. Forgetting to Set SRID

-- ❌ WRONG - no SRID
INSERT INTO locations (location) VALUES (ST_MakePoint(-122.4194, 37.7749));

-- ✅ CORRECT
INSERT INTO locations (location) VALUES (
  ST_SetSRID(ST_MakePoint(-122.4194, 37.7749), 4326)
);

2. Mixing GEOMETRY and GEOGRAPHY

-- ❌ WRONG - can't compare GEOMETRY with GEOGRAPHY
SELECT ST_Distance(geometry_col, geography_col);

-- ✅ CORRECT - cast to same type
SELECT ST_Distance(geometry_col::geography, geography_col);

3. Using ST_Distance Without Index

-- ❌ SLOW - ST_Distance doesn't use index
SELECT * FROM locations
WHERE ST_Distance(location, point) < 5000;

-- ✅ FAST - ST_DWithin uses index
SELECT * FROM locations
WHERE ST_DWithin(location, point, 5000);

---

Troubleshooting

Check if Extension is Enabled

SELECT * FROM pg_extension WHERE extname = 'postgis';

Verify SRID

SELECT ST_SRID(location) FROM locations LIMIT 1;
-- Should return 4326 for WGS84

Explain Query Plan

EXPLAIN ANALYZE
SELECT * FROM locations
WHERE ST_DWithin(location, ST_SetSRID(ST_MakePoint(-122.4194, 37.7749), 4326)::geography, 5000);

Look for: "Index Scan using idx_..." (good) vs "Seq Scan" (bad)

---

Resources

• PostGIS Documentation
• PostGIS Cheat Sheet
• Spatial Joins
• Performance Tuning