PostGIS Spatial Queries for Geodata Analysis: Practical Examples
PostGIS turns PostgreSQL into a full-featured geographic information system. Instead of loading data into QGIS or ArcGIS for every query, you can analyze geometry directly in SQL. It's faster, more reproducible, and more scalable.
Setup: Configuring PostGIS
Installation
-- In an existing PostgreSQL database
CREATE EXTENSION IF NOT EXISTS postgis;
CREATE EXTENSION IF NOT EXISTS postgis_topology;
-- Check version
SELECT PostGIS_Full_Version();
OpenStreetMap Data Structure in PostGIS
The standard osm2pgsql schema creates 4 main tables:
planet_osm_point -- Point features (POI, trees, lights)
planet_osm_line -- Line features (roads, rivers, fences)
planet_osm_polygon -- Polygon features (buildings, parks, lakes)
planet_osm_roads -- Major roads (for low-zoom rendering)
Geometry is stored in the way column in EPSG:3857 (Web Mercator) projection.
Beginner Level: Basic Queries
Finding Objects in an Area (Bounding Box)
-- All cafes in central Moscow
SELECT name, amenity, ST_AsText(way) AS geom
FROM planet_osm_point
WHERE amenity = 'cafe'
AND way && ST_Transform(
ST_MakeEnvelope(37.58, 55.74, 37.65, 55.76, 4326),
3857
);
The && operator is a fast bounding box intersection check. It uses the spatial index (GIST).
Counting Objects by Type
-- How many schools, hospitals, pharmacies in an area
SELECT amenity, COUNT(*) AS count
FROM planet_osm_point
WHERE amenity IN ('school', 'hospital', 'pharmacy')
AND way && ST_Transform(
ST_MakeEnvelope(37.5, 55.7, 37.7, 55.8, 4326),
3857
)
GROUP BY amenity
ORDER BY count DESC;
Building Area
-- Total building area in a district (square meters)
SELECT
COUNT(*) AS building_count,
ROUND(SUM(ST_Area(way))::numeric, 0) AS total_area_m2,
ROUND(AVG(ST_Area(way))::numeric, 0) AS avg_area_m2
FROM planet_osm_polygon
WHERE building IS NOT NULL
AND way && ST_Transform(
ST_MakeEnvelope(37.60, 55.74, 37.63, 55.76, 4326),
3857
);
Intermediate Level: Spatial Operations
Buffer Zones (ST_Buffer)
-- Buildings within 500 meters of a metro station
WITH metro_stations AS (
SELECT name, way
FROM planet_osm_point
WHERE railway = 'station'
AND station = 'subway'
AND name LIKE '%Arbatskaya%'
)
SELECT b.name, b.building,
ROUND(ST_Distance(b.way, m.way)::numeric, 0) AS distance_m
FROM planet_osm_polygon b, metro_stations m
WHERE b.building IS NOT NULL
AND ST_DWithin(b.way, m.way, 500)
ORDER BY distance_m;
ST_DWithin is more efficient than ST_Buffer + ST_Intersects because it uses the index.
Nearest Objects (KNN)
-- 10 nearest cafes to a given point
SELECT name,
ROUND(ST_Distance(
way,
ST_Transform(ST_SetSRID(ST_MakePoint(37.6176, 55.7558), 4326), 3857)
)::numeric, 0) AS distance_m
FROM planet_osm_point
WHERE amenity = 'cafe'
ORDER BY way <-> ST_Transform(
ST_SetSRID(ST_MakePoint(37.6176, 55.7558), 4326), 3857
)
LIMIT 10;
The <-> operator is KNN (K-Nearest Neighbors) using the GIST index. Runs in O(log n).
Layer Intersection (ST_Intersects)
-- Parks crossed by rivers
SELECT p.name AS park_name, r.name AS river_name,
ROUND(ST_Length(ST_Intersection(r.way, p.way))::numeric, 0) AS length_in_park_m
FROM planet_osm_polygon p
JOIN planet_osm_line r ON ST_Intersects(p.way, r.way)
WHERE p.leisure = 'park'
AND r.waterway = 'river'
AND p.way && ST_Transform(
ST_MakeEnvelope(37.4, 55.6, 37.8, 55.9, 4326), 3857
);
Object Density (Grid Analysis)
-- Restaurant density on a 500x500 m grid
WITH grid AS (
SELECT ST_SquareGrid(500, ST_Transform(
ST_MakeEnvelope(37.5, 55.7, 37.7, 55.8, 4326), 3857
)) AS cell
)
SELECT
ST_AsGeoJSON(ST_Transform((cell).geom, 4326)) AS grid_cell,
COUNT(p.*) AS restaurant_count
FROM grid
LEFT JOIN planet_osm_point p
ON ST_Intersects(p.way, (cell).geom)
AND p.amenity = 'restaurant'
GROUP BY (cell).geom
HAVING COUNT(p.*) > 0
ORDER BY restaurant_count DESC;
Advanced Level: Optimization and Complex Queries
Materialized Views for Repeated Queries
CREATE MATERIALIZED VIEW mv_buildings_with_floors AS
SELECT
osm_id, name, building,
"building:levels"::int AS floors,
ST_Area(way) AS area_m2, way
FROM planet_osm_polygon
WHERE building IS NOT NULL
AND "building:levels" IS NOT NULL
AND "building:levels" ~ '^\d+$'
WITH DATA;
CREATE INDEX idx_mv_buildings_way ON mv_buildings_with_floors USING GIST (way);
-- Refresh weekly
REFRESH MATERIALIZED VIEW CONCURRENTLY mv_buildings_with_floors;
Optimization with ST_SimplifyPreserveTopology
For large areas (city, region), full geometry is excessive. Simplification speeds queries 5-20x:
SELECT highway, name,
ST_SimplifyPreserveTopology(way, 10) AS simplified_way
FROM planet_osm_line
WHERE highway IN ('primary', 'secondary', 'tertiary')
AND way && ST_Transform(
ST_MakeEnvelope(37.3, 55.5, 37.9, 56.0, 4326), 3857
);
Parameter 10 is the tolerance in meters (for EPSG:3857). For 1:50000 scale, a 50m tolerance is usually imperceptible.
Point Clustering
SELECT
ST_ClusterDBSCAN(way, eps := 200, minpoints := 3)
OVER () AS cluster_id,
name, amenity, way
FROM planet_osm_point
WHERE amenity IS NOT NULL
AND way && ST_Transform(
ST_MakeEnvelope(37.58, 55.74, 37.65, 55.76, 4326), 3857
);
Isochrones (Reachability Zones)
Approximate walkability zone calculation (without road graph):
WITH center AS (
SELECT ST_Transform(
ST_SetSRID(ST_MakePoint(37.6176, 55.7558), 4326), 3857
) AS pt
)
SELECT
'walk_10min' AS zone,
ST_AsGeoJSON(ST_Transform(ST_Buffer(pt, 600), 4326)) AS geom,
(SELECT COUNT(*) FROM planet_osm_point p
WHERE p.amenity IS NOT NULL AND ST_DWithin(p.way, c.pt, 600)) AS poi_count
FROM center c;
Spatial Indexes
Creating a GIST Index
CREATE INDEX IF NOT EXISTS idx_point_way
ON planet_osm_point USING GIST (way);
-- Partial index for frequent queries
CREATE INDEX idx_point_amenity_cafe
ON planet_osm_point USING GIST (way)
WHERE amenity = 'cafe';
Checking Index Usage
EXPLAIN (ANALYZE, BUFFERS)
SELECT COUNT(*)
FROM planet_osm_point
WHERE amenity = 'cafe'
AND way && ST_Transform(
ST_MakeEnvelope(37.58, 55.74, 37.65, 55.76, 4326), 3857
);
If you see Seq Scan instead of Index Scan, the index isn't being used. Causes: outdated statistics (ANALYZE planet_osm_point), table too small, or query returns more than 10-20% of rows.
Common Mistakes
1. Forgot ST_Transform. OSM data is in EPSG:3857, but coordinates are usually in EPSG:4326. Comparing without transformation yields empty results.
2. ST_Intersects without index. For large tables, ST_Intersects without a spatial index is extremely slow. Always add && (bbox filter) or ensure a GIST index exists.
3. ST_Buffer instead of ST_DWithin. ST_Buffer(geom, 500) + ST_Intersects creates a temporary geometry. ST_DWithin(a, b, 500) does the same but uses the index and creates no intermediate objects.
4. Area in degrees. ST_Area for EPSG:4326 returns square degrees. For square meters, transform to UTM: ST_Area(ST_Transform(way, 32637)).
Conclusion
PostGIS gives full control over spatial data through SQL. From simple "find all cafes within radius" to complex analytics with clustering, grids, and materialized views — everything is expressed declaratively and executed at C++ speed. The key to performance is proper indexes and ST_DWithin instead of buffers.