Proof-of-concept in Shell for cutting a hole into polygon layers entirely

inside the database.

gisscripts/IntersectsRemoval.sh

@@ -0,0 +1,124 @@
+#!/bin/bash
+#
+PGHOST=geoscope.optiputer.net
+PGUSER=martin
+PGDATABASE=landcover
+PGSQL2SHPCONN="-h ${PGHOST} -u ${PGUSER} -g wkb_geometry -b -r ${PGDATABASE}"
+PSQL="psql -tA -h ${PGHOST} -U ${PGUSER} -d ${PGDATABASE} -c"
+SHAPE=`basename ${1} | cut -f 1 -d \. | uniq`
+CUTLAYER=tempcutout
+DIFFLAYER=tempdiff
+LAYERPREFIX=cs_
+LOGSCRIPT=${HOME}/cs_intersects.log
+WITHINSCRIPT=${HOME}/cs_withins.log
+WITHOUTSCRIPT=${HOME}/cs_withouts.log
+DUMPDIR=${HOME}/shp
+rm -f ${LOGSCRIPT} ${WITHINSCRIPT} ${WITHOUTSCRIPT}
+mkdir -p ${DUMPDIR}
+
+MODE="testing"  # production, testing
+
+InitCutoutTable () {
+    ${PSQL} "DROP TABLE ${CUTLAYER}"
+    ${PSQL} "DELETE FROM geometry_columns WHERE f_table_name LIKE '${CUTLAYER}'"
+    ${PSQL} "CREATE TABLE ${CUTLAYER} (ogc_fid serial NOT NULL, \
+                 wkb_geometry geometry, \
+                 CONSTRAINT enforce_dims_wkb_geometry CHECK (ndims(wkb_geometry) = 2), \
+                 CONSTRAINT enforce_geotype_wkb_geometry CHECK ((geometrytype(wkb_geometry) = 'POLYGON'::text) \
+                     OR  (geometrytype(wkb_geometry) = 'MULTIPOLYGON'::text)), \
+                 CONSTRAINT enforce_srid_wkb_geometry CHECK (srid(wkb_geometry) = 4326) \
+                 )"
+    ${PSQL} "INSERT INTO geometry_columns (f_table_catalog, f_table_schema, f_table_name, f_geometry_column, coord_dimension, srid, type) \
+                 VALUES ('', 'public', '${CUTLAYER}', 'wkb_geometry', 2, 4326, 'POLYGON')"
+    
+    ${PSQL} "ALTER TABLE ${CUTLAYER} ADD PRIMARY KEY (ogc_fid)"
+    ${PSQL} "ALTER TABLE ${CUTLAYER} ALTER COLUMN wkb_geometry SET STORAGE MAIN"
+}
+
+InitDiffTable () {
+    ${PSQL} "DROP TABLE ${DIFFLAYER}"
+    ${PSQL} "DELETE FROM geometry_columns WHERE f_table_name LIKE '${DIFFLAYER}'"
+    ${PSQL} "CREATE TABLE ${DIFFLAYER} (ogc_fid serial NOT NULL, wkb_geometry geometry NOT NULL, CONSTRAINT enforce_srid_wkb_geometry CHECK (srid(wkb_geometry) = 4326))"
+    ${PSQL} "INSERT INTO geometry_columns (f_table_catalog, f_table_schema, f_table_name, f_geometry_column, coord_dimension, srid, type) \
+                 VALUES ('', 'public', '${DIFFLAYER}', 'wkb_geometry', 2, 4326, 'POLYGON')"
+}
+
+ImportShape () {
+    rm -f ${SHAPE}\.dbf
+    ogr2ogr -f PostgreSQL PG:"host=${PGHOST} user=${PGUSER} dbname=${PGDATABASE}" -append -t_srs "EPSG:4326" ${SHAPE}\.shp -nln ${CUTLAYER} > /dev/null 2>&1
+    RETURN=${?}
+    if [ ${RETURN} -eq 0 ]; then
+        echo "Layer ${SHAPE} erfolgreich nach ${CUTLAYER} importiert"
+        ${PSQL} "CREATE INDEX ${CUTLAYER}_gindex ON ${CUTLAYER} USING GIST (wkb_geometry GIST_GEOMETRY_OPS)"
+    else
+        echo "Fehler beim Import von ${SHAPE} nach ${CUTLAYER} - exiting !!"
+        exit 1
+    fi
+}
+
+CheckWithin () {
+    GEOM=${1}
+    CSLAYER=`echo ${GEOM} | awk -F\# '{print $1}'`
+    OGC_FID=`echo ${GEOM} | awk -F\# '{print $2}'`
+    CHECKWITHINSTRING="SELECT ST_Within((SELECT wkb_geometry FROM ${CSLAYER} WHERE ogc_fid = ${OGC_FID}), (SELECT wkb_geometry FROM ${CUTLAYER}))"
+    RETURN=`${PSQL} "${CHECKWITHINSTRING}"`
+    if [ ${RETURN} == "t" ]; then
+        DELSTRING="DELETE FROM ${CSLAYER} WHERE ogc_fid = ${OGC_FID}"
+        if [ ${MODE} == "testing" ]; then
+            echo ${DELSTRING} >> ${WITHINSCRIPT}
+        fi
+    else
+        DIFFSTRING="SELECT ST_Difference((SELECT wkb_geometry FROM ${CSLAYER} WHERE ogc_fid = ${OGC_FID}), (SELECT wkb_geometry FROM ${CUTLAYER}))"
+        ${PSQL} "INSERT INTO ${DIFFLAYER} (wkb_geometry) ${DIFFSTRING}"
+        if [ ${MODE} == "testing" ]; then
+            echo ${GEOM} >> ${WITHOUTSCRIPT}
+        fi
+    fi
+}
+
+SingularizeDump () {
+    if [ `${PSQL} "SELECT COUNT(*) FROM ${DIFFLAYER}"` -gt 0 ]; then
+        for OGC_FID in `${PSQL} "SELECT ogc_fid FROM ${DIFFLAYER} WHERE NumGeometries(wkb_geometry) IS NOT NULL"`; do
+            ${PSQL} "INSERT INTO ${DIFFLAYER} (wkb_geometry) (SELECT GeometryN(wkb_geometry, generate_series(1, NumGeometries(wkb_geometry))) AS wkb_geometry FROM ${DIFFLAYER} WHERE ogc_fid = ${OGC_FID})"
+            ${PSQL} "DELETE FROM ${DIFFLAYER} WHERE ogc_fid = ${OGC_FID}"
+        done
+        CPSTRING="INSERT INTO ${CSLAYER} (wkb_geometry) (SELECT wkb_geometry FROM ${DIFFLAYER})"
+        if [ ${MODE} == "testing" ]; then
+            pgsql2shp -f ${DUMPDIR}/${CSLAYER}.shp ${PGSQL2SHPCONN} "SELECT * FROM ${DIFFLAYER}"
+        elif [ ${MODE} == "production" ]; then
+            ${PSQL} "${CPSTRING}"
+        fi
+    fi
+}
+
+CheckIntersects () {
+    SINGULAR=`${PSQL} "SELECT COUNT(*) FROM ${CUTLAYER}"`
+    if [ ${SINGULAR} -ne 1 ]; then
+        echo "Selection requires to contain exactly one single feature - exiting !"
+        exit 1
+    else
+        for CSLAYER in `${PSQL} "SELECT f_table_name FROM geometry_columns WHERE f_table_name LIKE '${LAYERPREFIX}%' AND type LIKE 'POLYGON' ORDER BY f_table_name"`; do
+            InitDiffTable
+            for OGC_FID in `${PSQL} "SELECT ogc_fid FROM ${CSLAYER} WHERE wkb_geometry && (SELECT wkb_geometry FROM ${CUTLAYER}) ORDER BY ogc_fid"`; do
+                RETURN=`${PSQL} "SELECT ST_Intersects((SELECT wkb_geometry FROM ${CUTLAYER}), (SELECT wkb_geometry FROM ${CSLAYER} WHERE ogc_fid = ${OGC_FID}))"`
+                if [ ${RETURN} == "t" ]; then
+                    DELSTRING="DELETE FROM ${CSLAYER} WHERE ogc_fid = ${OGC_FID}"
+                    if [ ${MODE} == "testing" ]; then
+                        echo ${DELSTRING} >> ${LOGSCRIPT}
+                    elif [ ${MODE} == "production" ]; then
+                        ${PSQL} "${DELSTRING}"
+                    fi
+                    CheckWithin "${CSLAYER}#${OGC_FID}"
+                fi
+            done
+            SingularizeDump
+        done
+    fi
+}
+
+
+InitCutoutTable
+ImportShape
+CheckIntersects
+
+# EOF