Fast PostGIS based raster statistics

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"