Working with Geospatial GIS Data in PostgreSQL with PostGIS: Storage, Indexing, and Queries

Database tutorial - IT technology blog
Database tutorial - IT technology blog

Context & Why: Solving the Proximity Problem

A few years ago, I built my first store locator for a retail chain. My initial approach was hopelessly naive. I stored coordinates in two standard DECIMAL(9,6) columns—latitude and longitude—and calculated distances with a handwritten Haversine SQL formula. It worked during development on 20 test locations. Then we deployed to staging with 120,000 storefronts. Response times blew past 1,800 ms per query, CPU usage spiked to 95%, and the database locked up under modest load.

Why did this fail? Standard B-tree indexes arrange data along a single scalar dimension. They handle numeric IDs, alphabetical strings, and timestamps with ease.

Real-world geography, however, demands two-dimensional filtering across an ellipsoidal surface. A B-tree index can filter latitude or longitude independently, but it cannot evaluate both simultaneously without scanning tens of thousands of candidate rows. Planar trigonometry also distorts real distances because it ignores Earth’s curvature.

PostGIS solves this fundamental limitation. As an open-source spatial extension for PostgreSQL, it adds dedicated spatial types: GEOMETRY for flat planar projections and GEOGRAPHY for spherical ellipsoidal calculations. It brings hundreds of optimized analytical functions alongside Generalized Search Tree (GiST) indexing. By delegating spatial calculations to PostGIS, you push heavy math into the database engine while trimming hundreds of lines of brittle math from your application layer.

In practice, geographic data rarely lands on your desk ready to ingest. Client locations frequently arrive as messy CSV dumps or unformatted spreadsheets. When I need to normalize CSV rows into structured JSON payloads for batch ingestion, I use toolcraft.app/en/tools/data/csv-to-json. It processes files directly in the browser without uploading client data to external servers, making it quick and secure for staging your imports.

Installation: Setting Up PostGIS

Because PostGIS runs as an in-database extension, you install it on top of your existing PostgreSQL cluster. On Ubuntu or Debian running PostgreSQL 16, install the official package via apt:

# Update package lists and install PostGIS for PostgreSQL 16
sudo apt-get update
sudo apt-get install -y postgresql-16-postgis-3

For containerized local workflows, spin up the official Docker image with PostGIS pre-configured:

docker run --name postgis-demo \
  -e POSTGRES_PASSWORD=secretpassword \
  -e POSTGRES_DB=gis_db \
  -p 5432:5432 \
  -d postgis/postgis:16-3.4

Next, connect to your database instance with psql and register the extension:

-- Switch to your target database
\c gis_db

-- Activate the PostGIS spatial engine
CREATE EXTENSION IF NOT EXISTS postgis;

-- Confirm successful installation
SELECT PostGIS_Full_Version();

This query outputs the active PostGIS build alongside linked GEOS, PROJ, and LibXML library versions. Seeing this output confirms your database is primed to process spatial data.

Configuration: Schemas, Data Storage, and Spatial Indexing

Configuring a spatial schema requires selecting the right data type and Spatial Reference System Identifier (SRID). For global consumer applications using standard GPS coordinates, SRID 4326 (WGS 84) serves as the industry standard.

PostGIS offers two core spatial primitives:

  • GEOMETRY: Maps coordinates to a flat Euclidean Cartesian plane. It delivers rapid calculations for localized regional maps, state-level projections (such as State Plane), and CAD applications.
  • GEOGRAPHY: Computes true great-circle distances directly over the WGS 84 ellipsoid in meters. It is ideal for GPS points, delivery radius lookups, and nationwide retail locators.

1. Creating the Table Schema

Here is a production-ready schema for coffee shop locations using GEOGRAPHY(Point, 4326):

CREATE TABLE coffee_shops (
    id SERIAL PRIMARY KEY,
    name VARCHAR(150) NOT NULL,
    address TEXT,
    location GEOGRAPHY(Point, 4326) NOT NULL,
    created_at TIMESTAMPTZ DEFAULT NOW()
);

2. Applying Spatial Indexing (GiST)

