I have 2 Polygon layers. 3.5 million records of all NJ Parcels and a 1 record state wide polygon shapefile of NJ. My goal is to emulate the ArcGIS erase analysis too where I would essentially erase all the areas in the NJ shapefile where the 3.5 million parcels intersect and what I would be left with is a statewide polygon of all NJ roadways.

Here is the query I have put together, it seems fairly straightforward but my question is threefold. 1. should I be using st_intersects instead of && 2. is my query constructed correctly so it will achieve my results? 3. is there a faster way for me to construct this query? indexes, etc.. because it has been running for about an hour already.

create table streetbyerase as
select st_difference(parcels.geom,nj.geom) as geom
from parcels,nj 
where parcels.geom && nj.geom
  • maybe union all the parcels first into one big geometry? – DPSSpatial Nov 10 '16 at 18:19
  • do you think that would make a difference in the speed of the query or better for the desired results? – ziggy Nov 10 '16 at 18:20
  • b/c i would think unionizing 3.5 million parcels would probably take even longer that the above query..? – ziggy Nov 10 '16 at 18:21
  • I don't have a dataset that big to test ... – DPSSpatial Nov 10 '16 at 19:21
  • Do you mean remove the roads that are outside the Joisey state boundary and then combine the result into one multilinestring? – John Powell Nov 10 '16 at 19:24

In the spirit of never saying "die", and given the negative attitude in my comment, here's an approach that might work, and does use indexes.

-- Turn NJ into a large number of small tractable areas
CREATE SEQUENCE nj_square_id;
CREATE TABLE nj_squares AS
    nextval('nj_square_id') AS nj_id, 
    ST_SubDivide(geom) AS geom
  FROM nj;

-- Index the squares for faster searching
CREATE INDEX nj_squares_x ON nj_squares USING GIST (geom);

-- Index parcels too in case you forgot
CREATE INDEX parcels_x ON parcels USING GIST (geom);

-- For each square, compute "bits that aren't parcels"
CREATE TABLE nj_not_parcels AS
WITH parcel_polys AS (
  SELECT nj.nj_id, ST_Union(p.geom) AS geom
  FROM nj_squares nj
  JOIN parcels p
  ON ST_Intersects(p.geom, nj.geom)
  GROUP BY nj.nj_id
SELECT nj_id,
  ST_Difference(nj.geom, pp.geom) AS geom
FROM parcel_polys pp 
JOIN nj_squares
USING (nj_id);

That should work passably well. If it's still too slow, using some gridding technique other than sub-divide to create more smaller squares as step one would be the trick.

  • I will try out this answer when I get back to work on monday – ziggy Nov 11 '16 at 17:54
  • your query has been running for 30 hours so far lol, is this expected given the large dateset I am working with? – ziggy Nov 15 '16 at 19:38
  • 1
    Well, maybe so, because the ST_Subdivide() may create some objectively very large polygons in the middle of the state that then result in very large unions of parcels. Regardless, you do have a computationally intensive question, so I'd expect it to take a while, I imagine there's quite a few parcels in NJ. You can multi-step it too, by creating the parcel_polys table as a table instead of a CTE. You can also run some smaller test queries, doing one nj_square at a time. Experiment, see how bad/good parts of it are. In the extreme, it's possible to parallelize it on nj_squares. – Paul Ramsey Nov 15 '16 at 20:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.