Top 5 Recent Articles

ARTICLES CATEGORIES

Gridding a geometry object (PostGIS)

A common question that comes up in the various database forums (PostGIS, SQL Server, Oracle) is how to “grid” a linear or polygonal object.

By “grid” one means work out the square pixels (rectangular polygons) that cover or define a vector geometry.

Here is some SQL for doing this for PostGIS.

Single Geometry

The following grids a single geometry into a single logical grid. The function Morton is provided in another article, while RegularGridXY is provided at the end of this article. Note that RegularGridXY is used in order to ensure that the gridding takes place according to a commonly defined MBR.

 WITH geomQuery AS (
 SELECT ST_XMIN(g.geom)::NUMERIC AS xmin, round(ST_XMAX(g.geom)::NUMERIC,2)::NUMERIC AS xmax,
        ST_YMIN(g.geom)::NUMERIC AS ymin, round(ST_YMAX(g.geom)::NUMERIC,2)::NUMERIC AS ymax,
        g.geom, 0.050::NUMERIC AS gridX, 0.050::NUMERIC AS gridY, 0::int4 AS loCol, 0::int4 AS loRow
   FROM (SELECT ST_SymDifference(ST_Buffer(a.geom,1.000::NUMERIC),ST_Buffer(a.geom,0.50::NUMERIC)) AS geom
           FROM (SELECT ST_GeomFromText('MULTIPOINT((09.25 10.00),(10.75 10.00),(10.00 10.75),(10.00 9.25))',0) AS geom ) AS a
        ) AS g
 )
 SELECT ROW_NUMBER() OVER (ORDER BY f.gcol, f.grow) AS tid,
        morton((f.gcol-f.loCol),(f.grow-f.loRow)) AS mkey,
        f.gcol,
        f.grow,
        COUNT(*) AS UnionedTileCount,
        ST_Union(f.geom) AS geom
   FROM (SELECT CASE WHEN ST_GeometryType(b.geom) IN ('ST_Polygon','ST_MultiPolygon')
                     THEN ST_Intersection(b.ageom,b.geom)
                     ELSE b.geom
                 END AS geom,
                b.gcol, b.grow, b.loCol, b.loRow
           FROM (SELECT a.geom AS ageom, a.loCol, a.loRow,
                        (RegularGridXYSQL(a.xmin,a.ymin,a.xmax,a.ymax,a.gridX,a.gridY,ST_Srid(a.geom))).*
                   FROM geomQuery AS a
                 ) AS b
          WHERE ST_Intersects(b.ageom,b.geom)
         ) AS f
  WHERE POSITION('POLY' IN UPPER(ST_AsText(f.geom))) > 0
  GROUP BY f.gcol, f.grow, f.loCol, f.loRow
  ORDER BY  2;

This is what this looks like.

Multiple Geometries

The following SQL grids multiple geometries into a single logical grid.

