[postgis-users] Topology: Simplifying Coastlines with Islands

[Michiel J. van Heek]

Quoting Sandro Santilli <strk at keybit.net>:

> On Sat, May 26, 2012 at 10:37:14PM +0200, Michiel J. van Heek wrote:
>
>> So, why not use ST_SimplifyPreserveTopology by default? Or is that
>> function generally a lot slower than ST_Simplify?
>
> Is a lot slower.
> But the script was just an example, you can tweak the heuristic to
> parse the error message and use *PreserveTopology for the 'non-simple'
> case, rather than reduce the tolerance. Also, you may want to reduce
> less than 1/2. You could have got away with a 0.45. Finally you could
> just let that edge go in the first iteration and run a second time.
> Some times just letting the function run for other edges automatically
> removes the obstacles for next iteration...

Sandro,

Based on your suggestions, I changed the script in three ways (see  
below for the result):

1. The number of iterations is now a parameter and running down the  
tolerance uses linear steps for each iteration (0.5, 0.45, 0.4, 0.35,  
etc.). The division by 2 solution (0.5, 0.25, 0.125, etc.) put a heavy  
weight on the lower tolerances, whereas it is often worthwile to also  
try higher tolerances (between initial tolerance and half of the  
initial tolerance).

2. If an iteration fails because of "not simple" the same iteration is  
run again, but with "PreserveTopology" now.

3. When all iterations fail, the function returns NULL now, and does  
not throw an exception. It is now possible to find "unwilling" edges  
by running:

SELECT *
FROM (
     SELECT
         edge_id,
         SimplifyEdgeGeom('countries_topology_4', edge_id, 0.5) AS tolerance
     FROM countries_topology_4.edge
) AS foo
WHERE (tolerance IS NULL);

Also, I removed the STABLE attribute, because the funtion *does* make  
changes to the database.

Michiel


CREATE OR REPLACE FUNCTION SimplifyEdgeGeom(topology_param text,  
edge_id_param integer, tolerance_param double precision,  
num_iterations_param integer DEFAULT 10) RETURNS double precision AS  
$BD$
DECLARE
     i integer := 1;
     xi text := '';
     tolerance_var double precision := tolerance_param;
     sql_var text;
BEGIN
     WHILE (i <= num_iterations_param) LOOP
         sql_var := 'SELECT topology.ST_ChangeEdgeGeom(' ||  
quote_literal(topology_param) || ', ' || edge_id_param
           || ', ST_Simplify' || xi || '(geom, ' || tolerance_var || ')) FROM '
           || quote_ident(topology_param) || '.edge WHERE edge_id = '  
|| edge_id_param;
         BEGIN
             RAISE DEBUG 'Running %', sql_var;
             EXECUTE sql_var;

             RETURN tolerance_var;
         EXCEPTION
             WHEN OTHERS THEN
                 RAISE WARNING 'Simplification of edge % failed in  
iteration %_%, with tolerance %: %',
                     edge_id_param, i, xi, tolerance_var, SQLERRM;

                 IF (SQLERRM = 'SQL/MM Spatial exception - curve not  
simple') AND (xi = '') THEN
                     xi := 'PreserveTopology';
                 ELSE
                     i := i + 1;
                     xi := '';
                     tolerance_var := tolerance_var - (tolerance_param  
/ num_iterations_param);
                 END IF;
             END;
         -- END EXCEPTION
     END LOOP;

     RETURN NULL;
END;
$BD$ LANGUAGE 'plpgsql' STRICT;