PostGIS Best Practices¶

Core Principles¶

Spatial Indexing: Use appropriate spatial indexes for performance
Coordinate Systems: Choose correct CRS for your use case
Data Types: Select optimal geometry vs geography types
Query Optimization: Optimize spatial queries for performance
Data Quality: Ensure spatial data integrity and validation

PostGIS Setup and Configuration¶

PostGIS Extension Installation¶

-- Enable PostGIS extension
CREATE EXTENSION IF NOT EXISTS postgis;
CREATE EXTENSION IF NOT EXISTS postgis_topology;
CREATE EXTENSION IF NOT EXISTS postgis_sfcgal;
CREATE EXTENSION IF NOT EXISTS fuzzystrmatch;
CREATE EXTENSION IF NOT EXISTS postgis_tiger_geocoder;

-- Check PostGIS version and capabilities
SELECT PostGIS_Version();
SELECT PostGIS_GEOS_Version();
SELECT PostGIS_Lib_Version();
SELECT PostGIS_Scripts_Build_Date();

-- Check available spatial reference systems
SELECT srid, auth_name, auth_srid, srtext 
FROM spatial_ref_sys 
WHERE auth_name = 'EPSG' 
ORDER BY auth_srid;

Spatial Reference System Configuration¶

-- Create tables with appropriate SRID
CREATE TABLE cities (
    id SERIAL PRIMARY KEY,
    name VARCHAR(100) NOT NULL,
    country VARCHAR(50) NOT NULL,
    geom GEOMETRY(POINT, 4326),              -- WGS84 for global data
    created_at TIMESTAMPTZ DEFAULT CURRENT_TIMESTAMP
);

CREATE TABLE buildings (
    id SERIAL PRIMARY KEY,
    name VARCHAR(200) NOT NULL,
    height_meters INTEGER,
    geom GEOMETRY(POLYGON, 3857),            -- Web Mercator for web mapping
    created_at TIMESTAMPTZ DEFAULT CURRENT_TIMESTAMP
);

CREATE TABLE administrative_boundaries (
    id SERIAL PRIMARY KEY,
    name VARCHAR(200) NOT NULL,
    level VARCHAR(50) NOT NULL,
    geom GEOMETRY(MULTIPOLYGON, 4326),       -- WGS84 for global boundaries
    created_at TIMESTAMPTZ DEFAULT CURRENT_TIMESTAMP
);

-- Create geography table for global calculations
CREATE TABLE global_events (
    id SERIAL PRIMARY KEY,
    event_name VARCHAR(200) NOT NULL,
    location GEOGRAPHY(POINT, 4326),         -- Geography for global distance calculations
    created_at TIMESTAMPTZ DEFAULT CURRENT_TIMESTAMP
);

Spatial Data Types and Storage¶

Geometry vs Geography¶

-- Use GEOMETRY for local/regional data
CREATE TABLE local_parcels (
    id SERIAL PRIMARY KEY,
    parcel_id VARCHAR(50) UNIQUE NOT NULL,
    owner_name VARCHAR(200),
    area_sq_meters NUMERIC(12,2),
    geom GEOMETRY(POLYGON, 4326),            -- Geometry for local operations
    created_at TIMESTAMPTZ DEFAULT CURRENT_TIMESTAMP
);

-- Use GEOGRAPHY for global data and distance calculations
CREATE TABLE global_sensors (
    id SERIAL PRIMARY KEY,
    sensor_id VARCHAR(50) UNIQUE NOT NULL,
    sensor_type VARCHAR(50) NOT NULL,
    location GEOGRAPHY(POINT, 4326),          -- Geography for global operations
    created_at TIMESTAMPTZ DEFAULT CURRENT_TIMESTAMP
);

-- Insert sample data
INSERT INTO cities (name, country, geom) VALUES
    ('New York', 'USA', ST_GeomFromText('POINT(-74.0059 40.7128)', 4326)),
    ('London', 'UK', ST_GeomFromText('POINT(-0.1276 51.5074)', 4326)),
    ('Tokyo', 'Japan', ST_GeomFromText('POINT(139.6917 35.6895)', 4326));

INSERT INTO global_sensors (sensor_id, sensor_type, location) VALUES
    ('SENSOR_001', 'temperature', ST_GeogFromText('POINT(-74.0059 40.7128)')),
    ('SENSOR_002', 'humidity', ST_GeogFromText('POINT(-0.1276 51.5074)')),
    ('SENSOR_003', 'pressure', ST_GeogFromText('POINT(139.6917 35.6895)'));

Spatial Data Validation¶

