Topology Simplification

00_Topology Simplification.md

Making Use of Topology

This runs on a restored dird_db from a Docker container, postgis:postgis/14-3.2; I run this as docker run -P ghcr.io/joshbrooks/openly_dird_prod

Pervious versions of postgis run into trouble on the topology building

First create the SimplifyEdgeGeom function: (https://gist.github.com/leplatrem/5729022) This wraps ST_ChangeEdgeGeom in error handling in cases of ugliness

01_SimplifyEdgeGeom.sql

Generate the source table. Here we're using the intersection of "area" and "altitude".

02_source_table.sql

Generate the topology requirements

03_topology_and_table.sql

Simplify the table, and replace the records in altitude_for_area, and truncate the altitude` table which we no longer require

05_finishing_touches.sql

01_SimplifyEdgeGeom.sql

-- This wraps `topology.ST_ChangeEdgeGeom`
-- with a retry-different-tolerance
CREATE OR REPLACE FUNCTION SimplifyEdgeGeom(atopo varchar, anedge int, maxtolerance float8)
RETURNS float8 AS $$
DECLARE
  tol float8;
  sql varchar;
BEGIN
  tol := maxtolerance;
  LOOP
    sql := 'SELECT topology.ST_ChangeEdgeGeom(' || quote_literal(atopo) || ', ' || anedge
      || ', ST_Simplify(geom, ' || tol || ')) FROM '
      || quote_ident(atopo) || '.edge WHERE edge_id = ' || anedge;
    BEGIN
      RAISE DEBUG 'Running %', sql;
      EXECUTE sql;
      RETURN tol;
    EXCEPTION
     WHEN OTHERS THEN
      RAISE WARNING 'Simplification of edge % with tolerance % failed: %', anedge, tol, SQLERRM;
      tol := round( (tol/2.0) * 1e8 ) / 1e8; -- round to get to zero quicker
      IF tol = 0 THEN raise
     --EXCEPTION 
     notice '%', SQLERRM; END IF;
    END;
  END LOOP;
end
$$ LANGUAGE 'plpgsql' STABLE STRICT;

02_source_table.sql

create table intersect_area_altitude(id serial, area_id int, altitude_level text, geom geometry(polygon, 4326));
WITH intersection AS (SELECT
    "area"."id" "area_id", "altitude"."altitude_level",
        ST_MULTI(ST_INTERSECTION("altitude"."geom", "area"."geom")) AS "geom"
    FROM simple_locations_area area, location_profile_altitude altitude, simple_locations_areatype areatype
        WHERE ST_INTERSECTS("altitude"."geom", "area"."geom")
        AND "area"."kind_id" = "areatype"."id"
        AND "areatype"."slug" = 'district'
    ),
    dump_geometries AS (
        SELECT "area_id", "altitude_level", (ST_DUMP(geom)).geom AS geom FROM intersection
)
INSERT INTO intersect_area_altitude (area_id, altitude_level, geom)
select area_id, altitude_level, geom from dump_geometries;

03_topology_and_table.sql

-- Create a temporary topology
-- and suitable table 
select topology.DropTopology('altitude_area_topology');
SELECT topology.CreateTopology('altitude_area_topology', 32756);
drop table if exists altitude_topology;
CREATE TABLE if not exists altitude_topology(intersect_area_altitude_id int, area_id int, altitude_level text);
-- The layer ID here should be 1
SELECT topology.AddTopoGeometryColumn('altitude_area_topology', 'public', 'altitude_topology', 'topo', 'POLYGON');

04_add_to_topology_table.sql

-- Loop over the intersected values
-- and try to add them to the topology table
DO $$DECLARE r record;
BEGIN
  FOR r IN SELECT * FROM intersect_area_altitude a where 
  a.id not in (select intersect_area_altitude_id from altitude_topology)
  LOOP
    begin
		INSERT INTO altitude_topology (intersect_area_altitude_id, area_id, altitude_level, topo)
		select
			id,
			area_id,
			altitude_level,
			topology.toTopoGeom(ST_Transform(geom, 32756), 'altitude_area_topology', 1, 10) from intersect_area_altitude
		where
			intersect_area_altitude.id = r.id;
    EXCEPTION
      WHEN OTHERS THEN
        RAISE WARNING 'Loading of record % failed: %', r.id, SQLERRM;
    END;
  END LOOP;
END$$;

04b_update_without_loop.sql

-- This should work if your postgis version is 3.2+
-- This updated 2163 rows in 1:37
INSERT INTO altitude_topology (intersect_area_altitude_id, area_id, altitude_level, topo)
		select
			id,
			area_id,
			altitude_level,
			topology.toTopoGeom(ST_Transform(geom, 32756), 'altitude_area_topology', 1, 10) from intersect_area_altitude

05_finishing_touches.sql

select SUM(SimplifyEdgeGeom('altitude_area_topology', edge_id, 200)) from altitude_area_topology.edge;
drop table if exists topology_altitude_200;
create table topology_altitude_200 as select area_id, altitude_level, topo::geometry from altitude_topology;


-- Replace the altitudeforarea records with the simplified ones
truncate location_profile_altitudeforarea;
insert into location_profile_altitudeforarea (altitude_level, geom, area_id)
select altitude_level, (ST_DUMP(ST_TRANSFORM(topo, 4326))).geom, area_id from topology_altitude_200;
-- We no longer need this table
truncate location_profile_altitude;
-- This zeroes out decimal points. No effect on table size on server, but 
-- reduces size for on-disk representation.
update location_profile_altitudeforarea set geom = ST_QuantizeCoordinates(geom, 3);

06_export_update_table.md

docker exec -it gracious_morse\
  pg_dump\
    -f location_profile_altitudeforarea.sql\
    -a\
    -t location_profile_altitudeforarea \
    --user dird dird_db

Fetch the output sql script

docker cp gracious_morse:location_profile_altitudeforarea.sql .

From your manage.py truncate the current tables

./manage.py dbshell

TRUNCATE location_profile_altitude;
TRUNCATE location_profile_altitudeforarea;

./manage.py dbshell < location_profile_altitudeforarea.sql

07_complete_script_for_rainfall

-- This wraps `topology.ST_ChangeEdgeGeom`
-- with a retry-different-tolerance
CREATE OR REPLACE FUNCTION SimplifyEdgeGeom(atopo varchar, anedge int, maxtolerance float8)
RETURNS float8 AS $$
DECLARE
  tol float8;
  sql varchar;
BEGIN
  tol := maxtolerance;
  LOOP
    sql := 'SELECT topology.ST_ChangeEdgeGeom(' || quote_literal(atopo) || ', ' || anedge
      || ', ST_Simplify(geom, ' || tol || ')) FROM '
      || quote_ident(atopo) || '.edge WHERE edge_id = ' || anedge;
    BEGIN
      RAISE DEBUG 'Running %', sql;
      EXECUTE sql;
      RETURN tol;
    EXCEPTION
     WHEN OTHERS THEN
      RAISE WARNING 'Simplification of edge % with tolerance % failed: %', anedge, tol, SQLERRM;
      tol := round( (tol/2.0) * 1e8 ) / 1e8; -- round to get to zero quicker
      IF tol = 0 THEN raise
     --EXCEPTION 
     notice '%', SQLERRM; END IF;
    END;
  END LOOP;
end
$$ LANGUAGE 'plpgsql' STABLE STRICT;

create table intersect_area_rainfall(id serial, area_id int, rainfall text, geom geometry(polygon, 4326));
WITH intersection AS (SELECT
    "area"."id" "area_id", "rainfall"."rainfall",
        ST_MULTI(ST_INTERSECTION("rainfall"."geom", "area"."geom")) AS "geom"
    FROM simple_locations_area area, location_profile_rainfall rainfall, simple_locations_areatype areatype
        WHERE ST_INTERSECTS("rainfall"."geom", "area"."geom")
        AND "area"."kind_id" = "areatype"."id"
        AND "areatype"."slug" = 'district'
    ),
    dump_geometries AS (
        SELECT "area_id", "rainfall", (ST_DUMP(geom)).geom AS geom FROM intersection
)
INSERT INTO intersect_area_rainfall (area_id, rainfall, geom)
select area_id, rainfall, geom from dump_geometries;

-- Create a temporary topology
-- and suitable table 

SELECT topology.CreateTopology('rainfall_area_topology', 32756);
drop table if exists rainfall_topology;
CREATE TABLE if not exists rainfall_topology(intersect_area_rainfall_id int, area_id int, rainfall text);
-- The layer ID here should be 1
SELECT topology.AddTopoGeometryColumn('rainfall_area_topology', 'public', 'rainfall_topology', 'topo', 'POLYGON');


-- This should work if your postgis version is 3.2+

INSERT INTO rainfall_topology (intersect_area_altitude_id, area_id, rainfall, topo)
		select
			id,
			area_id,
			rainfall,
			topology.toTopoGeom(ST_Transform(geom, 32756), 'rainfall_area_topology', 1, 10) from intersect_area_altitude



select SUM(SimplifyEdgeGeom('rainfall_area_topology', edge_id, 200)) from rainfall_area_topology.edge;
drop table if exists topology_rainfall_200;
create table topology_rainfall_200 as select area_id, rainfall, topo::geometry from rainfall_topology;


-- Replace the rainfallforarea records with the simplified ones
truncate location_profile_rainfallforarea;
insert into location_profile_rainfallforarea (rainfall, geom, area_id)
select rainfall, (ST_DUMP(ST_TRANSFORM(topo, 4326))).geom, area_id from topology_rainfall_200;
-- We no longer need this table
truncate location_profile_rainfall;
-- This zeroes out decimal points. No effect on table size on server, but 
-- reduces size for on-disk representation.
update location_profile_rainfallforarea set geom = ST_QuantizeCoordinates(geom, 3);



-- Clean up

--select topology.DropTopology('rainfall_area_topology');

08_complete_script_for_area.sql

SELECT topology.CreateTopology('area_topo', 32756);
drop table if exists area_topo_table;
CREATE TABLE if not exists area_topo_table(area_id int);
-- The layer ID here should be 1
SELECT topology.AddTopoGeometryColumn('area_topo', 'public', 'area_topo_table', 'topo', 'POLYGON');

INSERT INTO area_topo_table (area_id, topo)
    select
        id,
        topology.toTopoGeom(ST_Transform(geom, 32756), 'area_topo', 1, 10) from simple_locations_area;
       select * from area_topo_table

select SUM(SimplifyEdgeGeom('area_topo', edge_id, 100)) from area_topo.edge;
update simple_locations_area set geom = (select st_transform(topo::geometry, 4326) from area_topo_table
where area_geom_replace.area_id = simple_locations_area.id);
vacuum full verbose analyze simple_locations_area;

--SELECT topology.DropTopology('area_topo');

create table area_geom_replace as 
	select area_id,
	st_quantizeCoordinates(
		st_transform(topo::geometry, 4326),
		2
		)
 from  "area_topo_table";


--copy the sql to your ./manage.py directory

--docker exec {container name} pg_dump -t area_geom_replace --user dird dird_db -f area_geom_replace.sql
--docker cp {container name}:area_geom_replace.sql .
--./manage.py dbshell < area_geom_replace.sql

--then run this SQL from `./manage.py dbshell`
--UPDATE simple_locations_area a SET geom = (SELECT geom FROM area_geom_replace b WHERE id = area_id);
--DROP TABLE area_geom_replace;

-- then remove the sql
-- rm area_geom_replace.sql