Oct 30 2010

Tied up in knots with Spatial SQL

Published by at 2345h under Blog mentions,GIS in general,Spatial SQL

The rules for geometry validation are just different between the ESRI world and the OGC world as instantiated by Microsoft in SQL Server 2008 R2.
I’ve been pounding my head against the wall over what I thought were the challenges of taking fine examples from Alastair Aitchison (a.k.a. tanoshimi on the MSDN forums) that are presented in Beginning Spatial with SQL Server 2008, and trying to make these happen in useful ways on real GIS data.  That is to say, data that were not created by hand-entry or even well-known-text (WKT) entry into SSMS, which are often  simple single features.

I was particularly concerned reading the first paragraph in Chapter 11: “All of the techniques discussed in this chapter apply only to a single item of geography or geometry data, and generally do not require any parameters.”  As I understand now, Mr. Aitchison meant “All of the techniques as discussed…” and all the various Spatio-Temporal methods like STIntersects() work perfectly well on perfectly valid spatial tables.

I thought that these methods just weren’t going to work on my GIS data,  typically large tables with lots of objects, not the scalar objects of the book examples.  But as it turns out, I had a different problem.  My geometry, as transcribed via Shapefiles and the Shape2SQL tool from SharpGIS.net, was arriving in a form that was occasionally invalid for certain features.  One invalid feature in a spatial table makes the entire SQL statement crash.  Errors were reported with reasons that did not help me find my actual problem very easily.  Bad geometry meant that I needed to clean it up by running the MakeValid() method.

After a couple of hours goofing around with T-SQL syntax, I’ve made my first baby steps toward cleaning up my stuff.  Here’s how that is happening:

(select ID, PARCEL, Prop_ID, Shape_area, Shape_len, (CASE Shape.STIsValid() WHEN 1 THEN Shape ELSE Shape.MakeValid() END) as Shape INTO spatial_01.dbo.Parcel2 FROM spatial_01.dbo.parcel)

As so often with code, it looks stupid-simple to see it written out in the way that actually works; a major “duh” factor.
The pay-off is that with 100% valid geometry, table vs. table spatial queries are now working.  I’ve not yet broken any speed records in query performance vs. ArcGIS 10, and yet I have finally proven to myself that I can use  Spatial SQL to do powerful geoprocessing that, in some cases, will be easier to write and maintain in SQL than in ArcGIS / Model Builder / Python.

For context, I’ve struggled through the past four days trying to figure out where my Spatial SQL syntax was messed up.  As it turned out, exactly ZERO of my uploaded spatial data tables arrived with 100% valid GEOMETRY data.  Anyway, once I made everything OGC-valid, a query like this

SELECT NAME, b.PType INTO spatial_01.dbo.School_Geology
FROM spatial_01.dbo.school2 a, spatial_01.dbo.geology2 b
WHERE a.Shape.STIntersects(b.Shape) = 1
Is really all it takes to return the mapped geologic unit beneath the point for every school in the county.
For 120 school points, and 3700 geology polygons, the query above returns its result in about 30 seconds, with no tuning of the spatial indices and no SQL hints provided (yet).
"NAME","PType"
"TOMALES HIGH SCHOOL","Twg"
"BOLINAS-STINSON SCHOOL (STINSON)","Kfs"
"LAGUNA JOINT ELEM. SCHOOL","Qal"
"SAN GERONIMO VALLEY ELEM. SCHOOL","fsr"
"LAGUNITAS ELEM. SCHOOL","fsr"
"MARIN HORIZON SCHOOL","fsr"
"COUNTY COMMUNITY SCHOOL","Qal"
"ST RITA SCHOOL","Qal"
"ST PATRICK SCHOOL","fsr"
...

On a last note, this week I have validated the experience of using ESRI ArcSDE 9.3.1 settings to configure SDE (when running atop MS SQL Server 2008) to store its geometry in native OGC GEOMETRY binary form rather than the classic SDEBINARY form.  To do this, someone with db owner privileges can update the SDE_dbtune table with these directions to cause the GEOMETRY format to be used by default.  In principle, this should allow one to exploit the live feature classes in SDE for Spatial SQL queries—without the muss and fuss of exporting through Shapefiles and Shape2SQL (which is not intended for production-environment use, anyway), or need to code something in .NET to export a copy.

Wouldn’t it be nice to have ArcSDE running just like it always has, and then also be able to use Spatial SQL queries on the exact same feature classes that are being used throughout the enterprise?  Sounds great, in principle.  After tonight’s insight into the need to MakeValid() some perfectly functional feature classes, I do have some concerns that SDE feature classes could get broken by MakeValid().  On the other hand, it appears that one can take a Spatial SQL GEOMETRY-based table that meets SDE limitations: (1) only one geometry column, (2) all rows of the table must have the same geometry type, (3) all rows of the table must have the same Spatial Reference ID (SRID) defined—and any SDE feature class that was cleaned through MakeValid() should retain those characteristics—and then register the resulting spatial SQL table once again as an SDE feature class.  Hey, it sounds like it’s worth a try!

