Overlay Analysis

Overlay analysis is one of the most fundamental operations in GIS, allowing you to combine multiple spatial layers to understand how they interact, overlap, or differ. PostGIS provides powerful overlay functions that enable sophisticated spatial analysis for urban planning, environmental assessment, and decision-making. This comprehensive guide covers all essential overlay operations with practical applications.

Understanding Overlay Analysis

Overlay analysis involves comparing and combining multiple spatial layers to understand their spatial relationships. Think of it as stacking transparent maps and analyzing where they overlap, intersect, or differ. This technique is essential for:

• Finding areas of overlap between different land uses
• Identifying zones affected by multiple factors
• Calculating areas of intersection or difference
• Combining datasets to create new analytical layers
• Performing suitability analysis and site selection

Setting Up Sample Data

Let's create realistic datasets representing parks, flood zones, and development areas to demonstrate overlay analysis:

-- Create parks table
CREATE TABLE parks (
  id SERIAL PRIMARY KEY,
  name TEXT,
  park_type TEXT,
  area_hectares NUMERIC,
  geom GEOMETRY(POLYGON, 4326)
);

-- Create flood zones table
CREATE TABLE flood_zones (
  id SERIAL PRIMARY KEY,
  risk_level TEXT,
  return_period INTEGER, -- years
  geom GEOMETRY(POLYGON, 4326)
);

-- Create development zones table
CREATE TABLE development_zones (
  id SERIAL PRIMARY KEY,
  zone_type TEXT,
  development_status TEXT,
  geom GEOMETRY(POLYGON, 4326)
);

-- Insert sample parks (Delhi area coordinates)
INSERT INTO parks (name, park_type, area_hectares, geom) VALUES
('Central Park', 'Urban', 25.5, ST_GeomFromText('POLYGON((77.20 28.60, 77.25 28.60, 77.25 28.65, 77.20 28.65, 77.20 28.60))', 4326)),
('Riverside Park', 'Natural', 18.2, ST_GeomFromText('POLYGON((77.24 28.62, 77.30 28.62, 77.30 28.68, 77.24 28.68, 77.24 28.62))', 4326)),
('Community Garden', 'Community', 5.8, ST_GeomFromText('POLYGON((77.22 28.58, 77.26 28.58, 77.26 28.61, 77.22 28.61, 77.22 28.58))', 4326)),
('Sports Complex', 'Recreation', 12.3, ST_GeomFromText('POLYGON((77.26 28.64, 77.31 28.64, 77.31 28.67, 77.26 28.67, 77.26 28.64))', 4326));

-- Insert flood zones
INSERT INTO flood_zones (risk_level, return_period, geom) VALUES
('High', 10, ST_GeomFromText('POLYGON((77.23 28.61, 77.28 28.61, 77.28 28.66, 77.23 28.66, 77.23 28.61))', 4326)),
('Medium', 50, ST_GeomFromText('POLYGON((77.21 28.59, 77.32 28.59, 77.32 28.69, 77.21 28.69, 77.21 28.59))', 4326)),
('Low', 100, ST_GeomFromText('POLYGON((77.19 28.57, 77.33 28.57, 77.33 28.71, 77.19 28.71, 77.19 28.57))', 4326));

-- Insert development zones
INSERT INTO development_zones (zone_type, development_status, geom) VALUES
('Residential', 'Planned', ST_GeomFromText('POLYGON((77.25 28.59, 77.29 28.59, 77.29 28.63, 77.25 28.63, 77.25 28.59))', 4326)),
('Commercial', 'Under Construction', ST_GeomFromText('POLYGON((77.27 28.65, 77.32 28.65, 77.32 28.69, 77.27 28.69, 77.27 28.65))', 4326)),
('Industrial', 'Approved', ST_GeomFromText('POLYGON((77.18 28.62, 77.23 28.62, 77.23 28.66, 77.18 28.66, 77.18 28.62))', 4326));

-- Create spatial indexes
CREATE INDEX idx_parks_geom ON parks USING GIST (geom);
CREATE INDEX idx_flood_zones_geom ON flood_zones USING GIST (geom);
CREATE INDEX idx_development_zones_geom ON development_zones USING GIST (geom);

ST_Intersection - Finding Overlapping Areas

ST_Intersection returns the common area between geometries, essential for understanding where different layers overlap.

Basic Intersection Analysis

