Zonal Operations using PostGIS

Recently I have been reviewing training materials from the QGIS Changelog site. One particular lesson I liked was how to do Zonal Operations. I thought about how to replicate the output produced in the lesson using Cloud storage and the PostGIS database.

The basis of this tutorial assumes a user has the following:

1. Cloud Storage. This could be Minio or Amazon S3 buckets or any other storage providers.
2. A PostgreSQL database instance.
3. ogr_fdw extension installed in the PostgreSQL database.

For a local setup, I am using a docker-compose setup.

version: '3.9'

volumes:

  postgis-data:

  minio_data:

  


services:

  


  minio:

    image: quay.io/minio/minio

    environment:

      - MINIO_ROOT_USER=minio_admin

      - MINIO_ROOT_PASSWORD=secure_minio_secret

    entrypoint: /bin/bash

    command: -c 'minio server /data --console-address ":9001"'

    volumes:

      - minio_data:/mapproxy

    ports:

      - "9000:9000"

      - "9001:9001"

  


  db:

    image: kartoza/postgis:16-3.4

    volumes:

      - postgis-data:/var/lib/postgresql

    environment:

      - POSTGRES_DB=gis

      - POSTGRES_USER=docker

      - POSTGRES_PASS=docker

      - ALLOW_IP_RANGE=0.0.0.0/0

      # Add extensions you need to be enabled by default in the DB. Default are the five specified below

      - POSTGRES_MULTIPLE_EXTENSIONS=postgis,hstore,postgis_topology,postgis_raster,pgrouting,ogr_fdw

      - RUN_AS_ROOT=true

    ports:

      - "25432:5432"

    restart: on-failure

    depends_on:

      minio:

        condition: service_started

    healthcheck:

      test: "PGPASSWORD=docker pg_isready -h 127.0.0.1 -U docker -d gis"

1. Use the docker-compose specified above to bring the services up. docker-compose up -d
2. Navigate to the MinIO UI and create a bucket i.e. postgis.
3. Download sample datasets from the training tutorial and upload the datasets into the bucket you created in step 2.
4. Exec into the PostgreSQL container using the command docker-compose exec db bash.
5. Install postgis which will allow you to access the raster tool raster2pgsql

apt get update;apt install -y postgis

1. Load the raster data into the database

export PGPASSWORD="docker"

raster2pgsql -s 3857 -t 256x256 -I -R /vsicurl/http://minio:9000/postgis/kzn_pop_count.tif kzn_pop_count | psql -d gis -p 5432 -U docker -h localhost

1. Use GDAL to see if you can access the spatial data stored in the MinIO bucket.

ogrinfo -ro -al -so /vsicurl/http://minio:9000/postgis/districts.shp

1. Login to the database using psql command line.

psql -d gis -p 5432 -U docker -h localhost

1. Run he SQL commands to create a server and foreign table accessing the remote vector data.

-- Create remote server to interact with vector layer stored in minio bucket

CREATE SERVER remote_shp

	FOREIGN DATA WRAPPER ogr_fdw

	OPTIONS (

  datasource '/vsicurl/http://minio:9000/postgis/districts.shp',

  format 'ESRI Shapefile' );

	

-- create foreign table to map the vector layer into the DB	

CREATE FOREIGN TABLE "districts" (

	fid integer,

	"geom" geometry(MultiPolygon, 4326),

	"name" varchar ,

	"DISTRICT" varchar ,

	"PROVINCE" varchar)

	SERVER remote_shp

	OPTIONS (layer 'districts');

1. To speed up accessing the foreign table we create a materialized view of the vector layer.

--Also project the geom to match the CRS of the raster layer

CREATE materialized view mv_districts as

SELECT fid,name, "DISTRICT" as district, "PROVINCE" as province, st_transform(geom,3857) as geom

from districts;

1. Run the SQL for generating zonal statistics.

SELECT

    -- provides: count | sum | mean | stddev | min | max

    (ST_SummaryStats(ST_Clip(kzn_pop_count.rast, mv_districts.geom, TRUE))).*,

    mv_districts.name AS name,

    mv_districts.geom AS geometry

  


FROM

    kzn_pop_count, mv_districts;

1. If your raster data has multiple bands the SQL above will not work and should be adjusted in the following part uST_Clip(kzn_pop_count.rast,1, mv_districts.geom or alternatively use GDAL to translate the raster, extracting the band you need before uploading the raster layer in step 3, or pass the -b 1 option in raster2pgsql command .

There are multiple ways to tackle geospatial problems and various software provides different methods to use. Storing data on cloud storage providers allows multiple tools to access the same datasets and thereby allowing various analysis to be formulated.