Never query coordinates in production without a spatial index. PostGIS uses GiST (Generalized Search Tree) to construct hierarchical R-tree structures over spatial bounding boxes:

CREATE INDEX idx_coffee_shops_location 
ON coffee_shops 
USING GIST (location);

This bounding-box hierarchy allows PostgreSQL to prune irrelevant branches of the search tree immediately, avoiding full-table sequential scans.

3. Storing Coordinates

Keep this rule in mind: PostGIS follows Cartesian (Longitude, Latitude) order (X before Y). Most people intuitively say “Latitude, Longitude”. Inverting those values is the single most common mistake in GIS development, frequently sending New York coordinates into Antarctica.

Construct points with ST_SetSRID(ST_MakePoint(longitude, latitude), 4326)::geography:

-- Insert locations (Longitude first: X = Longitude, Y = Latitude)
INSERT INTO coffee_shops (name, address, location)
VALUES
    ('Central Perk', '90 Bedford St, New York', ST_SetSRID(ST_MakePoint(-74.0054, 40.7323), 4326)::geography),
    ('Midtown Roast', '350 5th Ave, New York', ST_SetSRID(ST_MakePoint(-73.9857, 40.7484), 4326)::geography),
    ('Brooklyn Beans', '100 Water St, Brooklyn', ST_SetSRID(ST_MakePoint(-73.9922, 40.7032), 4326)::geography),
    ('Uptown Espresso', '2880 Broadway, New York', ST_SetSRID(ST_MakePoint(-73.9645, 40.8082), 4326)::geography);

Verification & Monitoring: Practical Queries and Performance Tuning

With data indexed, you can run bounding-box and radius queries to evaluate real query performance.

1. Finding Places Within a Specific Radius (ST_DWithin)

Imagine a user stands outside the Empire State Building at coordinates (-73.9851, 40.7488) looking for coffee within a 2-kilometer (2,000 meters) walk. The ST_DWithin function evaluates this spatial boundary with high efficiency:

SELECT 
    name,
    ROUND(ST_Distance(location, ST_SetSRID(ST_MakePoint(-73.9851, 40.7488), 4326)::geography)::numeric, 1) AS distance_meters
FROM coffee_shops
WHERE ST_DWithin(
    location,
    ST_SetSRID(ST_MakePoint(-73.9851, 40.7488), 4326)::geography,
    2000 -- radius in meters
)
ORDER BY distance_meters ASC;

Because the column uses GEOGRAPHY, distances evaluate natively in meters along Earth’s spheroid. No manual degree-to-meter conversion is necessary.

2. Verifying Index Usage with EXPLAIN ANALYZE

Verify whether PostgreSQL navigates your GiST index or defaults to a full table scan:

EXPLAIN ANALYZE
SELECT id, name
FROM coffee_shops
WHERE ST_DWithin(
    location,
    ST_SetSRID(ST_MakePoint(-73.9851, 40.7488), 4326)::geography,
    2000
);

In a healthy system, you will see Index Scan using idx_coffee_shops_location or Bitmap Index Scan. On a dataset of 500,000 rows, a sequential scan often takes 600 ms to 900 ms. With the GiST index active, that same query finishes in 1.5 ms to 3 ms.

If you see Seq Scan on a large table, check whether an implicit type conversion is discarding your index, or ensure you haven’t wrapped the indexed column inside an unindexed function call.

3. Monitoring Index Health and Disk Footprint

Keep an eye on physical storage to avoid unexpected cache pressure as coordinate volumes grow:

SELECT
    relname AS object_name,
    pg_size_pretty(pg_total_relation_size(c.oid)) AS total_size,
    pg_size_pretty(pg_relation_size(c.oid)) AS table_or_index_size
FROM pg_class c
JOIN pg_namespace n ON n.oid = c.relnamespace
WHERE c.relname IN ('coffee_shops', 'idx_coffee_shops_location')
  AND n.nspname = 'public';

Run VACUUM ANALYZE coffee_shops; after large bulk imports. This refreshes PostgreSQL planner statistics so the cost-based optimizer reliably chooses GiST tree lookups over sequential scans.

Share: