[postgis-users] How to create a spatial index for a grid in a subquery?

[Tyler Erickson]

I am trying to create a query that intersects polygons with a regular grid. 
I would like to create the grid on-the-fly, so that any projection, extent,
and spacing could be used.

The query I've come up with is based off of example code found in the
PostGIS in Action book (Listing 5.3), which uses generate_series() and
ST_MakeBox2d() to create a grid in a subquery.  While this works for fine
for coarse grids, this approach quickly slows down for fine grids,
presumably because the grid created in the subquery does not have spatial
index.

What would be a good approach for creating and intersecting with a fine
grid?  Should I create the grid in a temp table, create a spatial index,
then intersect it with polygons?  Or is there a way to create a spatial
index within the a subquery?

Below are examples of my current query, and the example query from PostGIS
in Action.

- Tyler


-- Intersect polygon with grid (151 x 301 = 45451 cells; 14s execution time)
SELECT
  CAST(x AS text) || ' ' || CAST(y As text) As grid_xy,
  ST_AsText(ST_Intersection(g1.geom1, g2.geom2)) As intersect_geom
FROM
    (
    SELECT
      ST_GeomFromText('POLYGON((2 4.5,3 2.6,3 1.8,2 0,-1.5 2.2,0.056
3.222,-1.5 4.2,2 6.5,2 4.5))') As geom1
    ) As g1
INNER JOIN
    (
    SELECT 
      g0.x, g0.y,
      CAST(ST_MakeBox2d(ST_Point(-1.5 + g0.x, 0 + g0.y), ST_Point(-1.5 +
g0.x + 2, 0 + g0.y + 2)) As geometry) As geom2
    FROM
      (
      SELECT 
          x_raw/100.0 as x,
          y_raw/100.0 as y
      FROM generate_series(0,300,2) As x_raw 
      CROSS JOIN generate_series(0,600,2) As y_raw
      ) as g0
    ) As g2
ON ST_Intersects(g1.geom1,g2.geom2);


-- Intersect polygon with grid (2 x 4 = 8 cells;  12ms execution time)
-- Source: PostGIS in Action, Listing 5.2(3)
SELECT
  CAST(x AS text) || ' ' || CAST(y As text) As grid_xy,
  ST_AsText(ST_Intersection(g1.geom1, g2.geom2)) As intersect_geom
FROM
  (SELECT 
     ST_GeomFromText('POLYGON((2 4.5,3 2.6,3 1.8,2 0,-1.5 2.2,0.056
3.222,-1.5 4.2,2 6.5,2 4.5))') As geom1
   ) As g1
   INNER JOIN
   (SELECT
      x, y, CAST(ST_MakeBox2d(ST_Point(-1.5 + x, 0 + y),
      ST_Point(-1.5 + x + 2, 0 + y + 2)) As geometry) As geom2
    FROM generate_series(0,3,2) As x CROSS JOIN generate_series(0,6,2) As y
    ) As g2
ON ST_Intersects(g1.geom1,g2.geom2);
-- 
View this message in context: http://old.nabble.com/How-to-create-a-spatial-index-for-a-grid-in-a-subquery--tp28248823p28248823.html
Sent from the PostGIS - User mailing list archive at Nabble.com.