From 9ab51027fcdbae756d250e6120838098c60dedf4 Mon Sep 17 00:00:00 2001 From: abelvm Date: Wed, 18 Jan 2017 17:28:06 +0100 Subject: [PATCH 1/9] support multi --- src/pg/sql/13_PIA.sql | 34 +++++++++++++++++++++++++--------- 1 file changed, 25 insertions(+), 9 deletions(-) diff --git a/src/pg/sql/13_PIA.sql b/src/pg/sql/13_PIA.sql index d9a224d..d6caa10 100644 --- a/src/pg/sql/13_PIA.sql +++ b/src/pg/sql/13_PIA.sql @@ -96,27 +96,43 @@ $$ language plpgsql IMMUTABLE; -- signed distance point to polygon with holes -- negative is the point is out the polygon +-- rev 1. adding MULTIPOLYGON and GEOMETRYCOLLECTION support by @abelvm CREATE OR REPLACE FUNCTION _Signed_Dist( IN polygon geometry, IN point geometry ) RETURNS numeric AS $$ DECLARE + pols geometry[]; + pol geometry; i integer; + j integer; within integer; + w integer; holes integer; dist numeric; + d numeric; BEGIN dist := 1e999; - SELECT LEAST(dist, ST_distance(point, ST_ExteriorRing(polygon))::numeric) INTO dist; - SELECT CASE WHEN ST_Within(point,polygon) THEN 1 ELSE -1 END INTO within; - SELECT ST_NumInteriorRings(polygon) INTO holes; - IF holes > 0 THEN - FOR i IN 1..holes - LOOP - SELECT LEAST(dist, ST_distance(point, ST_InteriorRingN(polygon, i))::numeric) INTO dist; - END LOOP; - END IF; + pols := array_agg((ST_dump(polygon)).geom); + FOR j in 1..array_length(pols, 1); + LOOP + pol := pols[j]; + d := dist; + SELECT LEAST(dist, ST_distance(point, ST_ExteriorRing(pol))::numeric) INTO d; + SELECT CASE WHEN ST_Within(point,pol) THEN 1 ELSE -1 END INTO w; + SELECT ST_NumInteriorRings(pol) INTO holes; + IF holes > 0 THEN + FOR i IN 1..holes + LOOP + SELECT LEAST(d, ST_distance(point, ST_InteriorRingN(pol, i))::numeric) INTO d; + END LOOP; + END IF; + IF d < dist THEN + dist:= d; + within := w; + END IF; + END LOOP; dist := dist * within::numeric; RETURN dist; END; From e03c3eece2cfe2589759999d7412f8ef9cb9a52a Mon Sep 17 00:00:00 2001 From: abelvm Date: Wed, 18 Jan 2017 21:32:42 +0100 Subject: [PATCH 2/9] semi colon fix --- src/pg/sql/13_PIA.sql | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pg/sql/13_PIA.sql b/src/pg/sql/13_PIA.sql index d6caa10..ed3c76b 100644 --- a/src/pg/sql/13_PIA.sql +++ b/src/pg/sql/13_PIA.sql @@ -115,7 +115,7 @@ DECLARE BEGIN dist := 1e999; pols := array_agg((ST_dump(polygon)).geom); - FOR j in 1..array_length(pols, 1); + FOR j in 1..array_length(pols, 1) LOOP pol := pols[j]; d := dist; From d8604f3c9b6ab9d2c654a310919dad913ba66edd Mon Sep 17 00:00:00 2001 From: abelvm Date: Wed, 18 Jan 2017 21:40:36 +0100 Subject: [PATCH 3/9] agg set error fix --- src/pg/sql/13_PIA.sql | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pg/sql/13_PIA.sql b/src/pg/sql/13_PIA.sql index ed3c76b..2b90a95 100644 --- a/src/pg/sql/13_PIA.sql +++ b/src/pg/sql/13_PIA.sql @@ -114,7 +114,7 @@ DECLARE d numeric; BEGIN dist := 1e999; - pols := array_agg((ST_dump(polygon)).geom); + WITH collection as (SELECT (ST_dump(polygon)).geom as geom) SELECT array_agg(geom) into pols FROM collection; FOR j in 1..array_length(pols, 1) LOOP pol := pols[j]; From 47251daa5fb7c59406c77ae8be3a33dd5b2fd657 Mon Sep 17 00:00:00 2001 From: abelvm Date: Tue, 28 Mar 2017 13:02:37 +0200 Subject: [PATCH 4/9] fixed corner case centroid=PIA --- src/pg/sql/13_PIA.sql | 1 + 1 file changed, 1 insertion(+) diff --git a/src/pg/sql/13_PIA.sql b/src/pg/sql/13_PIA.sql index 2b90a95..bed2593 100644 --- a/src/pg/sql/13_PIA.sql +++ b/src/pg/sql/13_PIA.sql @@ -46,6 +46,7 @@ BEGIN SELECT array_agg(c) INTO cells FROM c1; -- 1st guess: centroid + best_c := polygon; best_d := cdb_crankshaft._Signed_Dist(polygon, ST_Centroid(Polygon)); -- looping the loop From c252c18adcd84c5e6f6225c4c314b833272d3092 Mon Sep 17 00:00:00 2001 From: abelvm Date: Tue, 28 Mar 2017 13:32:10 +0200 Subject: [PATCH 5/9] fixed indexing --- src/pg/sql/13_PIA.sql | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/pg/sql/13_PIA.sql b/src/pg/sql/13_PIA.sql index bed2593..38963fa 100644 --- a/src/pg/sql/13_PIA.sql +++ b/src/pg/sql/13_PIA.sql @@ -31,7 +31,7 @@ DECLARE sqr numeric; p geometry; BEGIN - sqr := |/2; + sqr := 0.5*(|/2.0); polygon := ST_Transform(polygon, 3857); -- grid #0 cell size @@ -46,7 +46,7 @@ BEGIN SELECT array_agg(c) INTO cells FROM c1; -- 1st guess: centroid - best_c := polygon; + -- best_c := polygon; best_d := cdb_crankshaft._Signed_Dist(polygon, ST_Centroid(Polygon)); -- looping the loop @@ -57,6 +57,7 @@ BEGIN EXIT WHEN i > n; cell := cells[i]; + i := i+1; -- cell side size, it's square @@ -64,13 +65,14 @@ BEGIN -- check distance test_d := cdb_crankshaft._Signed_Dist(polygon, ST_Centroid(cell)); + IF test_d > best_d THEN best_d := test_d; - best_c := cells[i]; + best_c := cell; END IF; -- longest distance within the cell - test_mx := test_d + (test_h/2 * sqr); + test_mx := test_d + (test_h * sqr); -- if the cell has no chance to contains the desired point, continue CONTINUE WHEN test_mx - best_d <= tolerance; @@ -95,6 +97,7 @@ END; $$ language plpgsql IMMUTABLE; + -- signed distance point to polygon with holes -- negative is the point is out the polygon -- rev 1. adding MULTIPOLYGON and GEOMETRYCOLLECTION support by @abelvm From 00327e6de269dc7cda2cdf04a820fcec22dc2908 Mon Sep 17 00:00:00 2001 From: abelvm Date: Tue, 28 Mar 2017 13:38:21 +0200 Subject: [PATCH 6/9] fixed indexing --- src/pg/test/expected/13_pia_test.out | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pg/test/expected/13_pia_test.out b/src/pg/test/expected/13_pia_test.out index 2367e20..56f32f9 100644 --- a/src/pg/test/expected/13_pia_test.out +++ b/src/pg/test/expected/13_pia_test.out @@ -2,6 +2,6 @@ SET client_min_messages TO WARNING; \set ECHO none st_astext ------------------------------------------- - POINT(-3.67484492582767 40.4395084885993) + POINT(-3.67484492582767 40.4394914243877) (1 row) From 7f5edb26b04b8a1b907c90b9a3ed8cd529830a88 Mon Sep 17 00:00:00 2001 From: abelvm Date: Tue, 28 Mar 2017 14:13:49 +0200 Subject: [PATCH 7/9] fixed corner case --- src/pg/sql/13_PIA.sql | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pg/sql/13_PIA.sql b/src/pg/sql/13_PIA.sql index 38963fa..02cafe7 100644 --- a/src/pg/sql/13_PIA.sql +++ b/src/pg/sql/13_PIA.sql @@ -46,7 +46,7 @@ BEGIN SELECT array_agg(c) INTO cells FROM c1; -- 1st guess: centroid - -- best_c := polygon; + best_c := polygon; best_d := cdb_crankshaft._Signed_Dist(polygon, ST_Centroid(Polygon)); -- looping the loop From 807a5373e8002066839e6a8326e1bde9e082c429 Mon Sep 17 00:00:00 2001 From: Andy Eschbacher Date: Wed, 10 Jan 2018 16:35:23 -0500 Subject: [PATCH 8/9] adds simple test --- src/pg/test/expected/13_pia_test.out | 5 +++++ src/pg/test/sql/13_pia_test.sql | 8 ++++++++ 2 files changed, 13 insertions(+) diff --git a/src/pg/test/expected/13_pia_test.out b/src/pg/test/expected/13_pia_test.out index 56f32f9..bdd9c73 100644 --- a/src/pg/test/expected/13_pia_test.out +++ b/src/pg/test/expected/13_pia_test.out @@ -5,3 +5,8 @@ SET client_min_messages TO WARNING; POINT(-3.67484492582767 40.4394914243877) (1 row) + st_astext +------------ + POINT(0 0) +(1 row) + diff --git a/src/pg/test/sql/13_pia_test.sql b/src/pg/test/sql/13_pia_test.sql index 8b37082..7516af3 100644 --- a/src/pg/test/sql/13_pia_test.sql +++ b/src/pg/test/sql/13_pia_test.sql @@ -5,3 +5,11 @@ with a as( select st_geomfromtext('POLYGON((-432540.453078056 4949775.20452642,-432329.947920966 4951361.232584,-431245.028163694 4952223.31516671,-429131.071033529 4951768.00415574,-424622.07505895 4952843.13503987,-423688.327170174 4953499.20752423,-424086.294349759 4954968.38274191,-423068.388925945 4954378.63345336,-423387.653225542 4953355.67417084,-420594.869840519 4953781.00230592,-416026.095299382 4951484.06849063,-412483.018546414 4951024.5410983,-410490.399661215 4954502.24032205,-408186.197521284 4956398.91417441,-407627.262358013 4959300.94633864,-406948.770061627 4959874.85407739,-404949.583326472 4959047.74518163,-402570.908447199 4953743.46829807,-400971.358683991 4952193.11680804,-403533.488084088 4949649.89857885,-406335.177028373 4950193.19571096,-407790.456731515 4952391.46015616,-412060.672398345 4950381.2389307,-410716.93482498 4949156.7509561,-408464.162289794 4943912.8940387,-409350.599394983 4942819.84896006,-408087.791091424 4942451.6711778,-407274.045613725 4940572.4807777,-404446.196589102 4939976.71501489,-402422.964843936 4940450.3670813,-401010.654464241 4939054.8061663,-397647.247369412 4940679.80737878,-395658.413346901 4940528.84765185,-395536.852462953 4938829.79565997,-394268.923462818 4938003.7277717,-393388.720249116 4934757.80596815,-392393.301362444 4934326.71675815,-392573.527618037 4932323.40974412,-393464.640141837 4931903.10653605,-393085.597275686 4931094.7353605,-398426.261165985 4929156.87541607,-398261.174361137 4926238.00816416,-394045.059966834 4925765.18668498,-392982.960705174 4926391.81893628,-393090.272694301 4927176.84692181,-391648.240010564 4924626.06386961,-391889.914625075 4923086.14787613,-394345.177314013 4923235.086036,-395550.878718795 4917812.79243978,-399009.463978251 4912927.7157945,-398948.794855767 4911941.91010796,-398092.636652078 4911806.57392519,-401991.601817112 4911722.9204501,-406225.972607907 4914505.47286319,-411104.994569885 4912569.26941163,-412925.513522316 4913030.3608866,-414630.148884835 4914436.69169949,-414207.691417276 4919205.78028405,-418306.141109809 4917994.9580478,-424184.700779621 4918938.12432889,-426816.961458921 4923664.37379373,-420956.324227126 4923381.98014807,-420186.661267781 4924286.48693378,-420943.411166194 4926812.76394433,-419779.45457046 4928527.43466337,-419768.767899344 4930681.94459216,-421911.668097113 4930432.40620397,-423482.386112205 4933451.28047252,-427272.814773717 4934151.56473242,-427144.908678797 4939731.77191996,-428982.125554848 4940522.84445172,-428986.133056516 4942437.17281266,-431237.792396792 4947309.68284815,-432476.889648814 4947791.74800037,-432540.453078056 4949775.20452642))', 3857) as g ) SELECT st_astext(cdb_crankshaft.CDB_PIA(g)) from a; + +-- square centered on 0,0 with sides of length 2 +-- expectation: point(0, 0) +WITH square AS ( + SELECT 'SRID=4326;POLYGON((-1 1, 1 1, 1 -1, -1 -1, -1 1))'::geometry as g +) +SELECT ST_AsText(cdb_crankshaft.CDB_PIA(g)) + FROM square From 628fd2b839db36fde9ebc1cf93231194343dcfa2 Mon Sep 17 00:00:00 2001 From: Andy Eschbacher Date: Wed, 10 Jan 2018 16:45:36 -0500 Subject: [PATCH 9/9] adds test on multipolygon --- src/pg/test/expected/13_pia_test.out | 5 +++++ src/pg/test/sql/13_pia_test.sql | 11 +++++++++++ 2 files changed, 16 insertions(+) diff --git a/src/pg/test/expected/13_pia_test.out b/src/pg/test/expected/13_pia_test.out index bdd9c73..2ccd544 100644 --- a/src/pg/test/expected/13_pia_test.out +++ b/src/pg/test/expected/13_pia_test.out @@ -10,3 +10,8 @@ SET client_min_messages TO WARNING; POINT(0 0) (1 row) + st_astext +------------ + POINT(0 0) +(1 row) + diff --git a/src/pg/test/sql/13_pia_test.sql b/src/pg/test/sql/13_pia_test.sql index 7516af3..a1c11b9 100644 --- a/src/pg/test/sql/13_pia_test.sql +++ b/src/pg/test/sql/13_pia_test.sql @@ -11,5 +11,16 @@ SELECT st_astext(cdb_crankshaft.CDB_PIA(g)) from a; WITH square AS ( SELECT 'SRID=4326;POLYGON((-1 1, 1 1, 1 -1, -1 -1, -1 1))'::geometry as g ) +SELECT ST_AsText(cdb_crankshaft.CDB_PIA(g)) + FROM square; + +-- MultiPolygon test +-- square centered on 0,0 with sides of length 2 +-- expectation: point(0, 0) +WITH square AS ( + SELECT + ST_Multi('SRID=4326;POLYGON((-1 1, 1 1, 1 -1, -1 -1, -1 1))'::geometry) as g +) SELECT ST_AsText(cdb_crankshaft.CDB_PIA(g)) FROM square +