Home

Current PostGIS Blog Articles

Search

Browse

Email me

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 11: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 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.

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 required T_Grid type and RegularGridXY and RegularGridSQL functions are:

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

### Comment [2]

Hi,

The Morton function dosen’t seem to be included in this post. Could you point me in the direction of the correct function I did have a look around but haven’t managed to track it down.

Thanks,

Mark

— Mark Davidson · 2 October 2012, 12:23 · #

I have modified the article to link to the correct page.
Thanks for pointing out the linkage problem.
regards
Simon

— Simon Greener · 3 October 2012, 02:00 · #

 Name Remember E-mail Message Textile Help