From 1f271c5cbd06e672932bc2edc33dd92b04f2e803 Mon Sep 17 00:00:00 2001 From: Ismael Bejarano Date: Thu, 30 May 2024 21:08:01 -0300 Subject: [PATCH 1/4] Use mapshaper to simplify a region geometry --- scripts/friction/load-friction-raster | 1 + scripts/population/load-regions | 19 ++++++++++++++++++- 2 files changed, 19 insertions(+), 1 deletion(-) diff --git a/scripts/friction/load-friction-raster b/scripts/friction/load-friction-raster index 542955af..20bc3e50 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/load-regions b/scripts/population/load-regions index a8f8e093..06c64608 100755 --- a/scripts/population/load-regions +++ b/scripts/population/load-regions @@ -33,13 +33,30 @@ 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 (see planar option below), +# otherwise uses the same units as the source data. +INTERVAL_DISTANCE=${INTERVAL_DISTANCE-1000} + 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 + 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}', From 2a7bfac0c804a275e6dd07e8c8608532b8702058 Mon Sep 17 00:00:00 2001 From: Ismael Bejarano Date: Fri, 31 May 2024 10:10:23 -0300 Subject: [PATCH 2/4] * Document changes in README * Add -f to force geometry updates * Print sizes of files --- scripts/population/README.md | 14 +++++- scripts/population/load-regions | 88 ++++++++++++++++++++++++++------- 2 files changed, 83 insertions(+), 19 deletions(-) diff --git a/scripts/population/README.md b/scripts/population/README.md index 8d5e38db..19f5c23a 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 06c64608..ea1ea91b 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,7 +41,7 @@ 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 @@ -37,9 +60,12 @@ fi # 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 (see planar option below), +# Uses meters when simplifying unprojected datasets in 3D space, # otherwise uses the same units as the source data. -INTERVAL_DISTANCE=${INTERVAL_DISTANCE-1000} +# 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 @@ -52,23 +78,49 @@ for i in $LEVELS; do 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 + if [ "$FORCE" = "1" ]; then + psql -q -d $POSTGRES_DB -U $POSTGRES_USER -h $POSTGRES_HOST -p $POSTGRES_PORT << SQL_SCRIPT + + WITH data AS (SELECT \$$`cat $OFILE`\$$::json AS fc) + UPDATE "regions" SET the_geom=f.the_geom + FROM ( + SELECT + '${COUNTRY_NAME}' as country, + feat#>>'{properties,name}' AS name, + ${i} as admin_level, + 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 dummy + ) AS f + WHERE regions.country=f.country AND regions.name=f.name AND regions.admin_level=f.admin_level; - 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; +SQL_SCRIPT + else + psql -q -d $POSTGRES_DB -U $POSTGRES_USER -h $POSTGRES_HOST -p $POSTGRES_PORT << SQL_SCRIPT + + 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; SQL_SCRIPT + fi # TODO: print errors echo " done!" done; From a54d4f530b992855bf9bfb154928b433b8c31fd2 Mon Sep 17 00:00:00 2001 From: Ismael Bejarano Date: Fri, 31 May 2024 10:40:06 -0300 Subject: [PATCH 3/4] Add postgis and pgrouting --- resources/planwise/plpgsql/isochrones.sql | 3 +++ resources/planwise/plpgsql/patches.sql | 3 +++ 2 files changed, 6 insertions(+) diff --git a/resources/planwise/plpgsql/isochrones.sql b/resources/planwise/plpgsql/isochrones.sql index 79d53050..0c68616d 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 da46f639..d77a5d4f 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 From 0aa5e8dde88c6bf7a1801a29f5a23fa26c0efb2f Mon Sep 17 00:00:00 2001 From: Ismael Bejarano Date: Wed, 12 Jun 2024 18:04:34 -0300 Subject: [PATCH 4/4] * 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 + scripts/population/load-regions | 40 +++++-------------- 3 files changed, 13 insertions(+), 29 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 00000000..3c38705d --- /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 00000000..e4f3df33 --- /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/scripts/population/load-regions b/scripts/population/load-regions index ea1ea91b..8ea19bb1 100755 --- a/scripts/population/load-regions +++ b/scripts/population/load-regions @@ -48,7 +48,7 @@ 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 @@ -77,38 +77,18 @@ for i in $LEVELS; do # Simplify using distance OFILE=${DATA_PATH}/${COUNTRY_CODE}/${COUNTRY_CODE}_adm${i}_${INTERVAL_DISTANCE}.geojson - mapshaper $FILE -simplify interval=$INTERVAL_DISTANCE -o $OFILE + mapshaper "$FILE" -simplify interval=$INTERVAL_DISTANCE -o "$OFILE" - SIZE=$(stat --printf="%s" $FILE) - OSIZE=$(stat --printf="%s" $OFILE) + SIZE=$(stat --printf="%s" "$FILE") + OSIZE=$(stat --printf="%s" "$OFILE") echo " -> Finished file size went from $SIZE to $OSIZE" fi - if [ "$FORCE" = "1" ]; then - psql -q -d $POSTGRES_DB -U $POSTGRES_USER -h $POSTGRES_HOST -p $POSTGRES_PORT << SQL_SCRIPT + psql -q -d $POSTGRES_DB -U $POSTGRES_USER -h $POSTGRES_HOST -p $POSTGRES_PORT << SQL_SCRIPT - WITH data AS (SELECT \$$`cat $OFILE`\$$::json AS fc) - UPDATE "regions" SET the_geom=f.the_geom - FROM ( - SELECT - '${COUNTRY_NAME}' as country, - feat#>>'{properties,name}' AS name, - ${i} as admin_level, - 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 dummy - ) AS f - WHERE regions.country=f.country AND regions.name=f.name AND regions.admin_level=f.admin_level; - -SQL_SCRIPT - else - psql -q -d $POSTGRES_DB -U $POSTGRES_USER -h $POSTGRES_HOST -p $POSTGRES_PORT << SQL_SCRIPT - - WITH data AS (SELECT \$$`cat $OFILE`\$$::json AS fc) - INSERT INTO "regions" (country, name, admin_level, the_geom) + 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, @@ -117,10 +97,12 @@ SQL_SCRIPT FROM ( SELECT json_array_elements(fc->'features') AS feat FROM data - ) AS f; + ) AS f + ON CONFLICT ON CONSTRAINT regions_name_index DO + UPDATE SET the_geom=excluded.the_geom; SQL_SCRIPT - fi + # TODO: print errors echo " done!" done;