-- Create function to validate spatial data
CREATE OR REPLACE FUNCTION validate_geometry(geom GEOMETRY)
RETURNS BOOLEAN AS $$
BEGIN
    -- Check if geometry is valid
    IF NOT ST_IsValid(geom) THEN
        RETURN FALSE;
    END IF;

    -- Check if geometry is not empty
    IF ST_IsEmpty(geom) THEN
        RETURN FALSE;
    END IF;

    -- Check if geometry has correct dimension
    IF ST_Dimension(geom) < 0 THEN
        RETURN FALSE;
    END IF;

    RETURN TRUE;
END;
$$ LANGUAGE plpgsql;

-- Add validation constraints
ALTER TABLE cities ADD CONSTRAINT cities_geom_valid 
    CHECK (validate_geometry(geom));

ALTER TABLE buildings ADD CONSTRAINT buildings_geom_valid 
    CHECK (validate_geometry(geom));

-- Create function to fix invalid geometries
CREATE OR REPLACE FUNCTION fix_invalid_geometry(geom GEOMETRY)
RETURNS GEOMETRY AS $$
BEGIN
    IF ST_IsValid(geom) THEN
        RETURN geom;
    ELSE
        RETURN ST_MakeValid(geom);
    END IF;
END;
$$ LANGUAGE plpgsql;

Spatial Indexing Strategies¶

GIST Indexes for Spatial Data¶

-- Create spatial indexes
CREATE INDEX idx_cities_geom_gist ON cities USING GIST (geom);
CREATE INDEX idx_buildings_geom_gist ON buildings USING GIST (geom);
CREATE INDEX idx_administrative_boundaries_geom_gist ON administrative_boundaries USING GIST (geom);
CREATE INDEX idx_global_sensors_location_gist ON global_sensors USING GIST (location);

-- Create partial spatial indexes
CREATE INDEX idx_cities_geom_gist_active ON cities USING GIST (geom) 
WHERE created_at > CURRENT_DATE - INTERVAL '1 year';

-- Create spatial indexes with specific operators
CREATE INDEX idx_buildings_geom_gist_ops ON buildings USING GIST (geom gist_geometry_ops_2d);

Spatial Index Performance¶

-- Analyze spatial index usage
SELECT 
    schemaname,
    tablename,
    indexname,
    idx_scan,
    idx_tup_read,
    idx_tup_fetch
FROM pg_stat_user_indexes
WHERE indexname LIKE '%_geom_gist%' OR indexname LIKE '%_location_gist%'
ORDER BY idx_scan DESC;

-- Check spatial index statistics
SELECT 
    schemaname,
    tablename,
    attname,
    n_distinct,
    correlation
FROM pg_stats
WHERE attname LIKE '%geom%' OR attname LIKE '%location%';

Spatial Query Optimization¶

Spatial Query Patterns¶

-- Point-in-polygon queries
SELECT 
    c.name as city_name,
    c.country,
    ab.name as boundary_name,
    ab.level
FROM cities c
JOIN administrative_boundaries ab ON ST_Contains(ab.geom, c.geom)
WHERE ab.level = 'country';

-- Distance-based queries
SELECT 
    s1.sensor_id as sensor1,
    s2.sensor_id as sensor2,
    ST_Distance(s1.location, s2.location) as distance_meters
FROM global_sensors s1
CROSS JOIN global_sensors s2
WHERE s1.sensor_id < s2.sensor_id
AND ST_DWithin(s1.location, s2.location, 1000)  -- Within 1km
ORDER BY distance_meters;

-- Spatial joins with performance optimization
SELECT 
    b.name as building_name,
    b.height_meters,
    c.name as city_name,
    ST_Distance(b.geom, c.geom) as distance_meters
FROM buildings b
JOIN cities c ON ST_DWithin(b.geom, c.geom, 10000)  -- Within 10km
WHERE b.height_meters > 50
ORDER BY distance_meters;

Spatial Aggregation¶

-- Spatial aggregation functions
SELECT 
    ab.name as boundary_name,
    ab.level,
    COUNT(c.id) as city_count,
    AVG(ST_Area(ab.geom)) as avg_area,
    ST_Union(c.geom) as city_union_geom
FROM administrative_boundaries ab
LEFT JOIN cities c ON ST_Contains(ab.geom, c.geom)
GROUP BY ab.id, ab.name, ab.level, ab.geom;

-- Spatial clustering
SELECT 
    ST_ClusterKMeans(geom, 5) as cluster_id,
    COUNT(*) as point_count,
    ST_Centroid(ST_Collect(geom)) as cluster_center
