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

resources/migrations/068-add_regions_name_index.down.sql

@@ -0,0 +1 @@
+ALTER TABLE regions DROP CONSTRAINT "regions_name_index";

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");

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 $$

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

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

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.
-

scripts/population/load-regions

@@ -4,12 +4,35 @@ set -euo pipefail
 export PGPASSWORD=$POSTGRES_PASSWORD;
 
 if [ $# -lt 2 ]; then
-    echo "Usage $0 <ISO> <Country Name>"
+    echo "Usage $0 [-f] <ISO> <Country Name>"
+    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,40 +41,68 @@ 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
   echo "   -> GeoJSON folder not found for country $COUNTRY_NAME"
   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;