Go to content Go to navigation and search

Home

Current PostGIS Blog Articles


Search

Browse

RSS / Atom

Email me

textpattern

Creative Commons License
All Blog Articles, Data Models and Free Source Code by Simon Greener, The SpatialDB Advisor is licensed under a Creative Commons Attribution-ShareAlike 3.0 Unported License.

Gridding a geometry object (PostGIS)

Tuesday December 20 2011 at 22:51

Keywordsgrid gridding raster vector chopping intersection boundary clip polygon line
Summary

A vector-to-raster algorithm generating tiles covering (multi-)line/polygon geometries. Boundary tiles are clipped forming non-square tiles.

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 and REGULARGRIDXY are provided in other blog articles or 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.

  1. WITH geomQuery AS (
  2. SELECT ST_XMIN(g.geom)::NUMERIC AS xmin, round(ST_XMAX(g.geom)::NUMERIC,2)::NUMERIC AS xmax,
  3.        ST_YMIN(g.geom)::NUMERIC AS ymin, round(ST_YMAX(g.geom)::NUMERIC,2)::NUMERIC AS ymax,
  4.        g.geom, 0.050::NUMERIC AS gridX, 0.050::NUMERIC AS gridY, 0::int4 AS loCol, 0::int4 AS loRow
  5.   FROM (SELECT ST_SymDifference(ST_Buffer(a.geom,1.000::NUMERIC),ST_Buffer(a.geom,0.50::NUMERIC)) AS geom
  6.           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
  7.        ) AS g
  8. )
  9. SELECT ROW_NUMBER() OVER (ORDER BY f.gcol, f.grow) AS tid,
  10.        morton((f.gcol-f.loCol),(f.grow-f.loRow)) AS mkey,
  11.        f.gcol,
  12.        f.grow,
  13.        COUNT(*) AS UnionedTileCount,
  14.        ST_Union(f.geom) AS geom
  15.   FROM (SELECT CASE WHEN ST_GeometryType(b.geom) IN ('ST_Polygon','ST_MultiPolygon')
  16.                     THEN ST_Intersection(b.ageom,b.geom)
  17.                     ELSE b.geom
  18.                 END AS geom,
  19.                b.gcol, b.grow, b.loCol, b.loRow
  20.           FROM (SELECT a.geom AS ageom, a.loCol, a.loRow,
  21.                        (RegularGridXYSQL(a.xmin,a.ymin,a.xmax,a.ymax,a.gridX,a.gridY,ST_Srid(a.geom))).*
  22.                   FROM geomQuery AS a
  23.                 ) AS b
  24.          WHERE ST_Intersects(b.ageom,b.geom)
  25.         ) AS f
  26.  WHERE POSITION('POLY' IN UPPER(ST_AsText(f.geom))) > 0
  27.  GROUP BY f.gcol, f.grow, f.loCol, f.loRow
  28.  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.

  1. WITH geomQuery AS (
  2. SELECT g.rid,
  3.        (MIN(ST_XMIN(g.geom))                   OVER (partition BY g.pid))::NUMERIC  AS xmin,
  4.        (MAX(round(ST_XMAX(g.geom)::NUMERIC,2)) OVER (partition BY g.pid))::NUMERIC AS xmax,
  5.        (MIN(ST_YMIN(g.geom))                   OVER (partition BY g.pid))::NUMERIC AS ymin,
  6.        (MAX(round(ST_YMAX(g.geom)::NUMERIC,2)) OVER (partition BY g.pid))::NUMERIC AS ymax,
  7.        g.geom, 0.050::NUMERIC AS gridX, 0.050::NUMERIC AS gridY, 0::int4 AS loCol, 0::int4 AS loRow
  8.   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
  9.           FROM (SELECT 1::int4 AS rid, ST_GeomFromText('POINT(09.50 10.00)',0) AS geom
  10.       UNION ALL SELECT 2::int4 AS rid, ST_GeomFromText('POINT(10.50 10.00)',0) AS geom
  11.       UNION ALL SELECT 3::int4 AS rid, ST_GeomFromText('POINT(10.00 10.50)',0) AS geom
  12.       UNION ALL SELECT 4::int4 AS rid, ST_GeomFromText('POINT(10.00 09.50)',0) AS geom ) a
  13.        ) g                        
  14. )
  15. SELECT ROW_NUMBER() OVER (ORDER BY f.gcol, f.grow) AS tid,
  16.        morton((f.gcol-f.loCol),(f.grow-f.loRow)) AS mkey,
  17.        f.gcol,
  18.        f.grow,
  19.        COUNT(*) AS UnionedTileCount,
  20.        ST_Union(f.geom) AS geom
  21.   FROM (SELECT CASE WHEN ST_GeometryType(b.geom) IN ('ST_Polygon','ST_MultiPolygon')
  22.                     THEN ST_Intersection(b.ageom,b.geom)
  23.                     ELSE b.geom
  24.                 END AS geom,
  25.                b.gcol, b.grow, b.loCol, b.loRow
  26.           FROM (SELECT a.geom AS ageom, a.loCol, a.loRow,
  27.                        (RegularGridXYSQL(a.xmin,a.ymin,a.xmax,a.ymax,a.gridX,a.gridY,ST_Srid(a.geom))).*
  28.                   FROM geomQuery a
  29.                 ) b
  30.          WHERE ST_Intersects(b.ageom,b.geom)
  31.         ) f
  32.  WHERE POSITION('POLY' IN UPPER(ST_AsText(f.geom))) > 0
  33.  GROUP BY f.gcol, f.grow, f.loCol, f.loRow
  34.  ORDER BY  2;

