[postgis-users] totopogeom resulting in SQL/MM Spatial exception

[Lars Aksel Opsahl]

>From: postgis-users <postgis-users-bounces at lists.osgeo.org> on behalf of Bo Guo <bo.guo at gisticinc.com>
>Sent: Wednesday, January 16, 2019 7:00 PM
>To: postgis-users at lists.osgeo.org
>Subject: Re: [postgis-users] totopogeom resulting in SQL/MM Spatial exception
>
>OK, I can understand the concept why this is the case now. But if there
>is any example or reference article on line that would be of great help.
>
>Bo
>
>On 1/16/19 10:37 AM, Sandro Santilli wrote:
>> On Wed, Jan 16, 2019 at 10:26:12AM -0700, Bo Guo wrote:
>>> strk,
>>>
>>> Here is how it break up the topo geom conversion (l_tolerance is 0.0000001)
>>>
>>>              LOOP
>>>                  WITH foo AS (SELECT grd_id FROM azgiv.roadcenterlines
>>>                                  WHERE grd_topo_geom IS NULL
>>>                                      LIMIT l_batch_size)
>>>                  UPDATE azgiv.roadcenterlines
>>>                                  SET grd_topo_geom =
>>> topology.totopogeom(grd_geom, 'azgiv_topo', l_topo_layer_id, l_tolerance)
>>>                                  FROM foo
>>>                                  WHERE foo.grd_id = roadcenterlines.grd_id;
>>>
>>>                  GET DIAGNOSTICS l_rowcount = ROW_COUNT;
>>>
>>>                  EXIT WHEN l_rowcount < l_batch_size;
>>>
>>>              END LOOP;
>> The goal of chunking was to get partial results rather than an
>> all-or-nothing behavior. If you use that plpgsql loop you'll want
>> to catch exceptions and set those TopoGeometries to NULL instead.
>> Then you'll be able to see what the loop was able to convert and
>> what not, and get back to the still-to-be-converted geoms, maybe
>> after cleaning up some of what you got converted already.
>>
>> --strk;


Hi


This is not a article but a short description on how we last year at Nibio used Postgis Topology to clean up messy vector data with a lot of overlap/gap. The test we did contained about 1.3 mill polygons with a total of around 90 mill points.


The input was simple feature polygons that we converted into Postgis Topology and then we then did line smoothing and point reduction on the border lines . There is attached a picture where the red lines are input and the green ones are an example of output at one stage.


The Postgis Topology worked quite good and better than other methods we tested. Yes we got some errors, but we also got very good help from Sandro and from examples on internet to resolve the problems we had, Thanks a lot. (When we where able to make a small test cases that failed it was certainly easier fix/identify problem and get help).


So here is a short description :


1) First we dived the input in many (more than a hundred cells) small cells by using a content based grid function .


2) Then we ran parallel processes on the cells. In each parallel process we did the following :

  *   Created a temporary topology layer with all the line strings from the polygons inside the cell.

  *   Then we cut the lines inside the cell so no lines was closer then one meter to the cell walls.

  *   Made a final result for each cell, where we snapped lines and smoothed them.

  *   Then we dumped this result into a master Postgis Topology table for all the cell.

  *   The lines parts in two the meter buffer zone on the cell border was just saved in a temporary table for later processing.


3) When all the parallel processes was done we had a master table where no cell had any connecting lines to another cell, because all the lines made in phase two was cut one meter from each cell border.

  *   When then started to connect the different cells using the border lines that was saved in the temporary table. The border lines had the original form and the point where line was inside cell had not been moved so it was possible to match them together. This lines was also quite short since we did a intersection on them before the where put into the temp table.

  *   Then we connected the border lines with  original lines and this was the more tricky part, but it seemed to work ok. This connection line part is the most unsure part of this  operation.

  *   Then we simplified the border lines and restored the missing cell border line polygons.


4) Then finally we create the simple feature polygons based on the Postgis Topology master polygon table where all cell now are connected.


This did actually work. As Sandro say the clue is to dived work into smaller pieces and do try catch. We will continue work more on this year hopefully.


I have added a example of some code below in this mail on how we loop and break line into single pieces and added each part separately as a last try. If it failed a the lowest level we logged the failing piece to to a log table ,so we have total control of failing lines.


If we get the time we may try to put out the code on a repo, but some parts of the util code that we depend on is out on github.


A) This repo contains code used by a simple client to edit Postgis Topology, but this is also used to by the Simple Feature cleaner described above. Here we also good from help Sandro.

(Some of this code is not generic and need much more work it if should be used by others)

https://github.com/NibioOpenSource/pgtopo_update_sql


At https://www.slideshare.net/laopsahl/postgis-topology-presentationfoss4goslo01092016

there is presentation about why we use Topology and how we use Postgis Topology at Nibio for online map editing.


B) This contains code for content based grids to be able work big datasets

https://github.com/larsop/content_balanced_grid .


If you look at page 18-23 at http://www.slideshare.net/laopsahl/foss4-g-topologyjuly162015

you find some info about content based grids.


C) Repo with different extended Chaikin smooth line methods written in psql used to smooth lines https://github.com/larsop/chaikin when cleaning up messy polygons in the process above.



Example code


CREATE OR REPLACE FUNCTION

topo_ar5_forest.create_nocutline_edge_domain_obj_retry(

json_feature text,

border_topo_info topo_update.input_meta_info,

server_json_feature text default null)

RETURNS TABLE(id integer) AS $$

DECLARE



-- holds dynamic sql to be able to use the same code for different

command_string text;


-- holde the computed value for json input reday to use