This aggregate ST_XMin etc SQL and RegularGridXY function ensures that the grids for one object align with the grids for the others.

 WITH geomQuery AS (
 SELECT g.rid,
        (MIN(ST_XMIN(g.geom))                   OVER (partition BY g.pid))::NUMERIC  AS xmin,
        (MAX(round(ST_XMAX(g.geom)::NUMERIC,2)) OVER (partition BY g.pid))::NUMERIC AS xmax,
        (MIN(ST_YMIN(g.geom))                   OVER (partition BY g.pid))::NUMERIC AS ymin,
        (MAX(round(ST_YMAX(g.geom)::NUMERIC,2)) OVER (partition BY g.pid))::NUMERIC AS ymax,
        g.geom, 0.050::NUMERIC AS gridX, 0.050::NUMERIC AS gridY, 0::int4 AS loCol, 0::int4 AS loRow
   FROM (SELECT 1::int4 AS pid, a.rid, ST_SymDifference(ST_Buffer(a.geom,1.000::NUMERIC),ST_Buffer(a.geom,0.750::NUMERIC)) AS geom
           FROM (SELECT 1::int4 AS rid, ST_GeomFromText('POINT(09.50 10.00)',0) AS geom
       UNION ALL SELECT 2::int4 AS rid, ST_GeomFromText('POINT(10.50 10.00)',0) AS geom
       UNION ALL SELECT 3::int4 AS rid, ST_GeomFromText('POINT(10.00 10.50)',0) AS geom
       UNION ALL SELECT 4::int4 AS rid, ST_GeomFromText('POINT(10.00 09.50)',0) AS geom ) a
        ) g
 )
 SELECT ROW_NUMBER() OVER (ORDER BY f.gcol, f.grow) AS tid,
        morton((f.gcol-f.loCol),(f.grow-f.loRow)) AS mkey,
        f.gcol,
        f.grow,
        COUNT(*) AS UnionedTileCount,
        ST_Union(f.geom) AS geom
   FROM (SELECT CASE WHEN ST_GeometryType(b.geom) IN ('ST_Polygon','ST_MultiPolygon')
                     THEN ST_Intersection(b.ageom,b.geom)
                     ELSE b.geom
                 END AS geom,
                b.gcol, b.grow, b.loCol, b.loRow
           FROM (SELECT a.geom AS ageom, a.loCol, a.loRow,
                        (RegularGridXYSQL(a.xmin,a.ymin,a.xmax,a.ymax,a.gridX,a.gridY,ST_Srid(a.geom))).*
                   FROM geomQuery a
                 ) b
         WHERE ST_Intersects(b.ageom,b.geom)
         ) f
  WHERE POSITION('POLY' IN UPPER(ST_AsText(f.geom))) > 0 
  GROUP BY f.gcol, f.grow, f.loCol, f.loRow
  ORDER BY  2;

This is what this looks like.

The required T_Grid type and RegularGridXY and RegularGridSQL functions are:

 DROP   TYPE T_Grid cascade;
 CREATE TYPE T_Grid AS
 (gcol  int4,
  grow  int4,
  geom  geometry);
 .
 CREATE OR REPLACE FUNCTION RegularGridXY( p_xmin       NUMERIC,
                                             p_ymin       NUMERIC,
                                             p_xmax       NUMERIC,
                                             p_ymax       NUMERIC,
                                             p_TileSizeX  NUMERIC,
                                             p_TileSizeY  NUMERIC,
                                             p_srid       int4)
 RETURNS SETOF T_Grid IMMUTABLE
 AS $$
 DECLARE
    v_loCol int4;
    v_hiCol int4;
    v_loRow int4;
    v_hiRow int4;
    v_geom  geometry;
    v_grid  t_grid;
 BEGIN
    v_loCol := trunc((p_XMIN / p_TileSizeX)::NUMERIC );
    v_hiCol := CEIL( (p_XMAX / p_TileSizeX)::NUMERIC ) - 1;
    v_loRow := trunc((p_YMIN / p_TileSizeY)::NUMERIC );
    v_hiRow := CEIL( (p_YMAX / p_TileSizeY)::NUMERIC ) - 1;
    FOR v_col IN v_loCol..v_hiCol Loop
      FOR v_row IN v_loRow..v_hiRow Loop
          v_geom := ST_SetSRID(ST_MakeBox2D(ST_Point(( v_col * p_TileSizeX),               (v_row * p_TileSizeY)),
                                            ST_Point(((v_col * p_TileSizeX)+p_TileSizeX), ((v_row * p_TileSizeY)+p_TileSizeY))),
                                            p_srid);
          SELECT v_col,v_row,v_geom INTO v_grid;
          -- SELECT v_col,v_row,ST_GeomFromText('POINT(' || v_col || ' ' || v_row ||')',0) INTO v_grid;
          RETURN NEXT v_grid;
      END Loop;
    END Loop;
 END;
 $$ LANGUAGE 'plpgsql';
 .
 CREATE OR REPLACE FUNCTION RegularGridXYSQL(p_xmin       NUMERIC,
                                               p_ymin       NUMERIC,
                                               p_xmax       NUMERIC,
                                               p_ymax       NUMERIC,
                                               p_TileSizeX  NUMERIC,
                                               p_TileSizeY  NUMERIC,
                                               p_srid       int4)
   RETURNS SETOF T_Grid IMMUTABLE
 AS
 $$
   SELECT * FROM RegularGridXY($1,$2,$3,$4,$5,$6,$7);
 $$
 LANGUAGE 'sql';

I hope this is useful to a PostGIS user.