More working through Aitchison’s book has yielded many successful uses of OGC methods on our GIS data:

-- SELECT Shape.STGeometryType() FROM [spatial_01].[dbo].[road]
-- SELECT Shape.STGeometryType() FROM [spatial_01].[dbo].[Parcel]
-- SELECT Shape.STGeometryType() FROM [spatial_01].[dbo].[school_pt]
-- SELECT Shape.STGeometryType() FROM [spatial_01].[dbo].[geology]

-- SELECT Shape.STDimension() FROM [spatial_01].[dbo].[geology]
-- SELECT Shape.STDimension() FROM [spatial_01].[dbo].[road]
-- SELECT Shape.STDimension() FROM [spatial_01].[dbo].[school_pt]

-- SELECT Shape.InstanceOf('Surface') FROM [spatial_01].[dbo].[geology]
-- SELECT Shape.InstanceOf('Polygon') FROM [spatial_01].[dbo].[geology]
-- SELECT Shape.InstanceOf('Curve') FROM [spatial_01].[dbo].[road]
-- SELECT Shape.InstanceOf('LineString') FROM [spatial_01].[dbo].[road]
-- SELECT Shape.InstanceOf('Point') FROM [spatial_01].[dbo].[school_pt]

-- SELECT Shape.STIsSimple() FROM [spatial_01].[dbo].[road]
-- SELECT Shape.STIsSimple() FROM [spatial_01].[dbo].[BldgFoot]

-- SELECT Shape.STIsClosed() FROM [spatial_01].[dbo].[BldgFoot]

-- SELECT Shape.STIsRing() FROM [spatial_01].[dbo].[road]

-- SELECT Shape.STNumPoints() FROM [spatial_01].[dbo].[road]
-- SELECT Shape.STNumPoints() FROM [spatial_01].[dbo].[geology]
-- SELECT Shape.STNumPoints() FROM [spatial_01].[dbo].[school_pt]

-- Find Centroids and blob them out big
-- SELECT Shape.STCentroid().STBuffer(3500) FROM [spatial_01].[dbo].[City]

-- Calculate Perimeter Lengths and sort them out
-- SELECT Shape.STLength() as Shape_len FROM [spatial_01].[dbo].[City] ORDER BY Shape_len

-- Calculate City areas and sort them out
-- SELECT NAME, CAST (Shape.STArea()/1000000 AS INTEGER) as Shape_km2
--   FROM [spatial_01].[dbo].[City] ORDER BY Shape_len
-- Reading the SRID
-- SELECT Shape.STSrid FROM [spatial_01].[dbo].[City]

-- Setting the SRID
-- UPDATE [spatial_01].[dbo].[City] SET Shape.STSrid = 32610
-- SELECT Shape.STSrid FROM [spatial_01].[dbo].[City]

-- oops, that's wrong, let's set it back to 2768
-- UPDATE [spatial_01].[dbo].[City] SET Shape.STSrid = 2768
-- SELECT Shape.STSrid FROM [spatial_01].[dbo].[City]

-- counting elements
-- SELECT SUM(Shape.STNumGeometries()) FROM [spatial_01].[dbo].[City]
-- SELECT SUM(Shape.STNumGeometries()) FROM [spatial_01].[dbo].[school_pt]
-- SELECT SUM(Shape.STNumGeometries()) FROM [spatial_01].[dbo].[geology]
-- SELECT SUM(Shape.STNumGeometries()) FROM [spatial_01].[dbo].[road]
-- SELECT SUM(Shape.STNumGeometries()) FROM [spatial_01].[dbo].[Parcel]
-- SELECT SUM(Shape.STNumGeometries()) FROM [spatial_01].[dbo].[BldgFoot]

-- Geometric Difference between cities and their own 3500-meter centroid blobs
-- SELECT Shape.STDifference(Shape.STCentroid().STBuffer(3500)) FROM [spatial_01].[dbo].[City]

-- Geometric Symmetric Difference between cities and their own 3500-meter centroid blobs
-- SELECT Shape.STSymDifference(Shape.STCentroid().STBuffer(3500)) FROM [spatial_01].[dbo].[City]

-- Creating 50-foot Buffer around roads
-- SELECT Shape.STBuffer(50) FROM [spatial_01].[dbo].[road]

-- Simplified Buffer around roads, tolerance 8 feet
-- SELECT Shape.BufferWithTolerance(50, 8, 'false') FROM [spatial_01].[dbo].[road]

-- Create convex hull around schools
--SELECT Shape.STConvexHull() FROM spatial_01.dbo.school_pt

-- buffer parcels
-- SELECT Shape.STBuffer(300) AS Shape INTO pcl_seed FROM parcel where Prop_ID = '001-001-01'
-- select * from pcl_seed  -- to verify the selection
-- SELECT Prop_ID, a.Shape FROM spatial_01.dbo.parcel a, spatial_01.dbo.pcl_seed b
--   WHERE a.Shape.STIntersects(b.Shape) = 1

No responses yet

Trackback URI | Comments RSS

Leave a Reply

You must be logged in to post a comment.