From 6da765a09481c6b86169dab585bbea71aabe120e Mon Sep 17 00:00:00 2001 From: Ismael Bejarano Date: Wed, 19 Jun 2024 11:15:23 -0300 Subject: [PATCH] Use mapshaper to simplify a region geometry (#732) Use mapshaper to simplify a region geometry * Document changes in README * Add -f to force geometry updates * Print sizes of files * Add postgis and pgrouting to failing tests * Add unique constraint to region's for (country, name, admin_level) * Use upsert to update region's geometry --- .../068-add_regions_name_index.down.sql | 1 + .../068-add_regions_name_index.up.sql | 1 + resources/planwise/plpgsql/isochrones.sql | 3 + resources/planwise/plpgsql/patches.sql | 3 + scripts/friction/load-friction-raster | 1 + scripts/population/README.md | 14 +++- scripts/population/load-regions | 81 +++++++++++++++---- 7 files changed, 88 insertions(+), 16 deletions(-) create mode 100644 resources/migrations/068-add_regions_name_index.down.sql create mode 100644 resources/migrations/068-add_regions_name_index.up.sql diff --git a/resources/migrations/068-add_regions_name_index.down.sql b/resources/migrations/068-add_regions_name_index.down.sql new file mode 100644 index 000000000..3c38705dc --- /dev/null +++ b/resources/migrations/068-add_regions_name_index.down.sql @@ -0,0 +1 @@ +ALTER TABLE regions DROP CONSTRAINT "regions_name_index"; diff --git a/resources/migrations/068-add_regions_name_index.up.sql b/resources/migrations/068-add_regions_name_index.up.sql new file mode 100644 index 000000000..e4f3df335 --- /dev/null +++ b/resources/migrations/068-add_regions_name_index.up.sql @@ -0,0 +1 @@ +ALTER TABLE regions ADD CONSTRAINT "regions_name_index" UNIQUE ("country", "name", "admin_level"); diff --git a/resources/planwise/plpgsql/isochrones.sql b/resources/planwise/plpgsql/isochrones.sql index 79d530505..0c68616dc 100644 --- a/resources/planwise/plpgsql/isochrones.sql +++ b/resources/planwise/plpgsql/isochrones.sql @@ -1,3 +1,6 @@ +CREATE EXTENSION IF NOT EXISTS postgis; +CREATE EXTENSION IF NOT EXISTS pgrouting; + -- find closest node to a point CREATE OR REPLACE FUNCTION closest_node (original geometry(point, 4326)) returns integer as $$ diff --git a/resources/planwise/plpgsql/patches.sql b/resources/planwise/plpgsql/patches.sql index da46f639f..d77a5d4f2 100644 --- a/resources/planwise/plpgsql/patches.sql +++ b/resources/planwise/plpgsql/patches.sql @@ -1,3 +1,6 @@ +CREATE EXTENSION IF NOT EXISTS postgis; +CREATE EXTENSION IF NOT EXISTS pgrouting; + -- Copied from pgRouting and modified -- Original file: src/alpha_shape/sql/alpha_shape.sql -- Original copyright notice follows diff --git a/scripts/friction/load-friction-raster b/scripts/friction/load-friction-raster index 542955af5..20bc3e50f 100755 --- a/scripts/friction/load-friction-raster +++ b/scripts/friction/load-friction-raster @@ -57,6 +57,7 @@ for region in $regions; do gdalwarp -co COMPRESS=DEFLATE \ -dstnodata -3.4e+38 \ -crop_to_cutline -csql "SELECT the_geom FROM regions WHERE id = ${region}" \ + -cblend 1 \ -cutline PG:"dbname=${POSTGRES_DB} host=${POSTGRES_HOST} port=${POSTGRES_PORT} user=${POSTGRES_USER} password=${POSTGRES_PASSWORD}" \ $raster_file $output_file else diff --git a/scripts/population/README.md b/scripts/population/README.md index 8d5e38db4..19f5c23ae 100644 --- a/scripts/population/README.md +++ b/scripts/population/README.md @@ -15,6 +15,19 @@ $ docker-compose run --rm app bash /app$ scripts/population/load-regions URY Uruguay ``` +### Geometry simplification + +The script simplifies the regions using mapshaper by a distance between points. +We can adjust this parameter with the environment `INTERVAL_DISTANCE`. + +For example to simplify with 120 meters distance + +```sh +$ INTERVAL_DISTANCE=120 ./load-regions -f URY Uruguay +``` + +Here `-f` is used to force updating an existing region. + ## Clip friction layer for new regions The friction raster is clipped to the country-level regions for quicker access @@ -46,4 +59,3 @@ INSERT INTO source_set (name, type, unit, raster_file) VALUES ('Uruguay PPP v2b The important bit of information is the last field. That should be a relative path from the `DATA_PATH` directory to the downloaded raster file. - diff --git a/scripts/population/load-regions b/scripts/population/load-regions index a8f8e0931..8ea19bb17 100755 --- a/scripts/population/load-regions +++ b/scripts/population/load-regions @@ -4,12 +4,35 @@ set -euo pipefail export PGPASSWORD=$POSTGRES_PASSWORD; if [ $# -lt 2 ]; then - echo "Usage $0 " + echo "Usage $0 [-f] " + echo "The option -f forces updating an existing country geometry" exit 1 fi -COUNTRY_CODE=$1 -COUNTRY_NAME=$2 +FORCE= +COUNTRY_CODE= +COUNTRY_NAME= + +while [ $# -gt 0 ]; do + case $1 in + -f|--force) + FORCE=1 + ;; + *) + if [ ! -z "$COUNTRY_NAME" ]; then + echo Unknown arguments specified + exit 1; + elif [ ! -z "$COUNTRY_CODE" ]; then + COUNTRY_NAME=$1 + else + COUNTRY_CODE=$1 + fi + ;; + esac + shift +done + + #LEVELS=${2:-2,4} LEVELS=(0,1) echo " -> Importing $COUNTRY_NAME at levels $LEVELS" @@ -18,14 +41,14 @@ echo " -> Importing $COUNTRY_NAME at levels $LEVELS" REGIONS=$(psql -d $POSTGRES_DB -U $POSTGRES_USER -h $POSTGRES_HOST -p $POSTGRES_PORT -t -A -c "SELECT id, country FROM regions WHERE country = '$COUNTRY_NAME';") REGIONS_LEN=$(echo $REGIONS | wc -w | tr -d '[[:space:]]') -if [ ${REGIONS_LEN} -gt 0 ]; then +if [ ${REGIONS_LEN} -gt 0 -a -z $FORCE ]; then echo " -> Regions for country $COUNTRY_NAME already imported." exit 0 fi DATA_PATH=${DATA_PATH:-/data}/geojson -mkdir -p ${DATA_PATH} +mkdir -p "${DATA_PATH}" # if the folder (still) doesn't exist then exit with error if [[ ! -e ${DATA_PATH}/$COUNTRY_CODE ]]; then @@ -33,25 +56,53 @@ if [[ ! -e ${DATA_PATH}/$COUNTRY_CODE ]]; then exit 1 fi +# Simplification parameter +# From mapshaper wiki https://github.com/mbloch/mapshaper/wiki/Command-Reference#-simplify +# +# interval= Specify simplification amount in units of distance. +# Uses meters when simplifying unprojected datasets in 3D space, +# otherwise uses the same units as the source data. +# In our case since our GeoJSON are not projected (ie. use WGS84 coordinates) +# then this means the resolution for the simplified shape will be in meters. + +INTERVAL_DISTANCE=${INTERVAL_DISTANCE-100} + IFS=','; for i in $LEVELS; do FILE=${DATA_PATH}/${COUNTRY_CODE}/${COUNTRY_CODE}_adm${i}.geojson echo -n " -> Processing $FILE ... " + + OFILE=$FILE + if [[ ! -z $INTERVAL_DISTANCE ]]; then + # Simplify using distance + OFILE=${DATA_PATH}/${COUNTRY_CODE}/${COUNTRY_CODE}_adm${i}_${INTERVAL_DISTANCE}.geojson + + mapshaper "$FILE" -simplify interval=$INTERVAL_DISTANCE -o "$OFILE" + + SIZE=$(stat --printf="%s" "$FILE") + OSIZE=$(stat --printf="%s" "$OFILE") + + echo " -> Finished file size went from $SIZE to $OSIZE" + fi + psql -q -d $POSTGRES_DB -U $POSTGRES_USER -h $POSTGRES_HOST -p $POSTGRES_PORT << SQL_SCRIPT - WITH data AS (SELECT \$$`cat $FILE`\$$::json AS fc) + WITH data AS (SELECT \$$`cat "$OFILE"`\$$::json AS fc) INSERT INTO "regions" (country, name, admin_level, the_geom) - SELECT - '${COUNTRY_NAME}', - feat#>>'{properties,name}' AS name, - ${i}, - ST_SetSRID(ST_CollectionExtract(ST_Multi(ST_GeomFromGeoJSON(feat->>'geometry')), 3), 4326) as the_geom - FROM ( - SELECT json_array_elements(fc->'features') AS feat - FROM data - ) AS f; + SELECT + '${COUNTRY_NAME}', + feat#>>'{properties,name}' AS name, + ${i}, + ST_SetSRID(ST_CollectionExtract(ST_Multi(ST_GeomFromGeoJSON(feat->>'geometry')), 3), 4326) as the_geom + FROM ( + SELECT json_array_elements(fc->'features') AS feat + FROM data + ) AS f + ON CONFLICT ON CONSTRAINT regions_name_index DO + UPDATE SET the_geom=excluded.the_geom; SQL_SCRIPT + # TODO: print errors echo " done!" done;