Archive for the 'GIS in general' Category

Oct 30 2010

Tied up in knots with 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, 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).

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

Oct 20 2010

Oh say, can you SQL (Spatial)?

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!

No responses yet

Sep 26 2010

Meandering into mobile map delivery

Published by under GIS in general,google earth

Opensim has slipped to a back burner for these past few weeks, as I button up a new edition of county terrain–this one fused topography and bathymetry.
Another focus has been to compile the map data layers necessary to publish a local edition of the large-scale topographic map, using the ESRI Community Maps Program template. This (Opensim-distracting) work forms a synergy, as the terrain drives revised hillshade and topo-bathy contour lines, as well as forming the basis for refined surface flow and stream mapping, all of which show up in the large-scale topo map. Also, I’ve compiled a great deal of web research into a single feature layer: classified building footprints, many of which have use categories [commercial, medical, educational, religious, etc.] and specific organization names included, for every structure on the associated parcel.

Why generate this large-scale topographic base map? It’s purpose is to provide local control over local detail, in a cloud-based service suitable for desktop, tablet, and mobile phone users who may consume applications that mash or overlay the base map with relatively sparse point or line features. While we’ll be able to serve up the base map through our own MarinMap servers, one of the key benefits of the conforming large-scale topographic map style (as made possible by the template) is that the caches of the four largest-scale views will be rendered and contributed for inclusion in the ArcGIS Online map services, as viewed at and consumable by many interactive map sites. The cached map tiles are essentially raster prints of the map template (which has styles customized at four local zoom levels) that can be served extremely fast.

Those fast maps can also be efficient, for mobile phones and tablets consuming base maps with no special rendering. Along that path, I’ve rebuilt my mobile phone this past week, taking the T-Mobile G1 and rooting it, then installing various new images until I found one that was up-to-date, stable, and attractive to use. As of this weekend, I’m now running Chromatic 4.5 with SetCPU 3.0.2 overclocking. Chromatic 4.5 is a packaging of CyanogenMod 6.0.0 with the ADW launcher and many other tunes. CyanogenMod 6.0.0 is an Android Open Source Project-based build of Android version 2.2 “Froyo”. So now, when I view maps and browser content, I can pinch to zoom out, spread to zoom in, and still push around the map to pan it on the screen. As part of sharing content with mobile phones, I’ve learned how to generate QR tags to provide access to map URL or geolocation tags.

For fun today, we visited the Silicon Valley area to show a guest some of the sights. I hadn’t visited the Googleplex in a couple of years, and many of the adjacent buildings that I recall as having different companies in them are now part of the Google campus. We found one interesting sight, too–and it appears that today may have been Google Inc.’s 12th birthday. We didn’t see any partying in the mostly-empty parking lots, which wasn’t surprising for a Sunday!

At the Googleplex 2010 09 26

No responses yet

Aug 27 2010

Marin County Topo-Bathy Surface “tbsm45cm” Released as 2010.08

Blemishes notwithstanding, nearly six months of back-burner work has reached a threshold of readiness and is outward bound to some engineering firms, flood mappers at FEMA, and interested parties within the county. A handful of known issues remain unresolved. Proper name is “tbsm45cm_20100823”, proper edition is “2010.08”.

This is the third version of the terrain. Second version was “2010.01” and included multiple LiDAR data sets, but fewer than presently used, and was a topographic model only. First version was “2009.09” and was mainly photogrammetry and FEMA LiDAR, and was the last version to be developed in California Coordinates. Once the massive NCALM LiDAR data sets were processed, it became easier to move everything into WGS84 UTM zone 10 north meters projection, WGS84 NAVD88 CONTUS Geoid 2003 vertical position.

The NOAA utility program VDatum, a brilliant Java-based application able to stream-process data sets of near-infinite size, brought the NCALM data to heel, and opened up decades of NOAA depth surveys to our use in integrated topographic-bathymetric surface modeling.

First-return NOAA ALACE LiDAR swaths were fused along the outer coast, as bare-earth filtered versions were not produced in 1997–2002; the benefits of LiDAR detail along the rocky coast do seem to outweigh the distracting appearance of structures near Rodeo Lagoon, Stinson Beach, and outboard Bolinas.

