Raster Data in PostGIS

What is Raster Data?

Raster data represents spatial information as a regular grid of cells (pixels), where each cell contains a value representing some phenomenon at that location. Unlike vector data which represents discrete features, rasters are especially useful for continuous phenomena that vary across space.

In PostGIS, rasters are stored using the raster data type. PostGIS raster support requires the postgis_raster extension, which is separate from the core PostGIS geometry functionality.

Key Properties of Rasters

Each raster in PostGIS includes essential metadata that defines its spatial characteristics:

Inspecting Raster Metadata

-- Basic raster properties
SELECT
  ST_Width(rast) AS width,
  ST_Height(rast) AS height,
  ST_NumBands(rast) AS bands,
  ST_SRID(rast) AS srid,
  ST_ScaleX(rast) AS pixel_width,
  ST_ScaleY(rast) AS pixel_height
FROM elevation_tiles
LIMIT 5;

-- Comprehensive metadata using ST_Metadata
SELECT
  (ST_Metadata(rast)).upperleftx,
  (ST_Metadata(rast)).upperlefty,
  (ST_Metadata(rast)).width,
  (ST_Metadata(rast)).height,
  (ST_Metadata(rast)).scalex,
  (ST_Metadata(rast)).scaley,
  (ST_Metadata(rast)).srid
FROM ndvi_data;

Creating and Loading Raster Data

Creating Raster Tables

-- Enable PostGIS raster extension
CREATE EXTENSION IF NOT EXISTS postgis_raster;

-- Create a table to store raster data
CREATE TABLE elevation_tiles (
  id SERIAL PRIMARY KEY,
  tile_name TEXT,
  rast RASTER
);

-- Create table for NDVI data
CREATE TABLE ndvi_data (
  id SERIAL PRIMARY KEY,
  acquisition_date DATE,
  satellite TEXT,
  rast RASTER
);

-- Create table for land cover classification
CREATE TABLE land_cover (
  id SERIAL PRIMARY KEY,
  year INTEGER,
  region TEXT,
  rast RASTER
);

Loading Raster Data with raster2pgsql

The raster2pgsql utility is the primary tool for loading raster files into PostGIS:

# Basic raster loading
raster2pgsql -s 4326 -I -C -M input.tif public.elevation_tiles > load_raster.sql
psql -d mydb -f load_raster.sql

# Load with tiling (break large rasters into smaller tiles)
raster2pgsql -s 4326 -I -C -M -t 256x256 large_raster.tif public.elevation_tiles > load_tiled.sql

# Load multiple bands
raster2pgsql -s 4326 -I -C -M -N multiband.tif public.satellite_imagery > load_multiband.sql

Raster Data Types and Applications

Working with Raster Values

Extracting Pixel Values

-- Get pixel value at a specific point
SELECT ST_Value(rast, ST_GeomFromText('POINT(77.59 12.97)', 4326)) AS elevation
FROM elevation_tiles
WHERE ST_Intersects(rast, ST_GeomFromText('POINT(77.59 12.97)', 4326));

-- Get values for multiple points
SELECT 
  p.name,
  ST_Value(r.rast, p.geom) AS ndvi_value
FROM points_of_interest p
JOIN ndvi_data r ON ST_Intersects(r.rast, p.geom);

-- Extract values along a line (profile)
SELECT 
  (ST_DumpValues(ST_Clip(rast, line_geom))).val AS elevation_profile
FROM elevation_tiles, 
     (SELECT ST_GeomFromText('LINESTRING(77.5 12.9, 77.6 13.0)', 4326) AS line_geom) AS line
WHERE ST_Intersects(rast, line_geom);

Raster Statistics

-- Calculate basic statistics for a raster
SELECT 
  (ST_SummaryStats(rast)).count AS pixel_count,
  (ST_SummaryStats(rast)).sum AS total_value,
  (ST_SummaryStats(rast)).mean AS average_value,
  (ST_SummaryStats(rast)).stddev AS standard_deviation,
  (ST_SummaryStats(rast)).min AS minimum_value,
  (ST_SummaryStats(rast)).max AS maximum_value
