[postgis-users] ST_Intersects of a diagonal buffer on a polygon grid

[Simon Greener]

Michael,

I had a quick look at this and decided to tackle the problem slightly differently. What I did was take your diagonal line and break it up into "vectors" (ie two point line segments) which I then buffered and segmentized and inserted into your table.  The changed section of your script is:

insert into throwaway_polybuff(the_geom)
select ST_Segmentize(ST_Buffer(ST_MakeLine(ST_MakePoint((startcoord).x,(startcoord).y), 
                               ST_MakePoint((endcoord).x,(endcoord).y)), 1000), 5000) as the_geom
  from GetVectorViaSQL('LINESTRING(200 403000,186000 367000,254000 134000,424000 23000,424000 23000)'::geometry) as p;

SELECT g.*
-- The Query --
-- Returns 2076 polygons but only 2095 with modified polybuff
-- Takes about 60 seconds on my computer (and 45 on mine)
-- but with the modified data it takes where now it take 3.3seconds
  FROM throwaway_grid g, 
       throwaway_polybuff p
 WHERE ST_Intersects(g.the_geom,p.the_geom);

SELECT g.*
-- Returns 2095 polygons 
-- With this one taking 17 seconds
  FROM throwaway_grid g, 
       (select ST_Segmentize(ST_Buffer(ST_MakeLine(ST_MakePoint((startcoord).x,(startcoord).y), 
                                       ST_MakePoint((endcoord).x,(endcoord).y)), 1000), 5000) as the_geom
          from GetVectorViaSQL('LINESTRING(200 403000,186000 367000,254000 134000,424000 23000,424000 23000)'::geometry) ) as p
 WHERE ST_Intersects(g.the_geom,p.the_geom );

The GetVector PgPLSQL function is:

CREATE OR REPLACE FUNCTION GetVectorViaSQL(p_geometry geometry)
  RETURNS SETOF VectorType IMMUTABLE  
AS $$
DECLARE
    v_GeometryType   varchar(1000);
    v_rec            RECORD;
    v_vector         VectorType;
    v_start          CoordType;
    v_end            CoordType;
    c_linestring CURSOR ( p_geom geometry ) 
    IS
      SELECT X(sp) as sx,Y(sp) as sy,Z(sp) as sz,M(sp) as sm,X(ep) as ex,Y(ep) as ey,Z(ep) as ez,M(ep) as em
       FROM ( SELECT pointn(p_geom, generate_series(1, npoints(p_geom)-1)) as sp,
                     pointn(p_geom, generate_series(2, npoints(p_geom)  )) as ep
            ) AS foo;
    c_multilinestring CURSOR ( p_geom geometry ) 
    IS
      SELECT X(sp),Y(sp),Z(sp),M(sp),X(ep),Y(ep),Z(ep),M(ep)
        FROM ( SELECT pointn(b.geom, generate_series(1, npoints(b.geom)-1)) as sp,
                      pointn(b.geom, generate_series(2, npoints(b.geom)  )) as ep
                 FROM (select geometryn(p_geom,generate_series(1,numgeometries(p_geom))) as geom) as b
             ) AS foo;
    c_polygon CURSOR ( p_geom geometry )
    IS
      SELECT X(sp),Y(sp),Z(sp),M(sp),X(ep),Y(ep),Z(ep),M(ep)
        FROM ( SELECT pointn(a.geom, generate_series(1, npoints(a.geom)-1)) as sp,
                      pointn(a.geom, generate_series(2, npoints(a.geom)  )) as ep
                 FROM ( SELECT exteriorring(p_geom) as geom
                        UNION ALL
                        SELECT interiorringn(p_geom,generate_series(1,numinteriorrings(p_geom))) as geom
                      ) a
             ) b;
    c_multipolygon CURSOR ( p_geom geometry )
    IS
      SELECT X(sp),Y(sp),Z(sp),M(sp),X(ep),Y(ep),Z(ep),M(ep)
        FROM ( SELECT pointn(a.geom, generate_series(1, npoints(a.geom)-1)) as sp,
                      pointn(a.geom, generate_series(2, npoints(a.geom)  )) as ep
                 FROM ( SELECT exteriorring(b.geom) as geom
                          FROM (select geometryn(p_geom,generate_series(1,numgeometries(p_geom))) as geom ) b    
                        UNION ALL
                        SELECT interiorringn(c.geom,generate_series(1,numinteriorrings(c.geom))) as geom
                          FROM (select geometryn(p_geom,generate_series(1,numgeometries(p_geom))) as geom ) c   
                      ) a
                    ) b;
    c_refcursor refcursor;
Begin
    If ( p_geometry is NULL ) Then
      return;
    End If;
    v_GeometryType := ST_GeometryType(p_geometry);
    If    ( v_GeometryType in ('ST_Point','ST_MultiPoint') ) Then
      return;
    End If;
    IF    ( v_GeometryType = 'ST_MultiLineString' ) THEN
      OPEN c_multilinestring(p_geometry);
      c_refcursor := c_multilinestring;
    ELSIF ( v_GeometryType = 'ST_LineString'      ) THEN
      OPEN c_linestring(p_geometry);
      c_refcursor := c_linestring;
    ELSIF ( v_GeometryType = 'ST_MultiPolygon'    ) THEN
      OPEN c_multipolygon(p_geometry);
      c_refcursor := c_multipolygon;
    ELSIF ( v_GeometryType = 'ST_Polygon'         ) THEN
      OPEN c_polygon(p_geometry);
      c_refcursor := c_polygon;
    ELSIF ( v_GeometryType = 'ST_Geometry'        ) THEN
      -- COuld be anything... use native PostGIS function
      v_GeometryType = GeometryType(p_geometry);
      IF ( v_geometryType = 'GEOMETRYCOLLECTION' ) THEN
         FOR v_geom IN 1..ST_NumGeometries(p_geometry) LOOP
            FOR v_rec IN SELECT * FROM GetVectorViaSQL(ST_GeometryN(p_geometry,v_geom)) LOOP
   	        RETURN NEXT v_rec;
   	    END LOOP;
         END LOOP;
      ELSE
         -- Probably CURVED something...
         RETURN;
      END IF;
    END IF;
    IF ( v_geometryType NOT IN ('ST_Geometry','GEOMETRYCOLLECTION') ) THEN
      LOOP
        FETCH c_refcursor INTO 
              v_start.x, v_start.y, v_start.z, v_start.m,
              v_end.x,   v_end.y,   v_end.z,   v_end.m;
        v_vector.startcoord := v_start;
        v_vector.endcoord   := v_end;
        EXIT WHEN NOT FOUND;
        RETURN NEXT v_vector;
      END LOOP;
      CLOSE c_refcursor;
    END IF;
end;
$$ LANGUAGE 'plpgsql';

regards
Simon
-- 
SpatialDB Advice and Design, Solutions Architecture and Programming,
Oracle Database 10g Administrator Certified Associate; Oracle Database 10g SQL Certified Professional
Oracle Spatial, SQL Server, PostGIS, MySQL, ArcSDE, Manifold GIS, Radius Topology and Studio Specialist.
39 Cliff View Drive, Allens Rivulet, 7150, Tasmania, Australia.
Website: www.spatialdbadvisor.com
  Email: simon at spatialdbadvisor.com
  Voice: +613 9016 3910
Mobile: +61 418 396391
Skype: sggreener
Longitude: 147.20515 (147° 12' 18" E)
Latitude: -43.01530 (43° 00' 55" S)
NAC:W80CK 7SWP3