FROM cities
GROUP BY ST_ClusterKMeans(geom, 5);

Coordinate Reference Systems¶

CRS Transformation¶

-- Transform between coordinate systems
SELECT 
    name,
    geom as original_geom,
    ST_Transform(geom, 3857) as web_mercator_geom,
    ST_Transform(geom, 32633) as utm_geom
FROM cities
WHERE name = 'New York';

-- Create function for CRS transformation
CREATE OR REPLACE FUNCTION transform_to_web_mercator(geom GEOMETRY)
RETURNS GEOMETRY AS $$
BEGIN
    RETURN ST_Transform(geom, 3857);
END;
$$ LANGUAGE plpgsql;

-- Use CRS transformation in queries
SELECT 
    name,
    ST_Area(ST_Transform(geom, 3857)) as area_sq_meters
FROM administrative_boundaries
WHERE level = 'state';

CRS Best Practices¶

-- Create tables with appropriate CRS
CREATE TABLE local_roads (
    id SERIAL PRIMARY KEY,
    road_name VARCHAR(200) NOT NULL,
    road_type VARCHAR(50) NOT NULL,
    geom GEOMETRY(LINESTRING, 4326),         -- WGS84 for global roads
    created_at TIMESTAMPTZ DEFAULT CURRENT_TIMESTAMP
);

CREATE TABLE local_parcels (
    id SERIAL PRIMARY KEY,
    parcel_id VARCHAR(50) UNIQUE NOT NULL,
    owner_name VARCHAR(200),
    geom GEOMETRY(POLYGON, 4326),            -- WGS84 for global parcels
    created_at TIMESTAMPTZ DEFAULT CURRENT_TIMESTAMP
);

-- Create function to get appropriate CRS for region
CREATE OR REPLACE FUNCTION get_regional_crs(geom GEOMETRY)
RETURNS INTEGER AS $$
DECLARE
    centroid POINT;
    longitude NUMERIC;
    latitude NUMERIC;
BEGIN
    centroid := ST_Centroid(geom);
    longitude := ST_X(centroid);
    latitude := ST_Y(centroid);

    -- Return appropriate UTM zone
    IF latitude >= 0 THEN
        RETURN 32600 + (FLOOR((longitude + 180) / 6) + 1);
    ELSE
        RETURN 32700 + (FLOOR((longitude + 180) / 6) + 1);
    END IF;
END;
$$ LANGUAGE plpgsql;

Spatial Analysis Functions¶

Geometric Operations¶

-- Create spatial analysis functions
CREATE OR REPLACE FUNCTION calculate_building_density(
    boundary_geom GEOMETRY,
    building_geom GEOMETRY
)
RETURNS NUMERIC AS $$
DECLARE
    building_count INTEGER;
    boundary_area NUMERIC;
BEGIN
    -- Count buildings within boundary
    SELECT COUNT(*) INTO building_count
    FROM buildings
    WHERE ST_Within(geom, boundary_geom);

    -- Calculate boundary area in square meters
    boundary_area := ST_Area(ST_Transform(boundary_geom, 3857));

    -- Return density (buildings per square kilometer)
    RETURN (building_count / (boundary_area / 1000000));
END;
$$ LANGUAGE plpgsql;

-- Use spatial analysis function
SELECT 
    ab.name as boundary_name,
    calculate_building_density(ab.geom, b.geom) as building_density
FROM administrative_boundaries ab
CROSS JOIN buildings b
WHERE ab.level = 'city';

Spatial Relationships¶

-- Spatial relationship queries
SELECT 
    c.name as city_name,
    b.name as building_name,
    ST_Distance(c.geom, b.geom) as distance_meters,
    CASE 
        WHEN ST_Contains(c.geom, b.geom) THEN 'Inside'
        WHEN ST_Intersects(c.geom, b.geom) THEN 'Intersects'
        WHEN ST_Touches(c.geom, b.geom) THEN 'Touches'
        ELSE 'Separate'
    END as spatial_relationship
FROM cities c
CROSS JOIN buildings b
WHERE ST_DWithin(c.geom, b.geom, 1000)
ORDER BY distance_meters;

Spatial Data Import/Export¶

Data Import Functions¶

-- Create function to import spatial data from GeoJSON
CREATE OR REPLACE FUNCTION import_geojson_features(
    geojson_data JSONB,
    target_table TEXT
)
RETURNS INTEGER AS $$
DECLARE
    feature JSONB;
    feature_count INTEGER := 0;
    geom_wkt TEXT;
    properties JSONB;