-- Find parks affected by flood zones
SELECT 
  p.name AS park_name,
  p.park_type,
  f.risk_level,
  f.return_period,
  ST_Intersection(p.geom, f.geom) AS affected_area,
  ST_Area(ST_Intersection(p.geom, f.geom)) AS affected_area_degrees,
  ST_Area(ST_Intersection(p.geom, f.geom)::geography) / 10000 AS affected_area_hectares,
  ROUND(
    (ST_Area(ST_Intersection(p.geom, f.geom)::geography) / ST_Area(p.geom::geography)) * 100, 
    1
  ) AS percent_affected
FROM parks p
JOIN flood_zones f ON ST_Intersects(p.geom, f.geom)
ORDER BY percent_affected DESC;

-- Development zones in flood risk areas
SELECT 
  d.zone_type,
  d.development_status,
  f.risk_level,
  ST_Area(ST_Intersection(d.geom, f.geom)::geography) / 10000 AS overlap_hectares,
  ROUND(
    (ST_Area(ST_Intersection(d.geom, f.geom)::geography) / ST_Area(d.geom::geography)) * 100, 
    1
  ) AS percent_in_flood_zone
FROM development_zones d
JOIN flood_zones f ON ST_Intersects(d.geom, f.geom)
WHERE f.risk_level IN ('High', 'Medium')
ORDER BY percent_in_flood_zone DESC;

Multi-Layer Intersection Analysis

-- Find areas where parks, development zones, and flood zones all intersect
SELECT 
  p.name AS park_name,
  d.zone_type,
  f.risk_level,
  ST_Area(
    ST_Intersection(
      ST_Intersection(p.geom, d.geom), 
      f.geom
    )::geography
  ) / 10000 AS triple_intersection_hectares
FROM parks p
JOIN development_zones d ON ST_Intersects(p.geom, d.geom)
JOIN flood_zones f ON ST_Intersects(ST_Intersection(p.geom, d.geom), f.geom)
WHERE ST_Area(ST_Intersection(ST_Intersection(p.geom, d.geom), f.geom)) > 0
ORDER BY triple_intersection_hectares DESC;

-- Create comprehensive intersection analysis
WITH intersections AS (
  SELECT 
    'Parks vs Flood Zones' AS analysis_type,
    COUNT(*) AS intersection_count,
    SUM(ST_Area(ST_Intersection(p.geom, f.geom)::geography)) / 10000 AS total_overlap_hectares
  FROM parks p
  JOIN flood_zones f ON ST_Intersects(p.geom, f.geom)
  
  UNION ALL
  
  SELECT 
    'Development vs Flood Zones' AS analysis_type,
    COUNT(*) AS intersection_count,
    SUM(ST_Area(ST_Intersection(d.geom, f.geom)::geography)) / 10000 AS total_overlap_hectares
  FROM development_zones d
  JOIN flood_zones f ON ST_Intersects(d.geom, f.geom)
  
  UNION ALL
  
  SELECT 
    'Parks vs Development' AS analysis_type,
    COUNT(*) AS intersection_count,
    SUM(ST_Area(ST_Intersection(p.geom, d.geom)::geography)) / 10000 AS total_overlap_hectares
  FROM parks p
  JOIN development_zones d ON ST_Intersects(p.geom, d.geom)
)
SELECT 
  analysis_type,
  intersection_count,
  ROUND(total_overlap_hectares, 2) AS overlap_hectares
FROM intersections
ORDER BY overlap_hectares DESC;

ST_Union - Merging Geometries

ST_Union combines multiple geometries into a single geometry, dissolving internal boundaries.

Basic Union Operations

-- Combine all parks into a single green space area
SELECT 
  'All Parks Combined' AS description,
  ST_Union(geom) AS combined_parks,
  ST_Area(ST_Union(geom)::geography) / 10000 AS total_park_area_hectares,
  COUNT(*) AS number_of_parks
FROM parks;

-- Combine flood zones by risk level
SELECT 
  risk_level,
  ST_Union(geom) AS combined_flood_zone,
  ST_Area(ST_Union(geom)::geography) / 10000 AS total_area_hectares,
  COUNT(*) AS zone_count
FROM flood_zones
GROUP BY risk_level
ORDER BY 
  CASE risk_level 
    WHEN 'High' THEN 1 
    WHEN 'Medium' THEN 2 
    WHEN 'Low' THEN 3 
  END;

-- Create comprehensive flood risk area
SELECT 
  'Combined Flood Risk' AS description,
  ST_Union(geom) AS all_flood_zones,
  ST_Area(ST_Union(geom)::geography) / 10000 AS total_flood_area_hectares
FROM flood_zones;

Advanced Union Analysis

-- Union with attribute aggregation
SELECT 
  park_type,
  ST_Union(geom) AS type_boundary,
  COUNT(*) AS park_count,
  SUM(area_hectares) AS total_area_hectares,
  AVG(area_hectares) AS avg_area_hectares,
  array_agg(name ORDER BY area_hectares DESC) AS park_names
