Category Archives: GIS

Geographic Information Systems

That’s “Code” Comfort

I’ve been working on a mobile app that will support adding vector-based mapping polygons. Using California and Massachusetts to do a test on MULTIPOLYGON states and a map extent spanning the lower 48, I imported the data, ran the code and when I saw the visual, thought something was messed up:
How could the code be putting California at the top of the image? I went back to my source map and sure enough, the northern tip of California really does go as high as central Massachusetts. It occurred to me that it could be because of the Lambert Conformal Conic projection for the contiguous U.S. that I’m using. But it turns out the northern boundary of California is 42° and so is the southern boundary of western Mass (or at least it seems that was the intention at the time, given the quality of surveying and certain conditions facing surveyors – see the Granby Notch story (http://www.nytimes.com/2001/01/26/nyregion/driveway-another-state-blunder-1642-creates-headaches-today-for-homeowners-who.html) .

File:Masscolony.png

Incidentally, several states share 42° as an official boundary:

File:42nd parallel US.svg
Sometimes our notions of things, based on certain conceptions, distort reality. I can see where my perceptions of California as being in the southwestern part of the country and being mostly warmer than the chilly, northeasterly Massachusetts would make it hard to believe they could have boundaries near the same latitude. At least I know my routines work. That’s “code comfort.” Still lots of work to do!

Burn Map

https://maps.google.com/maps/ms?msid=217066668344035014417.0004df0ab0046cd029a58&msa=0

The area burning is in a dense Ponderosa forest known as the Black Forest. I never went around there when living in the Springs from 1992-96. But remember it as a kid (much earlier). It was really pretty. But you look at the maps and can see why all the houses there have been destroyed. The map that shows lost houses (http://gazette.com/article/1502250) is somewhat marred spatially because they used lot centroids rather than pinpointing the structures. When it’s all over, I’m going to suggest to their fire marshals (as well as “fire community” people I know) that when pre-mapping lots they have their GIS people identify the center of the largest visible structure. Why? Because when analyzing the data post-mortem or trying to educate the public it will show which structures were the most vulnerable.

Of course it isn’t as simple as saying,  “Being right under the trees will guarantee destruction. Being in a clearing will not.” Embers travel for miles. And often it’s some structural or landscaping flaw, e.g. dry vegetation right up against the house and under the eaves that causes ignition. An informative public presentation would show time lapse fire and wind travel superimposed on aerial imagery, along with contextual images of damaged structures.

My condolences to anyone who has lost a home in that area.

The “Evil” GEOMETRYCOLLECTION

Recently we encountered a bug in SQL Server’s spatial aggregation algorithms. Sets of polygons that should have been combined strictly into WKT (well-known-text) POLYGON or MULTIPOLYGONs were being stored as GEOMETRYCOLLECTIONs instead. Turns out the algorithm was inserting one or two LINESTRINGs (of two points each) into the collection, forcing it to be stored this way. The bug was supposedly fixed a couple years ago, but I’m going to submit a new bug and sample dataset to Microsoft. Because I’ve double-checked we’re using the latest version of the DLL, etc.

So why do I think the GEOMETRYCOLLECTION is “evil” (or at least don’t think highly of it)? Primarily because it caused our customer to question both the integrity of our software and the use of SQL Server to store AND manipulate spatial data instead of relying on purely an ESRI ArcGIS solution. We were unable to open the SQL Server datasets via ArcCatalog or ArcMap (using the Database Connection in 10.1). On a broader note, I tend to side with ESRI in not supporting more than one data type in a layer, so I’m perplexed that the creators of the WKT standard even thought up such a datatype. I can see where it might make sense to group related geometric objects but believe there are better ways to do that.

I’m tempted (hey, in a post where the word evil is mentioned it seems appropriate word use) to say that this never would have happened in the first place if a GEOMETRYCOLLECTION wasn’t a possibility. But I know that even if the bug gets fixed, Microsoft and others will still adhere to the WKT specification for storing data and ESRI will continue to invalidate coverages that contain multiple datatypes. So we will improve on the hack below to ensure that what gets output in the end is valid by our customer’s definition.

// Well, semi-verified anyway. At least it won't contain any LINESTRINGs.
public SqlGeometry VerifiedPolygonGeometry(SqlGeometry inputGeometry)
{
 SqlGeometry outputGeometry = null;

try
 {
 var goodReservoirGeometriesList = new List<SqlGeometry>();

int numGeometries = (int)inputGeometry.STNumGeometries();

for (int i = 1; i <= numGeometries; i++)
 {

var geometry = inputGeometry.STGeometryN(i);

if (geometry.STGeometryType() != "LineString")
 {
 goodReservoirGeometriesList.Add(inputGeometry.STGeometryN(i));
 }
 }

var agg = new GeometryUnionAggregate();

agg.Init();

foreach (SqlGeometry geometry in goodReservoirGeometriesList)
 {
 agg.Accumulate(geometry);
 }

outputGeometry = agg.Terminate();
 }
 catch (Exception exception)
 {
 Debug.Write(exception.Message);
 }

return outputGeometry;
}

Lovin’ Me Some Cheap Libraries

In my MT State Library days, we worked with double metaphone in conjunction with the GNIS placenames data we had. It was remarkably accurate. I see now that the author of that algorithm has a new version, which is available as a C#, C++ or Java library for $40…quite reasonable for what it does. One of the problems I’ve had with open source is that it sometimes precludes decent code from being propagated because the developer has no way to at least receive some modest compensation. On the other hand, I come across insane prices for code libraries all the time.

What I like is something such as this where it makes sense to purchase because it works and it is so cheap. Recently, I paid $35 for a DLL that does datum shifts (moving lat/lon coordinates from one geodetic model to another). The code was available (after much digging) in Java. But I would have spent 3-4X that amount converting and testing it.

No affiliation with this company: http://www.amorphics.com/

Using SMO to Create a Spatial Index

With the advent of ArcGIS 10.1, you can now direct-connect in a read-only fashion to spatial data on SQL Server. That has worked well and while the read-only aspect means no editing, it also means no locks on the data, which is a problem with File Geodatabases. I’ve found that I can programmatically refresh or rebuild a table in SQL Server and just by panning or zooming have ArcMap redraw it. I’d never be able to leave an FGDB open in ArcMap while trying to do the same.

However, I ran into one problem with really large datasets brought into ArcMap. They couldn’t be exported to FGDB or Shape format without a spatial index on the table in SQL Server. Understandable and really it’s better for performance if a spatial index exists. Since it’s an application that creates the tables and populates them, it fell to the application to also create a spatial index.

I could, of course, have used a parameterized SQL command to do this, but since I’m already using SMO to create a tabular index thought I’d try it for a spatial geometry column. Of course, as I soon realized, there’s a catch…one that almost makes it just as well to use parameterized SQL. But if you like being able to specify things like:

 SpatialIndexType = SpatialIndexType.GeometryAutoGrid 

as needed, then there is some case to be made for mixing SMO with SQL…and mix you have to. Because here’s the catch, you can’t create a spatial index if you don’t already know the bounding box (or envelope) coordinates. The only way I know to do this is via a SQL query to the server first. So here’s a code sample:

var serverConn = new ServerConnection("servername");

serverConn.ConnectTimeout = 180;

// provide appropriate login credentials here

var srv = new Server(serverConn);

Database db = srv.Databases["tablename"];
var tb = db.Tables["tableName"];
try
{
 if (db != null)
 {
 // Perform spatial query to get the bounding box
 var sql = String.Format(@"SELECT
 geometry::EnvelopeAggregate(GEOM).STPointN(1).STX AS MinX,
 geometry::EnvelopeAggregate(GEOM).STPointN(1).STY AS MinY,
 geometry::EnvelopeAggregate(GEOM).STPointN(3).STX AS MaxX,
 geometry::EnvelopeAggregate(GEOM).STPointN(3).STY AS MaxY
FROM {0}", "tableName");

var dataSet = db.ExecuteWithResults(sql);

if ((dataSet != null) && (dataSet.Tables.Count > 0) && (dataSet.Tables[0].Rows.Count > 0))
 {
 var boundingBoxXMin = (Double)dataSet.Tables[0].Rows[0]["MinX"];
 var boundingBoxYMin = (Double)dataSet.Tables[0].Rows[0]["MinY"];
 var boundingBoxXMax = (Double)dataSet.Tables[0].Rows[0]["MaxX"];
 var boundingBoxYMax = (Double)dataSet.Tables[0].Rows[0]["MaxY"];
 //spatial index
 var spatialIndex = new Index(tb, "Spatial" + tableName)
 {
 SpatialIndexType = SpatialIndexType.GeometryAutoGrid,
 BoundingBoxXMax = boundingBoxXMax,
 BoundingBoxXMin = boundingBoxXMin,
 BoundingBoxYMax = boundingBoxYMax,
 BoundingBoxYMin = boundingBoxYMin,
 CellsPerObject = 16,
 PadIndex = false,
 CompactLargeObjects = false,
 DisallowPageLocks = false,
 SortInTempdb = false,
 OnlineIndexOperation = false,
 DisallowRowLocks = false
 };

spatialIndex.IndexedColumns.Add(new IndexedColumn(spatialIndex, "GEOM"));

spatialIndex.Create();

return true;
 }
 }
}
catch (Exception err)
{

Debug.WriteLine(err.message)
}

What’s in a map projection?

If you take a look at the map below of the contiguous U.S. you might ask, “I see that Maine has the easternmost point in the lower 48, but shouldn’t it also have the northernmost point?”

Projections can sometimes be misleading. Well, they’re all misleading in some aspect or another. Take, for example, the Mercator projection, which is used in mapping services like those provided by Bing or Google. (There are sound technical reasons for this BTW). It greatly distorts areas the farther north or south you go. However, in this case, it does make it easy to see that the northernmost point of Maine is south of quite a bit of the U.S.-Canadian border:

So, how far north does Maine go? Approximately 47.46°, which puts it south of that large stretch of the border that runs along the 49th parallel. Incidentally, Sumas, WA – because of a surveying error – extends just slightly past the 49th parallel (49.002389,-122.261111).