Opened 9 years ago

Closed 9 years ago

#3108 closed defect (fixed)

ST_SubDivide leaves pockets

Reported by: robe Owned by: pramsey
Priority: high Milestone: PostGIS 2.2.0
Component: postgis Version: master
Keywords: Cc:

Description

It might be a misunderstanding of mine, but I assume ST_SubDivide should just return my geometries copped into bits.

I'll try to get this down to something a bit more digestible, but I downloaded my favorite test layer from here - http://www.mass.gov/anf/research-and-tech/it-serv-and-support/application-serv/office-of-geographic-information-massgis/datalayers/townsurvey.html

Loaded the one called TOWNSSURVEY_POLYM into a table called

towns:

and then ran this query:

-- after_subdivide
SELECT gid, ST_SubDivide(geom)
FROM towns;

-- before_subdivide
SELECT gid, ST_SubDivide(geom)
FROM towns;

I'll isolate this to something smaller and more digestable, but putting here so I don't forget to complain.

Attachments (10)

before_subdivide.png (140.6 KB ) - added by robe 9 years ago.
Before
after_subdivide.png (167.6 KB ) - added by robe 9 years ago.
after
mb_before_subdivide.png (19.7 KB ) - added by robe 9 years ago.
mb before subdivde
mb_after_subdivide.png (20.1 KB ) - added by robe 9 years ago.
mb after subdivide
subdivide_test.sql (25.7 KB ) - added by robe 9 years ago.
problem town
mb_after_clipbybox2d.png (9.8 KB ) - added by robe 9 years ago.
mb_after_subdivide_450.png (13.9 KB ) - added by robe 9 years ago.
mb_after_subdivide_100.png (12.8 KB ) - added by robe 9 years ago.
mb_before_subdivide_5.png (8.2 KB ) - added by robe 9 years ago.
mb_after_subdivide_5.png (8.8 KB ) - added by robe 9 years ago.

Download all attachments as: .zip

Change History (20)

by robe, 9 years ago

Attachment: before_subdivide.png added

Before

by robe, 9 years ago

Attachment: after_subdivide.png added

after

comment:1 by robe, 9 years ago

Oops forgot running on windows — POSTGIS="2.2.0dev r13480" GEOS="3.5.0dev-CAPI-1.9.0 r4054" SFCGAL="1.0.5" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.11.1, released 2014/09/24" LIBXML="2.7.8" LIBJSON="0.12" RASTER

I'll check the vintage of my GEOS to make sure I'm running latest.

by robe, 9 years ago

Attachment: mb_before_subdivide.png added

mb before subdivde

by robe, 9 years ago

Attachment: mb_after_subdivide.png added

mb after subdivide

by robe, 9 years ago

Attachment: subdivide_test.sql added

problem town

comment:2 by robe, 9 years ago

Okay I took one of the problem towns and simplified to just the point before the problem disappears with this command:

CREATE TABLE subdivide_test AS
SELECT gid, town, ST_Simplify(ST_SnapToGrid(ST_UnaryUnion(geom),5),2)::geometry(POLYGON,26986) As geom
FROM towns
WHERE gid = 93;

Attached is the generated table:

-- mb_after_subdivide

SELECT ST_Subdivide(geom)
FROM subdivide_test;

Give me my land back :)

comment:3 by robe, 9 years ago

and I checked my version of GEOS and it's as fresh as it can be. r4054 was the last commit

comment:4 by robe, 9 years ago

At strk's urging, I tried a similar exercise with ST_ClipByBox and can't get it to fail despite tweaking my cuts.

WITH grid As (
SELECT ST_Tile(
 ST_MakeEmptyRaster(
   ceiling( (ST_XMax(ext) - ST_Xmin(ext))/10 )::integer, 
   ceiling((ST_YMax(ext) - ST_Ymin(ext))/10)::integer, 
   ST_XMin(ext), ST_YMax(ext), 10, -10,0,0, 26986), 700,800
    )::geometry As b
FROM (SELECT ST_Extent(geom) As ext from subdivide_test) As f )
SELECT t.town, ST_Multi(ST_ClipByBox2D(t.geom, grid.b::box2d) )
FROM subdivide_test AS t INNER JOIN grid ON ST_Intersects(t.geom, grid.b) ;

See attached image.

by robe, 9 years ago

Attachment: mb_after_clipbybox2d.png added

comment:5 by robe, 9 years ago

Hmm

-- okay but below 450 I start loosing land -- 
SELECT ST_SubDivide(geom,450)
from subdivide_test;

-- really bad
SELECT ST_SubDivide(geom,100)
from subdivide_test

by robe, 9 years ago

Attachment: mb_after_subdivide_450.png added

by robe, 9 years ago

Attachment: mb_after_subdivide_100.png added

comment:6 by robe, 9 years ago

Priority: mediumhigh

comment:7 by strk, 9 years ago

The code counts the number of vertices that fall within the clipping box. In your case, there are NO vertices in the box, right ? In that case, the code checks if the polygon contains one corner of the box and if it does, it returns a fully filled box, otherwise it returns nothing.

Well, in your cases the box contains NO vertices but it still cuts the polygon, so it's neither a full solid nor a missing. If this analysis is correct, the code is bogus.

It will be a fun exercise at finding the smallest input that can reproduce the problem.

comment:8 by robe, 9 years ago

This test is small enough we can fit in our regress :)

SELECT ST_SubDivide('POLYGON((249045 857495,252985 856270,256785 851705,257130 847260,261225 840210,250170 838400,248445 847565,242075 850775,241470 853295,243165 854590,245930 853740,249045 857495))'::geometry,5);

See attached before and after

by robe, 9 years ago

Attachment: mb_before_subdivide_5.png added

by robe, 9 years ago

Attachment: mb_after_subdivide_5.png added

comment:10 by pramsey, 9 years ago

Resolution: fixed
Status: newclosed

OK, fixed at r13489

Note: See TracTickets for help on using tickets.