This is what this looks like.

The RegularGridXY and RegularGridSQL functions are:

  1. CREATE OR REPLACE FUNCTION RegularGridXY( p_xmin       NUMERIC,
  2.                                           p_ymin       NUMERIC,
  3.                                           p_xmax       NUMERIC,
  4.                                           p_ymax       NUMERIC,
  5.                                           p_TileSizeX  NUMERIC,
  6.                                           p_TileSizeY  NUMERIC,
  7.                                           p_srid       int4)
  8. RETURNS SETOF T_Grid IMMUTABLE
  9. AS $$
  10. DECLARE
  11.    v_loCol int4;
  12.    v_hiCol int4;
  13.    v_loRow int4;
  14.    v_hiRow int4;
  15.    v_geom  geometry;
  16.    v_grid  t_grid;
  17. BEGIN
  18.    v_loCol := trunc((p_XMIN / p_TileSizeX)::NUMERIC );
  19.    v_hiCol := CEIL( (p_XMAX / p_TileSizeX)::NUMERIC ) - 1;
  20.    v_loRow := trunc((p_YMIN / p_TileSizeY)::NUMERIC );
  21.    v_hiRow := CEIL( (p_YMAX / p_TileSizeY)::NUMERIC ) - 1;
  22.    FOR v_col IN v_loCol..v_hiCol Loop
  23.      FOR v_row IN v_loRow..v_hiRow Loop
  24.          v_geom := ST_SetSRID(ST_MakeBox2D(ST_Point(( v_col * p_TileSizeX),               (v_row * p_TileSizeY)),
  25.                                            ST_Point(((v_col * p_TileSizeX)+p_TileSizeX), ((v_row * p_TileSizeY)+p_TileSizeY))),
  26.                                            p_srid);
  27.          SELECT v_col,v_row,v_geom INTO v_grid;
  28.          -- SELECT v_col,v_row,ST_GeomFromText('POINT(' || v_col || ' ' || v_row ||')',0) INTO v_grid;
  29.          RETURN NEXT v_grid;
  30.      END Loop;
  31.    END Loop;
  32. END;
  33. $$ LANGUAGE 'plpgsql';
  34. .
  35. CREATE OR REPLACE FUNCTION RegularGridXYSQL(p_xmin       NUMERIC,
  36.                                             p_ymin       NUMERIC,
  37.                                             p_xmax       NUMERIC,
  38.                                             p_ymax       NUMERIC,
  39.                                             p_TileSizeX  NUMERIC,
  40.                                             p_TileSizeY  NUMERIC,
  41.                                             p_srid       int4)
  42.   RETURNS SETOF T_Grid IMMUTABLE
  43. AS
  44. $$
  45.   SELECT * FROM RegularGridXY($1,$2,$3,$4,$5,$6,$7);
  46. $$
  47. LANGUAGE 'sql';
I hope this is useful to a PostGIS user.

Creative Commons License

post this at del.icio.uspost this at Diggpost this at Technoratipost this at Redditpost this at Farkpost this at Yahoo! my webpost this at Windows Livepost this at Google Bookmarkspost this to Twitter

Comment

Article Navigation:   Previous