Monday, 9 July 2018

postgis - How do I work around floating-point errors using ST_SubDivide?


I'm using ST_SubDivide() to speed up some complex calculations to answer a question. Yet I'm noticing there is a different boundary on the geometry that results from ST_SubDivide(). Essentially, I have a query that works reasonably well until ST_SubDivide(). Then it intersects with other polygons that it wouldn't otherwise intersect with. I'm running an query and joining with ST_Intersects in my actual app.


# SELECT gid, ST_Area(ST_Union(geom)) - ST_Area(state.geog::geometry) 
FROM census.state
JOIN (
SELECT
gid,

ST_SubDivide(geog::geometry,10) AS geom
FROM census.state
) AS t
USING (gid)
GROUP BY gid
HAVING NOT ST_Equals(state.geog::geometry, ST_Union(geom));

gid | ?column?
-----+-----------------------
34 | 1.03028696685215e-13

43 | -1.18529221992958e-08
25 | -3.64934146901419e-07
32 | 8.88178419700125e-16
8 | -8.51020365288946e-09
12 | 5.55111512312578e-16
1 | -7.74611432774464e-08
10 | -1.11022302462516e-14
26 | -5.37694003810429e-06
... contd.


For reference, I'm trying this with 2016 TIGER state date. I loaded it with


shp2pgsql -DIGc -s 4326 census/tiger_state/*.shp census.state | psql`

I got the idea from a video cast by Regina Obe



For any geometry you give it, it will chop it up into smaller bits. And this is useful for performance reasons. For example, if you have a huge state it's often useful to chop that into smaller polygons. So when you do intersects the number of points that it has to be dealt with are a lot fewer. So here is an example where I take four states, so I have four rows and then I pass it through the SubDivide function. [...] So you have four states converted to 186 polygons.





No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...