FROM elevation_tiles;

-- Statistics for specific band
SELECT 
  (ST_SummaryStats(rast, 1)).mean AS band1_mean,
  (ST_SummaryStats(rast, 2)).mean AS band2_mean,
  (ST_SummaryStats(rast, 3)).mean AS band3_mean
FROM satellite_imagery;

Spatial Analysis with Rasters

Raster Clipping and Masking

-- Clip raster by polygon boundary
SELECT ST_Clip(rast, boundary_geom) AS clipped_raster
FROM elevation_tiles, administrative_boundaries
WHERE ST_Intersects(rast, boundary_geom);

-- Create mask from polygon
SELECT ST_AsRaster(
  boundary_geom, 
  rast, 
  '8BUI', 
  1, 
  0
) AS mask_raster
FROM administrative_boundaries, elevation_tiles;

Raster Algebra

-- Calculate NDVI from red and near-infrared bands
-- NDVI = (NIR - Red) / (NIR + Red)
SELECT ST_MapAlgebra(
  red_band, 1,
  nir_band, 1,
  '([rast2] - [rast1]) / ([rast2] + [rast1])'::text
) AS ndvi_raster
FROM satellite_data;

-- Reclassify elevation data
SELECT ST_Reclass(
  rast, 1,
  '0-100:1, 100-500:2, 500-1000:3, 1000-2000:4, 2000-9999:5',
  '8BUI'
) AS elevation_classes
FROM elevation_tiles;

Raster Indexing and Performance

Creating Spatial Indexes

-- Create spatial index on raster column
CREATE INDEX idx_elevation_rast_gist ON elevation_tiles USING GIST(ST_ConvexHull(rast));

-- Create index on raster envelope
CREATE INDEX idx_ndvi_envelope ON ndvi_data USING GIST(ST_Envelope(rast));

-- Add raster constraints for better performance
SELECT AddRasterConstraints('public', 'elevation_tiles', 'rast');

Performance Optimization Tips

Raster Output and Visualization

Exporting Rasters

-- Export raster as GeoTIFF
COPY (
  SELECT ST_AsGDALRaster(rast, 'GTiff') 
  FROM elevation_tiles 
  WHERE id = 1
) TO '/path/to/output.tif';

-- Export as PNG for visualization
SELECT ST_AsPNG(rast) FROM ndvi_data WHERE id = 1;

-- Export with color map
SELECT ST_AsPNG(
  ST_ColorMap(rast, 1, 'elevation')
) FROM elevation_tiles;

Real-World Applications

Environmental Monitoring

-- Monitor vegetation changes over time
SELECT 
  year,
  AVG((ST_SummaryStats(rast)).mean) AS avg_ndvi
FROM ndvi_time_series
GROUP BY year
ORDER BY year;

-- Identify areas of vegetation loss
SELECT 
  ST_Union(ST_Clip(rast, study_area)) AS degraded_areas
FROM ndvi_data
WHERE (ST_SummaryStats(rast)).mean < 0.3;

Climate Analysis

-- Calculate temperature trends
SELECT 
  extract(month from date) AS month,
  AVG((ST_SummaryStats(rast)).mean) AS avg_temperature
FROM temperature_data
WHERE date >= '2020-01-01'
GROUP BY extract(month from date)
ORDER BY month;

Best Practices

1. Choose appropriate tile sizes: 256x256 or 512x512 pixels work well for most applications
2. Use correct data types: Match pixel data type to your data range (8BUI for 0-255, 32BF for floating point)
3. Create spatial indexes: Always index raster tables for better query performance
4. Add constraints: Use AddRasterConstraints() to improve query optimization
5. Consider storage: Large rasters can consume significant disk space—use compression when appropriate
6. Validate data: Check for NoData values and ensure proper SRID assignment