Creating a Grid of Geographic Tiles

Over the years I have made available functions that generate the tiles of a grid given user input.

These have been made available in Oracle, PostGIS and SQL Server Spatial (geometry).

Here are the type signatures of these functions:

Oracle Object

CREATE TYPE T_Grid AS OBJECT (
  gcol  number, -- Column Reference 
  grow  number, -- Row Reference
  geom  mdsys.sdo_geometry -- SDO_GEOMETRY coded as Optimized Rectangle.
);

CREATE OBJECT T_GEOMETRY....

   Member Function ST_Tile(p_Tile_X    In number,
                          p_Tile_Y    In Number,
                          p_grid_type in varchar2 Default 'TILE',
                          p_option    in varchar2 default 'TOUCH',
                          p_unit      in varchar2 default null)
            Return T_Grids Pipelined,

Oracle Package

CREATE TYPE t_grid AS OBJECT (
   gcol  number,
   grow  number,
   geom  mdsys.sdo_geometry
);

CREATE PACKAGE TESSELATE .....

   TYPE t_GridSet   IS TABLE OF &&defaultSchema..T_Grid;

   Function RegularGridXY(p_xmin      in number,
                          p_ymin      in number,
                          p_xmax      in number,
                          p_ymax      in number,
                          p_TileSizeX in number,
                          p_TileSizeY in number,
                          p_srid      in pls_integer,
                          p_point     in pls_integer Default 0)
     return TESSELATE.t_GridSet pipelined;

   Function RegularGrid( p_LL       In MDSYS.SDO_Point_Type,
                         p_UR       In MDSYS.SDO_Point_Type,
                         p_TileSize In MDSYS.SDO_Point_Type,
                         p_srid     In Number Default Null,
                         p_point    in pls_integer default 0)
            Return TESSELATE.t_GridSet Pipelined;

   Function RegularGridPoint(p_LL       In MDSYS.SDO_POINT_TYPE,
                             p_UR       In MDSYS.SDO_POINT_TYPE,
                             p_TileSize In MDSYS.SDO_POINT_TYPE,
                             p_srid     In Number Default Null )
            Return TESSELATE.t_GridSet Pipelined;

   -- Version that uses a 2003/3003 Optimized Rectangle geometry defining the required extent (eg from SDO_AGGR_UNION)
   -- p_point = 1 means add SDO_POINT to SDO_GEOMETRY; = 0 means don't!
   Function RegularGrid( p_geometry In MDSYS.SDO_Geometry,
                         p_TileSize In MDSYS.SDO_Point_Type,
                         p_point    in pls_integer Default 0)
            Return TESSELATE.t_GridSet Pipelined;

   Function RegularGrid( p_geometry In MDSYS.SDO_Geometry,
                         p_divisor  In Number,
                         p_point    in pls_integer := 0)
            Return TESSELATE.t_GridSet Pipelined;

PostGIS

CREATE TYPE spdba.T_Grid AS (
 gcol  int4,
 grow  int4,
 geom  geometry
);

CREATE FUNCTION spdba.ST_GridGeometry(
  p_geometry    geometry,
  p_TileSizeX   numeric,
  p_TileSizeY   numeric,
  p_rotateX     numeric,
  p_rotateY     numeric,
  p_rotateAngle numeric,
  p_point       boolean -- Return point rather than polygon
)
Returns SETOF spdba.T_Grid IMMUTABLE

CREATE FUNCTION spdba.ST_GridFromXY(
  p_xmin        numeric,
  p_ymin        numeric,
  p_xmax        numeric,
  p_ymax        numeric,
  p_TileSizeX   numeric,
  p_TileSizeY   numeric,
  p_rotateX     numeric,
  p_rotateY     numeric,
  p_rotateAngle numeric,
  p_srid        int4,
  p_point       boolean -- Return point rather than polygon
)
Returns SETOF spdba.T_Grid IMMUTABLE

CREATE FUNCTION spdba.ST_GridFromPoint(
  p_geometry    geometry,
  p_TileSizeX   numeric,
  p_TileSizeY   numeric,
  p_NumTilesX   integer,
  p_NumTilesY   integer,
  p_rotateX     numeric,
  p_rotateY     numeric,
  p_rotateAngle numeric,
  p_point       boolean -- Return point rather than polygon
)
Returns SETOF spdba.T_Grid IMMUTABLE

SQL Server Geometry

CREATE FUNCTION [spdba].[STTileByNumGrids] (
  @p_geometry  geometry,
  @p_NumGridsX integer,
  @p_NumGridsY integer, 
  @p_rPoint    geometry, -- Point(rx,ry)
  @p_rAngle    float
) Returns @table table (
  col  Int,
  row  Int,
  geom geometry 
)

