Skip to content

Instantly share code, notes, and snippets.

@yellowcap

yellowcap/pgraster-statistics.sh

Created September 5, 2024 10:45
Show Gist options
  • Save yellowcap/cedd25dd19b1641bbcd45dc885ad9495 to your computer and use it in GitHub Desktop.
Save yellowcap/cedd25dd19b1641bbcd45dc885ad9495 to your computer and use it in GitHub Desktop.
Download ZIP
Fast PostGIS based raster statistics
Raw
pgraster-statistics.sh
# Freguesias from
# https://dados.gov.pt/pt/datasets/freguesias-de-portugal/
# Landcover from
# aws s3 cp s3://esa-worldcover/v200/2021/map/ESA_WorldCover_10m_2021_v200_N39W009_Map.tif ~/Desktop/wcpt/
# Create new database
psql -d postgres -c "DROP DATABASE IF EXISTS wc_stats"
psql -d postgres -c "CREATE DATABASE wc_stats"
# Install postgis and postgis-raster into this new database
psql -d wc_stats -c "CREATE EXTENSION postgis"
psql -d wc_stats -c "CREATE EXTENSION postgis_raster"
# Describe the new database (list tables)
psql -d wc_stats -c "\d"
# Load raster into postgis.
/Applications/Postgres.app/Contents/Versions/16/bin/raster2pgsql \
-t 500x500 -I /Users/tam/Desktop/wcpt/ESA_WorldCover_10m_2021_v200_N39W009_Map.tif wc_raster | psql -d wc_stats
# Load freguesias into postgis.
ogr2ogr -nln freguesias -nlt MULTIPOLYGON -t_srs EPSG:4326 -f PostgreSQL "PG:user=tam dbname=wc_stats" /Users/tam/Desktop/cont-aad-caop2017/Cont_AAD_CAOP2017.shp
# Ensure nodata value is set on all rasters.
psql -d wc_stats -c "UPDATE wc_raster SET rast = ST_SetBandNoDataValue(rast, 0)"
# Confirm spatial index (automatically created on import)
psql -d wc_stats -c "SELECT tablename, indexname, indexdef FROM pg_indexes WHERE schemaname = 'public' ORDER BY tablename, indexname"
# Compute zonal statistics.
psql -d wc_stats -c "CREATE TABLE unique_count_stats AS
SELECT Freguesia, (ST_ValueCount(ST_Union(ST_Clip(rast, wkb_geometry)))).*
FROM wc_raster, freguesias
WHERE ST_Intersects(rast, wkb_geometry)
GROUP BY Freguesia"
# Show zonal stats.
psql -d wc_stats -c "SELECT * FROM unique_count_stats LIMIT 25"
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment