This post is an extension of ”SQL Server spatial query, part 1”. The first post primarily focused on how the basic query and geometries are put together and the types of results you should expect to see. In part 2, the goal is to wrap some of these techniques into a SQL procedure to evaluate if an input feature intersects a reference layer.
Suppose we wanted downloaded some weather data from NOAA and store in our SQL Server database. In many cases this data is available in a vector format (point or polygon), or it’s available in grid/raster format. An efficient way to store this weather information is to store vector geometry only once (with a unique ID), and every time we download the data we store the weather value and the related geometry ID instead of storing the geometry every time we download. This way we can also store all historical data without the bloat of storing the redundant geometry and still be able to perform quick statistical analysis at any location(s) directly in query builder or with a SQL View. Since the ID of each feature from NOAA is unknown (in most cases), or can change, we can use a spatial overlay to determine an ID from our own lookup table (spatial reference layer).
Side note: This also works great with raster/grid data even thought it isn’t directly supported in SQL Server 2008 - I just convert the raster cells to points or polys before running the analysis using a Python script (something for another post).
For this scenario, let’s assume we downloaded a set of polygons in a gridded polygon format and we want to determine if each polygon intersects our reference point layer. If it does - we will store the new value and related reference point ID.
The SQL Server spatial query
Skipping over the many ways you can download data from NOAA and processing into many different formats (Shapefile, WKT, etc), let’s just jump forward to the part where we have thousands of WKT polygons with attributes: a date attribute and current weather parameters. If you need some help getting to this stage, you might want to checkout NOAA’s ”Weather and Climate Toolkit” and use a little Python magic.
Similar to part one we want to compare an input geometry, a WKT polygon, with our reference point layer. As we loop through each polygon record we will send the data to the following SQL procedure:
|-- Name: SQL Server Procedure with STIntersection|
|-- Author: Bryan McIntosh|
|-- Description: Populated data table from SQL Server Spatial query|
|CREATE PROCEDURE [dbo].[prc_NoaaDataLoad]|
|@dt datetime, --DateTime from NOAA header info|
|@val float, --Value of weather data (mm of rain, temperature in F/C, etc)|
|@wType varchar(20), -- Type of data stored (precipitation, temperature, etc)|
|@wkt varchar(8000) --polygon shape (geography type)|
|SET NOCOUNT ON; --added to prevent extra result sets from interfering with SELECT statements.|
|DECLARE @polygeo GEOGRAPHY;|
|SET @polygeo = GEOGRAPHY::STGeomFromText(@wkt,4269); --4269 is NAD83 GCS|
|INSERT INTO tbl_NoaaData (PointID, LogDate, DataValue)|
|SELECT POINTID, @dt, @val, @wType|
|WHERE tbl_ReferencePointsGCS.shape.STIntersects(@polygeo) = 1|
The procedure includes a “@wType” parameter just in case we want to download multiple types of NOAA data and still be able to use the same proc. Since I pre-process the NOAA download using python, I ca send data directly to the SQL procedure in Python (w/ PYODBC):
MyStatement = "exec prc_NoaaDataLoad '" + str(dt) + "', " + str(val) + ", '" + "Temperature" + "', '" + wktGeom + "'" cursor.execute(MyStatement)
This code is executed for each row in the NOAA layer, and stores the data for each record (in tbl_NoaaData) when the polygon from NOAA intersects the reference point table. The WHERE clause determines if there is an intersect (=1) versus no intersect (=0). So if the result of the overlap = 0, the row isn’t loaded and we move on to the next. Seems so… simple - and it is.
Even though I’m using Python for this process, I should mention that I’m not using any ArcPy or other GIS module. Just straight python calling a sql server spatial query. This also has the added benefit of running extremely fast (< 6 seconds), and not requiring an ArcPy/ArcGIS license.