json_input_structure topo_update.json_input_structure;


-- holde the computed value for json input reday to use

json_input_structure_tmp topo_update.json_input_structure;


counter integer:=0;

g Geometry;

inGeom Geometry;

tolerance real = 10;

start_tolerance integer = 10;


start_time timestamp with time zone;

done_time timestamp with time zone;

used_time real;

apoint geometry[] DEFAULT '{}';

spltt_line geometry;

failed_to_insert boolean ;

feat json;


rec record;

BEGIN


start_time := clock_timestamp();


feat := json_feature::json;

json_input_structure.input_geo := ST_GeomFromGeoJSON(feat->>'geometry');


start_tolerance = border_topo_info.snap_tolerance;


BEGIN


RAISE NOTICE 'work start with %, containing % points', json_input_structure.input_geo, ST_NumPoints(json_input_structure.input_geo);



--test remove reptead points

inGeom := ST_RemoveRepeatedPoints(json_input_structure.input_geo,(start_tolerance*3));

json_input_structure.input_geo := inGeom;


IF ST_NumPoints(json_input_structure.input_geo) < 1000 THEN

perform topo_ar5_forest.create_nocutline_edge_domain_try_one( border_topo_info, json_input_structure);

ELSE

RAISE NOTICE 'work start, to big:% border_layer_id %, with a line containing % points', start_time, border_topo_info.border_layer_id, ST_NumPoints(json_input_structure.input_geo);

json_input_structure_tmp := topo_update.handle_input_json_props(json_feature::json,server_json_feature::json,border_topo_info.srid);

inGeom := json_input_structure.input_geo;

LOOP

counter := counter + 1000;

-- some computations

IF counter > (ST_NPoints(inGeom)-1) THEN

EXIT; -- exit loop

ELSE

apoint := array_append(apoint, ST_PointN(inGeom,counter));

--ST_PointN(g,counter);

END IF;

END LOOP;

spltt_line := ST_Split(inGeom,ST_Multi(ST_Collect(apoint)));

drop table if exists line_list_tmp;

create temp table line_list_tmp as (select (ST_Dump(spltt_line)).geom AS line_part);

FOR rec IN

SELECT *

FROM line_list_tmp

LOOP

RAISE NOTICE 'rec %', ST_Length(rec.line_part);

json_input_structure_tmp.input_geo = rec.line_part;

perform topo_ar5_forest.create_nocutline_edge_domain_try_one( border_topo_info, json_input_structure_tmp);

END LOOP;

END IF;


done_time := clock_timestamp();

EXCEPTION WHEN OTHERS THEN

RAISE NOTICE 'failed ::::::1 %', border_topo_info.border_layer_id;

counter := 0;


json_input_structure_tmp := topo_update.handle_input_json_props(json_feature::json,server_json_feature::json,border_topo_info.srid);

inGeom := json_input_structure.input_geo;

FOR i IN 1..(ST_NPoints(inGeom)-1)

LOOP

border_topo_info.snap_tolerance := start_tolerance;


counter:=counter+1;

g := ST_MakeLine(ST_PointN(inGeom,i),ST_PointN(inGeom,i+1));

perform ST_setSrid(g,border_topo_info.srid);

BEGIN

json_input_structure_tmp.input_geo = g;

perform topo_ar5_forest.create_nocutline_edge_domain_try_one( border_topo_info, json_input_structure_tmp);

-- catch EXCEPTION

EXCEPTION WHEN OTHERS THEN

RAISE NOTICE 'failed ::::::2 % num %, border_topo_info.snap_tolerance %', border_topo_info.border_layer_id, i, border_topo_info.snap_tolerance;


tolerance := border_topo_info.snap_tolerance;


WHILE tolerance > 0 LOOP

failed_to_insert := false;

IF tolerance = 1 THEN

tolerance := 0.01;

ELSE

tolerance := tolerance - 1;

END IF;


IF tolerance < 0 THEN

tolerance := 0;

END IF;

BEGIN

border_topo_info.snap_tolerance = tolerance ;

perform topo_ar5_forest.create_nocutline_edge_domain_try_one( border_topo_info, json_input_structure_tmp);

exit;

EXCEPTION WHEN OTHERS THEN

RAISE NOTICE 'failed with with % : %', ST_AsText(g), tolerance;

failed_to_insert := true;

END;

END LOOP;

IF failed_to_insert THEN

done_time := clock_timestamp();

used_time := (EXTRACT(EPOCH FROM (done_time - start_time)));

RAISE NOTICE 'ERROR failed to use %, length: % tolerance : %', ST_AsText(g), ST_length(g), tolerance;

insert into topo_ar5_forest.no_cut_line_failed(error_info,geo)

values('Failed with exception time used ' || used_time::varchar || ' length ' || ST_length(g), g);

END IF;

END;

END LOOP;


END;


done_time := clock_timestamp();

used_time := (EXTRACT(EPOCH FROM (done_time - start_time)));


IF used_time > 10 THEN

RAISE NOTICE 'very long single line % time with geo for % ', used_time, json_input_structure.input_geo;

insert into topo_ar5_forest.long_time_log1(execute_time,info,geo)

values(used_time,'long ' || used_time::varchar || ' num points ' || ST_NumPoints(json_input_structure.input_geo), json_input_structure.input_geo);

END IF;

return;


END;

$$ LANGUAGE plpgsql;


Lars

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20190118/a9a57d02/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: segm_topo4_smaapolydetalj.png
Type: image/png
Size: 38935 bytes
Desc: segm_topo4_smaapolydetalj.png
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20190118/a9a57d02/attachment.png>