PostGIS erase logic and speed

Question

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

Answer 1

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
  SELECT 
    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.