CREATE FUNCTION [spdba].[STTileGeom] (
  @p_geometry geometry,
  @p_TileX    float,
  @p_TileY    float,
  @p_rx       float,
  @p_ry       float,
  @p_rangle   float 
) returns @table table Int,
  row  Int,
  geom geometry
)

CREATE FUNCTION [spdba].[STTileGeomByPoint] (
  @p_point   geometry,
  @p_numTileX integer,
  @p_numTileY integer,
  @p_TileX      float,
  @p_TileY      float,
  @p_rAngle     float
) returns @table table (
  col  Int,
  row  Int,
  geom geometry
)

Requirement

A customer has the requirement to generation grids over archaeological and construction sites which are small in nature.

The work crosses country boundaries and national grids so that the solution has to be in in the geographic space as a common denominator.

The only way they present the grids is via Google Maps using Leaflet (EPSG 4326).

SQL Server Geography

All of the above tiling functions operate in projected space. Generating a grid in Oracle Spatial or PostGIS can be done by creating the tiles in projected space then transforming the result into geographic space via SDO_CS.TRANFORM (Oracle) or ST_Transform (PostGIS).

This option is not available in SQL Server Spatial due to the lack of a native/built-in projection engine.

Thus generation must be done in geographic/geodetic space.

This is now possible via a new function which I have added to the SQL Server functions (which is available via a donation on this website).

This function generates tiles in geographic (ellipsoidal) space.

Here is its function declaration:

CREATE FUNCTION [spdba].[STTileGeogByPoint] (
  @p_point  geography,
  @p_numTileX integer,
  @p_numTileY integer,
  @p_TileX      float,
  @p_TileY      float,
  @p_rAngle     float
) returns @table table (
  col  Int,
  row  Int,
  geom geography
)

The function depends on another function called [cogo].[STDirectVincenty] which, using an ellipsoidal earth model, computes a destination point given a start point, and initial bearing, and a distance. The calculations use the direct solution of geodesics on the ellipsoid devised by Thaddeus Vincenty.

CREATE FUNCTION [cogo].[STDirectVincenty] (
  @p_point    geography,
  @p_startBearing float,
  @p_distance     float
) 
Returns geography

STDirectVincenty is SRID aware. It uses @p_point.STSrid) to query and extract the SRID’s ellipsoid parameters from sys.spatial_reference_systems.

Example

select col,row, geom.STAsText() as tileGeom
  from [dbo].[STTileGeogByPoint] ( 
               geography::Point(55.634269978244,12.051864414446,4326),
               2,2,
               10.0, 15.0,
               22.5
          ) as t;
GO

col row tileGeom
0   0   POLYGON ((12.052084452911 55.634218419749, 12.052304491086 55.634166861254, 12.052365253079 55.634249843067, 12.052145214983 55.634301401561, 12.052084452911 55.634218419749))
0   1   POLYGON ((12.052145214983 55.634301401561, 12.052365253624 55.634249843067, 12.052426015745 55.634332824878, 12.052205977184 55.634384383372, 12.052145214983 55.634301401561))
1   0   POLYGON ((12.052304490797 55.63416686086, 12.052524528684 55.634115302364, 12.052585290597 55.634198284177, 12.05236525279 55.634249842672, 12.052304490797 55.63416686086))
1   1   POLYGON ((12.05236525279 55.634249842672, 12.052585291142 55.634198284177, 12.052646053184 55.63428126599, 12.052426014912 55.634332824484, 12.05236525279 55.634249842672))

Increasing the number of tiles from 2,2 to 10,10 delivers the following (laid over OpenStreetMap)

The following SQL demonstrates that the generated geographic tiles sides are indeed correct even when rotated.

select ROUND(v.length,3) as Side_length, v.geom.STAsText() as tileGeom
  from [dbo].[STTileGeogByPoint] ( 
         geography::Point(55.634269978244,12.051864414446,4326),
         1,1, 10.0,15.0, 22.5) as t
       cross apply
	   [dbo].[STVectorize] (
	     [dbo].[STToGeometry] (t.geom,4326)
       ) as v;
 GO

SideLength tileGeom
15         LINESTRING (12.052084452911 55.634218419749, 12.052304491086 55.634166861254)
10         LINESTRING (12.052304491086 55.634166861254, 12.052365253079 55.634249843067)
15         LINESTRING (12.052365253079 55.634249843067, 12.052145214983 55.634301401561)
10         LINESTRING (12.052145214983 55.634301401561, 12.052084452911 55.634218419749)

This function can be accessed by making a donation and downloading the SQL Server Spatial package of TSQL functions.