[postgis-users] A bit of advice would be greatly appreciated

[James G Wilkinson]

All,

I have been monitoring this forum for some time, and it has been
very helpful in my efforts.  I have run into a problem that I simply
cannot seem to get a handle on.  I have two spatial data sets:
regular_grid (gid integer, area double, icell integer, jcell integer,
the_geom geometry) and plant_area (gid integer, area double,
area_fraction double, the_geom geometry).  regular_grid represents
a 36 km by 36 km set of grid cells (polygons) that covers the
continental United States and is about 10 MB in size.  plant_area is
a polygon coverage that covers North America which contains
the area covered by a single plant species and is about 6 GB in
size (each polygon in plant_area is roughly on the order of
4,000,000 square meters).  Both coverages are in LAT/LON
coordinates.  I am trying to perform the following:

CREATE TABLE intersect_results AS
   SELECT a.icell, a.jcell, SUM(AREA(INTERSECTION(a.the_geom,b.the_geom)) /
                            b.area * area_fraction *
                            AREA(TRANSFORM(b.the_geom,32768)))) as area
   FROM regular_grid a, plant_area b
   WHERE b.the_geom && a.the_geom AND
                  DISTANCE(a.the_geom,b.the_geom) = 0
   GROUP BY a.icell, a.jcell
   ORDER BY a.icell, a.jcell;

The TRANSFORM to 32768 gives me the area in square meters
in a Lambert Conformal Conic projection that I routinely use in my
work.  Both tables are indexed using GIST on the geometries.  I have
run VACUUM ANALYZE.  I aborted the query on my Opteron
machine after 24 hours of CPU time (it was the only thing running).

I have not spent anytime trying to diagnose the query plan results
as I am almost completely ignorant of how to use them -- and maybe
this is the time to fix that gap in my knowledge.  Regardless, I am holding
out some hope that some expert on this list can provide some advice on
how better to render this query.

Regards,
terrakit