Geocoding Points Against a Grid of Cells/Tiles With/out Spatial Operators
Southern Tasmania
Let’s now construct an more realistic example that uses a rotated grid in a real place near home!
We will endeavor to discover the comparative performance of these two methods, though straight comparison depends on real situations in which each one solution may be best but the other not. (See conclusion.)
Now let’s use a different function to create an angled grid in a real world location.
/****f* TILING/STTileGeomByPoint (2008) * NAME * STTileGeomByPoint -- Creates mesh of tiles in geometric space anchored to a supplied origin. * SYNOPSIS * Function STTileGeomByPoint ( * @p_point geometry, * @p_origin varchar(2), * @p_numTileX integer, * @p_numTileY integer, * @p_TileX float, * @p_TileX float, * @p_rAngle float, * @p_AsPoint bit * ) * Returns @table table * ( * col Int, * row Int, * geom geometry * ) * DESCRIPTION * This function generates a mesh (grid) of tiles anchored to the supplied anchor point. * The mesh of tiles is controlled by five parameters: * 1. An anchor point * 2. Whether the anchor point is at the LowerLeft (LL), UpperRight (UR), UpperLeft (UL)L or LowerRight (LR) of the grid * 3. The number of tiles in X and Y direction; * 4. XY tile size in meters; * 5. Optional rotation angle (around origin/anchor point) * INPUTS * @p_point (geometry) -- Starting Point for grid (Upper Left) * @p_origin (varchar) -- Position of point wrt grid: LL,UL,LR,UR * @p_numTileX (integer) -- Number of tiles in X (longitude) direction * @p_numTileY (integer) -- Number of tiles in Y (latitude) direction * @v_TileX (float) -- Size of a Tile's X dimension in real world units along parallel of Latitude (ie X distance) * @v_TileY (float) -- Size of a Tile's Y dimension in real world units along meridian of Longitude (ie Y distance) * @p_rAngle (float) -- Optional rotation angle from North. * @p_AsPoint (bit) -- Return Tile as point or polygon * RESULT * A Table of the following is returned * ( * col Int -- The column reference for a tile * row Int -- The row reference for a tile * geom geometry -- The polygon geometry covering the area of the Tile. * )
Create and Store a Grid
So, let’s create, and store, an angled grid, with a Morton key as its primary key, and whose anchor point is the LL of the grid.
drop table if exists dbo.grid_cells
go
with data as (
select -100.0 as rAngle,
'LL' as anchor_point,
1000.0 as tileX,
1000.0 as tileY,
15 as numTiles, -- 8x8 = 64
64 as PointsPerTile, -- 8*8*64 = 4096
geometry::Point(537038,5215350,28355) as origin -- LL
)
select [dbo].[STMorton](t.col,t.row) as morton_id, col, row, geom
into dbo.grid_cells
from data as a
cross apply
[dbo].[STTileGeomByPoint] (
a.origin,
/*p_origin*/ a.anchor_point,
/*@p_numTileX*/ a.numTiles,
/*@p_numTileY*/ a.numTiles, /* 64 cells */
/*@p_TileX */ a.tileX,
/*@p_TileY */ a.tileY,
/*@p_rAngle */ a.rAngle,
/*@p_AsPoint*/ 0
) as t;
alter table dbo.grid_cells alter column morton_id int not null;
go
alter table dbo.grid_cells add constraint grid_cells_pk primary key clustered (morton_id) ;
go
CREATE SPATIAL INDEX [grid_cells_spx] ON [dbo].[grid_cells]([geom])
USING GEOMETRY_AUTO_GRID
WITH (BOUNDING_BOX =(200000, 5100000, 500000, 5500000))
ON [PRIMARY]
GO
The generated grid looks like this.
Create And Store Sample Points – Fish, Rocks, Boats?
We can now will generate a set of points for our angled grid which can be used to compare the two methods.
Note that we cannot use the cells to directly create sample points because the grid is angled (it is easy if it is not). The approach the table creation code below uses is to take all the grid’s cells and merge them into a single object; then compute the single polygon’s envelope and then use it to generate a set of candidate points. Once the candidate points are created, merge them back to the merged single polygon keeping only those that fall within it.
drop table if exists dbo.grid_points;
go
with data as (
select -100.0 as rAngle,
'LL' as anchor_point,
1000.0 as tileX,
1000.0 as tileY,
8 as numTiles, -- 8x8 = 64
96 as PointsPerTile, -- 8*8*96 = 6144
geometry::Point(537038,5215350,28355) as origin -- LL
)
, mbr as (
-- Generate single tile to create area for creating points
select mbr.minx, mbr.maxx, mbr.miny, mbr.maxy
from (select geometry::UnionAggregate(a.geom) as geom
from dbo.grid_cells as a
) as t
cross apply
[dbo].[STGeometry2MBR] (t.geom) as mbr
) -- Now create points in this aggregate polygon (will remove that which are not needed later)
, points as (
SELECT gs.IntValue as point_id,
geometry::Point(
[dbo].[STRandomBetween](m.minx, m.maxx),
[dbo].[STRandomBetween](m.miny, m.maxy),
a.origin.STSrid
) as geom
FROM data as a CROSS APPLY [dbo].[generate_series](1,a.numTiles * a.numTiles * a.PointsPerTile,1) as gs,
mbr as m
)
SELECT f.point_id, f.geom
INTO dbo.grid_points
FROM (SELECT p.point_id,
row_number() over (partition by p.point_id order by p.point_id) as low_point_id,
p.geom
FROM points as p
inner join
dbo.grid_cells as t
on (t.geom.STIntersects(p.geom) = 1) -- select only point that fall within single polygon
) as f
WHERE low_point_id = 1
GO
alter table dbo.grid_points alter column point_id int not null;
go
alter table dbo.grid_points add constraint grid_points_pk primary key clustered (point_id) ;
go
CREATE SPATIAL INDEX [grid_points_spx] ON [dbo].[grid_points]([geom])
USING GEOMETRY_AUTO_GRID
WITH (BOUNDING_BOX =(200000, 5100000, 500000, 5500000))
ON [PRIMARY]
GO
select count(*) from dbo.grid_points;
GO
-- 3263 out of 6144
These points can be visualised over the top of our grid as follow.
Note that there is a reduction in the number of points. They were generated for a square polygon over the angled grid; that mismatch causes the point reduction.
Comparison
We will now compare the two methods with our real world grid.
First the SQL Server Spatial (point in polygon) method.
Note that the result will include the two morton key values side by side, so we can visually compare the two morton keys (or we could filter as well).
-- Note, that the table has a spatial index.
select t.morton_id,
[dbo].[STMorton](t.col,t.row) as morton_point_id,
p.geom
from dbo.grid_points as p
inner join
dbo.grid_cells as t
on (t.geom.STContains(p.geom) = 1)
go
-- 3263 Rows in 3 seconds.
Next the direct reference approach.
with pointsWithRef as (
-- Create Morton Key from generated col/row references.
SELECT [dbo].[STMorton](r.col,r.row) as morton_point_id, r.col, r.row, f.geom
FROM dbo.grid_points as f
CROSS APPLY
[dbo].[STTileColRowRef4GeomPoint](
f.geom,
geometry::Point(537038,5215350,28355),
1000.0,
1000.0,
-100.0
) as r
)
select t.morton_id, morton_point_id
from pointsWithRef as p
inner join
dbo.grid_cells as t
on (t.morton_id = p.morton_point_id)
go
-- 3263 Rows in 7 seconds
Even with STContains’s best case use, SQL Server Spatial’s operator was still 2.3 times (7 sec / 3 secs) or 233% faster.
If we look at this from the point of view of speed / geometry we get:
- 7 / 3263 = 0.002 seconds per feature, and
- 3 / 3263 = 0.0009 seconds per feature.
Both are pretty quick.
But is the difference really that bad?