Darb Dabney's 3d GIS – OpenSim

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

"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

(With apologies to Bill Ham and ZZ Top)

Some progress, and apparently some shifting ice floes underfoot.  Some blog posts say that the SQL Server Management Studio (SSMS) 2010 Express does not feature a “Spatial Results” tab for queries.  In fact, at least this week, it does.  It does seem like full-bore database development tools are withheld  by Microsoft until the Visual Studio Premium 2010 tier, and yet I’m  presently unable to locate  an IntelliSense add-in for the OGC-compliant methods on Geometry instances.

Here they are, 43 OGC methods on geometry instances,
plus 16 OGC static geometry methods,
along with the typical MS-style augmentation
11 augmented MS-extended methods on geometry,
and 4 augmented MS-extended static geometry methods.

Another interesting topic is how to merge certain feature classes for use by both ArcSDE and SQL Spatial.  In this regard, since ArcGIS Server 9.3 it appears that one can have ArcSDE store feature classes with SQL Server 2008 GEOMETRY data type.  Here’s an ESRI documentation page on the subject.

Sometimes, it’s the simple stuff that one stumbles over (at least I do).  Here’s some of the bootstrapping stuff for editing SQL in the Visual Studio 2010 environment.  Plus, although I’m personally comfortable with programming languages like structured C, and have been trained in (but not worked professionally with) Objective C++ and Visual C++ (at 6.0), I’m not finding immediate transparency with C#; most of the Visual Studio Web Developer  tutorial examples use Visual Basic.

1) How to Start the T-SQL Editor
2) The intro videos for C# development
3)

Many interesting projects have been happening, so many that projects have been backed up while workstations grind through their days-long work flows. This creates opportunities to update systems while we wait, and so the past week has seen massive updating of Windows system stuff and deleting old development tools so that the decks were cleared for MS SQL Server 2008 R2 Management Studio. With a bit of shape2sql from SharpGIS and today we’ve had our first spatial tables loaded and fledgling spatial SQL queries made.  Oh, at the end of the day we got the third of five workstations updated to ArcGIS Desktop 10 and the full complement of seats moved over to ArcGIS 10 licensing.

What’s been holding things up has been some final processing of cartographic-grade contours, one of the key new feature classes developed for support of the ESRI Community Maps Program, where local jurisdictions grind their own cache tiles for large-scale topographic mapping.  Since it’s part of a worldwide seamless map, of course the contours and spot elevations must be metric—which meant generating new contours.  The bathymetry was done on one-meter interval, with some half-meter supplemental contours in shallower waters.  The topography was done at quarter-meter intervals up through 25 meters elevation, where we’d included all available LiDAR data, the half-meter intervals through 50 meters, and one meter up to the summits.  All non-integer meter contours were flagged as supplemental, and contour index attributes were also calculated for 2, 5, 10, 20, and 50-meter intervals.   Spot elevations were derived from VARGIS photogrammetric spot elevations that were screened down from about 75,000 points to just 440 points, using neighborhood focal statistics.

For performance’s sake, after index attributes were calculated, the contours were chopped into segments not-to-exceed 500 meters shape length.  When all the chopped segments from all the various elevation ranges had been merged to a single polyline feature class in a feature dataset in a file geodatabase, spatial indices were built for the cache tile scales about 1:1200, 1:4800, and 1:19000.  The purpose is to make the contours as fast as possible not just for web app interaction, but for rendering the large-scale cache tiles.  In the end, the contours in this one feature class have about 1.4 million polyline features that fill the file geodatabase to about 6.8 GB of data.

Meanwhile, some great brainstorming has taken place with regard to hydro-enforcement of the terrain model, so that accurate synthetic flow lines can be generated county-wide.  Based on prototype work performed by Evan K. Babb before summertime, we knew how to deal with our terrain features and manually generate hydro-enforcement features.  Now, with a funded project to do the work county-wide, we’ve needed to devise a spatial analysis technique using the approach described by Poppenga and others in 2009.  Based on some correspondence and personal communication with Poppenga, we have sketched a workflow that should automate the creation of most simple (single road-crossing) hydrologic enforcement features, especially in the absence of accurate or complete spatial features for culverts.

So with the contours completed, it should be possible to have some progress made on automation of hydro-enforcement, completing more of the data development for large-scale topographic base mapping, help my colleagues continue the ArcGIS 10 migration, and follow up with some demonstration spatial SQL queries.  It’s an exciting couple of weeks!

There’s also been a thread of thought over the past six years or so that supports the workflow of using Visual Studio to compose SQL queries, in preference to the direct SSMS (SQL Server Management Studio, Tool Formerly Known As “Enterprise Manager”) ways of composing a query, because the typed pattern matching (MS: Intellisense) pulls up all the various ST__  spatial methods that are available in SQL Server 2008 spatial, but saves one needing to look-it-up and type-it-in correctly.  Such productivity gains really help boot up the learning experience a bit faster!