When ArcGIS 9.4 beta 2 reached its limit in ability to render the terrain dataset into 45cm grid over the full extent, the clipping quadrants created to resolve this problem ended up chopping a very small portion of Sonoma county that drains into Estero Americano; the full watershed remains intact in the 1-meter version of the terrain grid under analysis for county-wide hydrology. Likewise, the tighter clipping quadrants lost a few hundred meters of San Pablo Bay bathymetry just west of where Marin, Sonoma, Contra Costa, and Solano counties meet. Also, tighter clipping quadrants snipped a portion of the San Francisco Bar southerly of San Francisco’s Seal Rock that was intended to be part of the model. All of these areas exist in the 100cm grid, and will be part of drainage analysis.

Happily, we have updated the workstation to ArcGIS 10, and have been enjoying such great speed gains with Spatial Analyst that our ERDAS use has been noticably reduced. Finally, Spatial Analyst is often showing performance nearly on par with ERDAS. Thank goodness that the Raster Calculator survived the transition to version 10 ArcGIS!

Painfully, the existence of unutilized bathymetric data sets for upper broad-channel Corte Madera Creek and Bolinas Lagoon have been revealed this week. Hey, there’s already something to look forward to for the next build!

The new terrain is getting some immediate use in support of an effort to participate in ESRI’s Community Maps Program for large-scale topographic mapping. The Program provides a template geodatabase with 36 vector feature classes and two raster, into which local agencies may pour their data. Once tucked into a conforming schema, a template multi-scale map document is provided with 120 layers—30 at each of four large scales that correspond to Google Maps and Bing Maps projection and cache tiling schema. The difference is that the template document makes use of ESRI tools to allow much more local detail to be packed into a map designed with notably more sophisticated cartography than either Google or Bing maps now have. The Community Maps Program concept is that local agencies may publish their local detailed content in a fairly uniform style, while retaining a world-wide seamless context for their surrounding area.

Qualitatively, the effect is that, when viewing the topo map alongside either Google or Bing maps (on two monitors, with comparison made at the same scale), the map looks to be a larger scale. It isn’t, and I’ve measured the size of features to convince myself, but my mind insists that I’m zoomed farther in on the map for some reason. My guess is that it is a perceptual effect of the much greater amount of information that is cleanly displayed in the map versus the much sparser Google and Bing content at these large zoom levels. Try it out—it’s like a carto version of an optical illusion!

The 120 layers in the template large-scale topographic base map from the ESRI Community Maps Program are arranged to provide four precise cartographic designs for Google/Bing map cache levels 16 through 19, which correspond to these display scales
1:15000–1:6001 (level 16, a.k.a. ~9k)
1:6000–1:3501 (level 17, a.k.a. ~4.5k)
1:3000–1:1501 (level 18, a.k.a. ~2k)
1:1500–1:501 (level 19, a.k.a. ~1k)
One of the most attractive areas currently online is Toronto, ON where at levels 18 and 19, individual building outlines are graced with street addresses.

Anyhow, the new tbsm45cm model will serve County of Marin’s effort at large-scale topographic mapping several ways. First, it has made possible a very detailed hillshade that helps emphasize the grading around each hillside structure in the county. Second, it helps us to create the required metric topographic contours. These are necessary to meet world-wide mapping standards, and throughout this weekend, contours are being generated from a related (smoothed version) of the terrain on 50cm vertical interval. Needless to say, most of these won’t get used in the map renderings, but the ESRI cartographers have shared a very clever indexing scheme that will help us use this single set of metric contours to support the requirements for all four of our topographic map scales.

No responses yet

Feb 04 2010

A visit to ScienceSim Geography regions – OpenSim with turbo boost

I don’t have much to say about these regions that hasn’t been written already, and my views have been less aesthetic than Shenlei’s.
But in the interest of boosting the bandwidth by which I can share OpenSim, I’ve invested in a much newer Adobe Premiere Elements than I’d been using for the past five or so years. It’s a gas to have it multi-thread while rendering, and I have direct-to-FLV write. Trying to share as much of the motion and fidelity via YouTube as I possibly can, I’ve crafted a video resolution that is a multiple of my Hippo / SL viewer screen. The FRAPS video direct to AVI (sorry, it’s Win XP) is 1600 x 1140 @ 10 fps. Yup, those are video frames. In the interest of surviving an upload, I’ve rendered them highly time-compressed, with output at 1515 x 1080 @ 15 fps. As of tonight there’s no sound, no intertitles, just the rushes.
oops, if I read the YouTube Instructions for best formats, I should have trimmed the width to 1440, which is a multiple of 16.
Also, I have more direct upload options now with Premiere 8 than I had with my (recently demised) copy Premiere 1.0. Go Figure ;^)

