ST_SubDivide leaves pockets

[robe]

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.

[robe]

Before

[robe]

after

[robe]

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.

[robe]

mb before subdivde

[robe]

mb after subdivide

[robe]

problem town

[robe]

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

[robe]

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

[robe]

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.

[robe]

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

[strk]

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.

[robe]

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

[pramsey]

OK, fixed at r13489