FROM parks
GROUP BY park_type
ORDER BY total_area_hectares DESC;

-- Create development priority zones by combining overlapping areas
WITH development_overlaps AS (
  SELECT 
    d1.zone_type || ' + ' || d2.zone_type AS combined_type,
    ST_Intersection(d1.geom, d2.geom) AS overlap_geom
  FROM development_zones d1
  JOIN development_zones d2 ON d1.id < d2.id AND ST_Intersects(d1.geom, d2.geom)
)
SELECT 
  combined_type,
  ST_Union(overlap_geom) AS priority_zones,
  ST_Area(ST_Union(overlap_geom)::geography) / 10000 AS priority_area_hectares
FROM development_overlaps
GROUP BY combined_type;

ST_Difference - Finding Exclusive Areas

ST_Difference returns the area that exists in the first geometry but not in the second.

Basic Difference Operations

-- Find park areas outside flood zones (safe areas)
SELECT 
  p.name AS park_name,
  p.park_type,
  ST_Difference(p.geom, ST_Union(f.geom)) AS safe_area,
  ST_Area(p.geom::geography) / 10000 AS total_park_hectares,
  ST_Area(ST_Difference(p.geom, ST_Union(f.geom))::geography) / 10000 AS safe_area_hectares,
  ROUND(
    (ST_Area(ST_Difference(p.geom, ST_Union(f.geom))::geography) / ST_Area(p.geom::geography)) * 100, 
    1
  ) AS percent_safe
FROM parks p
CROSS JOIN (SELECT ST_Union(geom) AS geom FROM flood_zones) f
WHERE ST_Area(ST_Difference(p.geom, f.geom)) > 0
ORDER BY percent_safe DESC;

-- Development zones outside high-risk flood areas
SELECT 
  d.zone_type,
  d.development_status,
  ST_Difference(d.geom, f.geom) AS developable_area,
  ST_Area(ST_Difference(d.geom, f.geom)::geography) / 10000 AS developable_hectares,
  ROUND(
    (ST_Area(ST_Difference(d.geom, f.geom)::geography) / ST_Area(d.geom::geography)) * 100, 
    1
  ) AS percent_developable
FROM development_zones d
CROSS JOIN (
  SELECT ST_Union(geom) AS geom 
  FROM flood_zones 
  WHERE risk_level = 'High'
) f
ORDER BY developable_hectares DESC;

Complex Difference Analysis

-- Find areas suitable for development (outside parks and high flood risk)
WITH restricted_areas AS (
  SELECT ST_Union(geom) AS geom FROM parks
  UNION ALL
  SELECT ST_Union(geom) AS geom FROM flood_zones WHERE risk_level = 'High'
),
combined_restrictions AS (
  SELECT ST_Union(geom) AS restricted_geom FROM restricted_areas
)
SELECT 
  d.zone_type,
  d.development_status,
  ST_Area(d.geom::geography) / 10000 AS total_zone_hectares,
  ST_Area(ST_Difference(d.geom, cr.restricted_geom)::geography) / 10000 AS suitable_hectares,
  ROUND(
    (ST_Area(ST_Difference(d.geom, cr.restricted_geom)::geography) / ST_Area(d.geom::geography)) * 100, 
    1
  ) AS percent_suitable
FROM development_zones d
CROSS JOIN combined_restrictions cr
ORDER BY suitable_hectares DESC;

ST_SymDifference - Finding Non-Overlapping Areas

ST_SymDifference returns areas that are in either geometry but not in both (exclusive or).

-- Find areas that are either parks or development zones but not both
SELECT 
  'Parks XOR Development' AS analysis_type,
  ST_SymDifference(p_union.geom, d_union.geom) AS exclusive_areas,
  ST_Area(ST_SymDifference(p_union.geom, d_union.geom)::geography) / 10000 AS exclusive_area_hectares
FROM 
  (SELECT ST_Union(geom) AS geom FROM parks) p_union
CROSS JOIN 
  (SELECT ST_Union(geom) AS geom FROM development_zones) d_union;

-- Compare different flood risk zones
SELECT 
  f1.risk_level || ' XOR ' || f2.risk_level AS comparison,
  ST_Area(ST_SymDifference(f1.geom, f2.geom)::geography) / 10000 AS exclusive_area_hectares
FROM flood_zones f1
CROSS JOIN flood_zones f2
WHERE f1.risk_level < f2.risk_level
ORDER BY exclusive_area_hectares DESC;

Complex Overlay Scenarios

Multi-Criteria Site Suitability Analysis

