From c352ece6f1ed40866bc5d017c31e2427480194d6 Mon Sep 17 00:00:00 2001 From: Martin Spott Date: Mon, 22 Nov 2010 07:51:48 -0800 Subject: [PATCH] Proof-of-concept in Shell for cutting a hole into polygon layers entirely inside the database. --- gisscripts/IntersectsRemoval.sh | 124 ++++++++++++++++++++++++++++++++ 1 file changed, 124 insertions(+) create mode 100755 gisscripts/IntersectsRemoval.sh diff --git a/gisscripts/IntersectsRemoval.sh b/gisscripts/IntersectsRemoval.sh new file mode 100755 index 00000000..d2d40e57 --- /dev/null +++ b/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