While the Windows box grinds out the video print, I’m over here on Ubuntu blogging in a tab of 64-bit Chrome and it is fine & fast.

For these videos, I visited ScienceSim Geography22_44 region and set the view to wide angle, then sat up at about 500 meters and watched the regions rez their terrain. For some folks, it will rank right up there with watching lead-based paint dry. For geography folks I’m hoping that these few minutes of sped-up video will convey, by dogged repetition, the primacy of regions in the provision of virtual environment simulators.

By the way, I’ve got a task: I need to find a better buzz word for the GIS community. I’ve been advised by some serious and well-intentioned (not to mention well-informed) folk that terms like “virtual” and “immersive” are actually boring to GIS’ers. So I’ll need to think about how to convey the concepts of “Mirror World”, “Multiuser Virtual Environment”, “Immersive Connected Experience”, “Third-Person Virtual World”, and related concepts into a catchy moniker. Hopefully, one that is not presently trademarked, either!

I’m trying to remain serious about this, but some of the options are treacherous. Geography in Social Media has a possibly awkward acronym; maybe it can be saved in recognizable form as “GIS for Social Media” or “Geography for Social Media”: GFSM
The term “3D Map with Me” is terse, slightly ambiguous

Here is the video chopped as it was when uploaded with 1515 x 1080 resolution. Problem with that is that by not preserving dimensions at a multiple of 16, and saving my viewer’s aspect ratio rather than the (standard since 16mm film) 4:3 aspect, my upload is clobbered into something perhaps suitable for a smartphone. So please consider this the Smartphone Version of last night’s rushes:

Then, once again with feeling, or at least with a little more rest, there is what I hope to be an HD-friendly moving vision of OpenSim, as it appears on the ScienceSim Geography regions. Yes! After it ripened on the YouTube servers for a few hours, I now see all the higher-res versions available. At 1440 x 1080, this is pretty close to what I see on my screen with a live Hippo viewer.

And after a day’s cogitation: anyone care to comment on the term: “Social Immersive Media GIS” as a moniker? Oops — I used “immersive” 8^(

No responses yet

May 24 2009

Some thoughts on geography

updated 2009 05 26

I’ve been waiting for some property boundary issues to resolve in SL, and it’s sort of pitiful to see how long that can take.  It’s with ever more regret that I find myself on the Mainland.  But that hasn’t kept lots of real-world interesting stuff from taking shape.

The following video is not new.  In fact it’s about a year old, but somehow I hadn’t seen it until tonight and I found it somewhat encouraging. Thanks for O’Reilly and Where 2.0 for bringing these two on stage together!

And the following pean to Google Earth did inspire me, personally. Hey, I was reading road maps at 5, covered my wall completely with National Geographic maps at 10, learned to navigate with nautical charts at 12, read aeronautical charts and completed an urban planning project at 14. Sometimes, it’s fun in rare moments when it’s dark overcast and I’m in an exotic place for the first time and I don’t know the way north; more often, I’ll savor the feeling of knowing which way is north while dreaming.

Meanwhile, back at the lab, the global set of county terrain is being compiled into an ESRI Terrain Dataset. This will include over 360 million masspoints, merging both interpolated 2-foot interval contour vertices together with FEMA LiDAR mass points, plus break lines and waterlines from photogrammetry. The goal is to use the ESRI Terrain data as a format to stage everything together to produce 30cm grid interval DEM in the urban areas. With luck, we’ll have that ready about the time that the latest photo mosaic finally gets loaded into ArcSDE successfully. Maybe grids from the Terrain can help create very detailed 3D county models. Hey Wei – we still have inverted terrain in Google Earth at the quarry on San Pedro Point! ;^)

No responses yet

« Prev