-- Find optimal locations for new facilities
-- Criteria: Outside high flood risk, near parks, away from industrial zones
WITH suitability_analysis AS (
  SELECT 
    -- Start with medium flood risk areas (some risk acceptable)
    f.geom AS base_area,
    f.risk_level,
    -- Areas within 500m of parks (buffer analysis)
    ST_Intersects(f.geom, ST_Buffer(p.geom::geography, 500)::geometry) AS near_parks,
    -- Areas not overlapping with industrial zones
    NOT ST_Intersects(f.geom, d.geom) AS away_from_industrial
  FROM flood_zones f
  CROSS JOIN parks p
  LEFT JOIN development_zones d ON d.zone_type = 'Industrial' AND ST_Intersects(f.geom, d.geom)
  WHERE f.risk_level IN ('Medium', 'Low')
)
SELECT 
  risk_level,
  COUNT(*) AS potential_sites,
  COUNT(CASE WHEN near_parks AND away_from_industrial THEN 1 END) AS suitable_sites,
  ROUND(
    (COUNT(CASE WHEN near_parks AND away_from_industrial THEN 1 END)::NUMERIC / COUNT(*)) * 100, 
    1
  ) AS suitability_percentage
FROM suitability_analysis
GROUP BY risk_level
ORDER BY suitability_percentage DESC;

Environmental Impact Assessment

-- Assess environmental impact of development on parks
WITH impact_analysis AS (
  SELECT 
    p.name AS park_name,
    p.park_type,
    d.zone_type AS development_type,
    d.development_status,
    ST_Intersection(p.geom, ST_Buffer(d.geom::geography, 200)::geometry) AS impact_zone,
    ST_Area(p.geom::geography) / 10000 AS park_area_hectares
  FROM parks p
  JOIN development_zones d ON ST_DWithin(p.geom::geography, d.geom::geography, 200)
)
SELECT 
  park_name,
  park_type,
  development_type,
  development_status,
  park_area_hectares,
  ST_Area(impact_zone::geography) / 10000 AS impacted_area_hectares,
  ROUND(
    (ST_Area(impact_zone::geography) / (park_area_hectares * 10000)) * 100, 
    1
  ) AS percent_impacted,
  CASE 
    WHEN (ST_Area(impact_zone::geography) / (park_area_hectares * 10000)) * 100 > 50 THEN 'High Impact'
    WHEN (ST_Area(impact_zone::geography) / (park_area_hectares * 10000)) * 100 > 25 THEN 'Medium Impact'
    ELSE 'Low Impact'
  END AS impact_level
FROM impact_analysis
ORDER BY percent_impacted DESC;

Visualization and Export

-- Create summary geometries for visualization
CREATE TABLE overlay_results AS
SELECT 
  'Park-Flood Intersection' AS layer_type,
  p.name || ' - ' || f.risk_level AS label,
  ST_Intersection(p.geom, f.geom) AS geom,
  ST_Area(ST_Intersection(p.geom, f.geom)::geography) / 10000 AS area_hectares
FROM parks p
JOIN flood_zones f ON ST_Intersects(p.geom, f.geom)

UNION ALL

SELECT 
  'Safe Park Areas' AS layer_type,
  p.name || ' - Safe Zone' AS label,
  ST_Difference(p.geom, ST_Union(f.geom)) AS geom,
  ST_Area(ST_Difference(p.geom, ST_Union(f.geom))::geography) / 10000 AS area_hectares
FROM parks p
CROSS JOIN (SELECT ST_Union(geom) AS geom FROM flood_zones WHERE risk_level = 'High') f
WHERE ST_Area(ST_Difference(p.geom, f.geom)) > 0;

-- Export results as GeoJSON for web mapping
SELECT 
  json_build_object(
    'type', 'FeatureCollection',
    'features', json_agg(
      json_build_object(
        'type', 'Feature',
        'geometry', ST_AsGeoJSON(geom)::json,
        'properties', json_build_object(
          'layer_type', layer_type,
          'label', label,
          'area_hectares', ROUND(area_hectares, 2)
        )
      )
    )
  ) AS geojson_output
FROM overlay_results;

Performance Optimization

Best Practices

1. Validate geometries: Use ST_IsValid() and ST_MakeValid() before overlay operations
2. Use appropriate coordinate systems: Project to local systems for accurate area calculations
3. Index spatial columns: Create GIST indexes for better performance
4. Handle edge cases: Check for empty results and null geometries
5. Optimize query order: Filter with bounding boxes before expensive operations
6. Document methodology: Record overlay logic and assumptions
7. Test with known data: Verify results with simple, known cases

Common Applications