BEGIN
    -- Process each feature in the GeoJSON
    FOR feature IN SELECT * FROM jsonb_array_elements(geojson_data->'features') LOOP
        -- Extract geometry as WKT
        geom_wkt := ST_AsText(ST_GeomFromGeoJSON(feature->'geometry'));

        -- Extract properties
        properties := feature->'properties';

        -- Insert into target table
        EXECUTE format('INSERT INTO %I (geom, properties) VALUES (ST_GeomFromText(%L, 4326), %L)',
                      target_table, geom_wkt, properties);

        feature_count := feature_count + 1;
    END LOOP;

    RETURN feature_count;
END;
$$ LANGUAGE plpgsql;

Data Export Functions¶

-- Create function to export spatial data to GeoJSON
CREATE OR REPLACE FUNCTION export_to_geojson(
    table_name TEXT,
    geom_column TEXT DEFAULT 'geom',
    where_clause TEXT DEFAULT '1=1'
)
RETURNS JSONB AS $$
DECLARE
    result JSONB;
    features JSONB;
BEGIN
    -- Build GeoJSON features
    EXECUTE format('
        SELECT jsonb_build_object(
            ''type'', ''FeatureCollection'',
            ''features'', jsonb_agg(
                jsonb_build_object(
                    ''type'', ''Feature'',
                    ''geometry'', ST_AsGeoJSON(%I)::jsonb,
                    ''properties'', row_to_json(t) - ''%I''
                )
            )
        )
        FROM (
            SELECT * FROM %I WHERE %s
        ) t
    ', geom_column, geom_column, table_name, where_clause) INTO result;

    RETURN result;
END;
$$ LANGUAGE plpgsql;

Spatial Data Quality¶

Data Quality Validation¶

-- Create spatial data quality validation
CREATE OR REPLACE FUNCTION validate_spatial_data_quality(
    geom GEOMETRY
)
RETURNS TABLE (
    validation_check TEXT,
    result BOOLEAN,
    message TEXT
) AS $$
BEGIN
    -- Check if geometry is valid
    RETURN QUERY SELECT 'valid_geometry', ST_IsValid(geom), 
        CASE WHEN ST_IsValid(geom) THEN 'Geometry is valid' ELSE 'Geometry is invalid' END;

    -- Check if geometry is not empty
    RETURN QUERY SELECT 'not_empty', NOT ST_IsEmpty(geom),
        CASE WHEN NOT ST_IsEmpty(geom) THEN 'Geometry is not empty' ELSE 'Geometry is empty' END;

    -- Check if geometry has correct dimension
    RETURN QUERY SELECT 'correct_dimension', ST_Dimension(geom) >= 0,
        CASE WHEN ST_Dimension(geom) >= 0 THEN 'Geometry has correct dimension' ELSE 'Geometry has incorrect dimension' END;

    -- Check if geometry is not too complex
    RETURN QUERY SELECT 'not_too_complex', ST_NPoints(geom) < 10000,
        CASE WHEN ST_NPoints(geom) < 10000 THEN 'Geometry is not too complex' ELSE 'Geometry is too complex' END;
END;
$$ LANGUAGE plpgsql;

Spatial Data Cleaning¶

-- Create function to clean spatial data
CREATE OR REPLACE FUNCTION clean_spatial_data(
    geom GEOMETRY
)
RETURNS GEOMETRY AS $$
BEGIN
    -- Remove duplicate points
    geom := ST_RemoveRepeatedPoints(geom);

    -- Simplify geometry if too complex
    IF ST_NPoints(geom) > 1000 THEN
        geom := ST_Simplify(geom, 0.0001);
    END IF;

    -- Fix invalid geometries
    IF NOT ST_IsValid(geom) THEN
        geom := ST_MakeValid(geom);
    END IF;

    -- Remove empty geometries
    IF ST_IsEmpty(geom) THEN
        RETURN NULL;
    END IF;

    RETURN geom;
END;
$$ LANGUAGE plpgsql;

Spatial Performance Monitoring¶

Spatial Query Performance¶

# monitoring/spatial_monitor.py
import psycopg2
import json
from datetime import datetime
import logging

class SpatialMonitor:
    def __init__(self, connection_params):
        self.conn_params = connection_params
        self.setup_logging()

    def setup_logging(self):
        """Setup logging."""
        logging.basicConfig(level=logging.INFO)
        self.logger = logging.getLogger(__name__)

    def get_spatial_index_usage(self):
        """Get spatial index usage statistics."""
        conn = psycopg2.connect(**self.conn_params)

        try:
            with conn.cursor() as cur:
                cur.execute("""
                    SELECT 
                        schemaname,
                        tablename,
                        indexname,
                        idx_scan,
                        idx_tup_read,
                        idx_tup_fetch,
                        pg_size_pretty(pg_relation_size(indexrelid)) as index_size
                    FROM pg_stat_user_indexes
                    WHERE indexname LIKE '%_geom_gist%' OR indexname LIKE '%_location_gist%'
                    ORDER BY idx_scan DESC
                """)

                spatial_indexes = cur.fetchall()
                return spatial_indexes

        except Exception as e:
            self.logger.error(f"Error getting spatial index usage: {e}")
            return []
        finally:
            conn.close()

    def analyze_spatial_query_performance(self):
        """Analyze spatial query performance."""
        conn = psycopg2.connect(**self.conn_params)

        try:
            with conn.cursor() as cur:
                # Test point-in-polygon query performance
                start_time = datetime.now()
                cur.execute("""
                    SELECT COUNT(*) 
                    FROM cities c
                    JOIN administrative_boundaries ab ON ST_Contains(ab.geom, c.geom)
                """)
                pip_time = (datetime.now() - start_time).total_seconds()

                # Test distance query performance
                start_time = datetime.now()
                cur.execute("""
                    SELECT COUNT(*) 
                    FROM global_sensors s1
                    CROSS JOIN global_sensors s2
                    WHERE ST_DWithin(s1.location, s2.location, 1000)
                """)
                distance_time = (datetime.now() - start_time).total_seconds()

                return {
                    'point_in_polygon_time': pip_time,
                    'distance_query_time': distance_time
                }

        except Exception as e:
            self.logger.error(f"Error analyzing spatial query performance: {e}")
            return {}
        finally:
            conn.close()

    def get_spatial_data_statistics(self):
        """Get spatial data statistics."""
        conn = psycopg2.connect(**self.conn_params)

        try:
            with conn.cursor() as cur:
                cur.execute("""
                    SELECT 
                        table_name,
                        column_name,
                        data_type,
                        is_nullable,
                        column_default
                    FROM information_schema.columns
                    WHERE column_name LIKE '%geom%' OR column_name LIKE '%location%'
                    ORDER BY table_name, column_name
                """)

                spatial_columns = cur.fetchall()
                return spatial_columns

        except Exception as e:
            self.logger.error(f"Error getting spatial data statistics: {e}")
            return []
        finally:
            conn.close()

    def generate_spatial_report(self):
        """Generate comprehensive spatial data report."""
        spatial_indexes = self.get_spatial_index_usage()
        query_performance = self.analyze_spatial_query_performance()
        spatial_columns = self.get_spatial_data_statistics()

        report = {
            'report_timestamp': datetime.now().isoformat(),
            'spatial_indexes': spatial_indexes,
            'query_performance': query_performance,
            'spatial_columns': spatial_columns
        }

        return report

# Usage
if __name__ == "__main__":
    monitor = SpatialMonitor({
        'host': 'localhost',
        'database': 'production',
        'user': 'monitor_user',
        'password': 'monitor_password'
    })

    report = monitor.generate_spatial_report()
    print(json.dumps(report, indent=2))

TL;DR Runbook¶

Quick Start¶

-- 1. Enable PostGIS extension
CREATE EXTENSION IF NOT EXISTS postgis;

-- 2. Create spatial tables
CREATE TABLE cities (
    id SERIAL PRIMARY KEY,
    name VARCHAR(100) NOT NULL,
    geom GEOMETRY(POINT, 4326)
);

-- 3. Create spatial indexes
CREATE INDEX idx_cities_geom_gist ON cities USING GIST (geom);

-- 4. Insert spatial data
INSERT INTO cities (name, geom) VALUES
    ('New York', ST_GeomFromText('POINT(-74.0059 40.7128)', 4326));

-- 5. Query spatial data
SELECT name, ST_AsText(geom) as geometry
FROM cities
WHERE ST_DWithin(geom, ST_GeomFromText('POINT(-74.0059 40.7128)', 4326), 1000);

Essential Patterns¶

# Complete PostGIS spatial data setup
def setup_postgis_spatial_data():
    # 1. PostGIS extension setup
    # 2. Spatial data types and storage
    # 3. Spatial indexing strategies
    # 4. Spatial query optimization
    # 5. Coordinate reference systems
    # 6. Spatial analysis functions
    # 7. Data import/export
    # 8. Spatial data quality and monitoring

    print("PostGIS spatial data setup complete!")