Calculating distance in C# (with SqlServer.Types)

WiredPrairie's post on Calculating the distance between two Zip Codes in C# starts with the line

Unfortunately, there are many classes built into the .NET framework, but a zip code distance calculator isn’t included.

which is true -- however with the release of SQL Server 2008 Microsoft has made available the redistributable Microsoft.SqlServer.Types assembly that contains the SqlGeography and SqlGeometry types. Within SQL Server 2008 this enables the new geospatial functions, but these types can also be used in C#. I thought I'd add the SQL Server types to his example, just for fun...

I've modified WiredPrairie's code to use Microsoft.SqlServer.Types.SqlGeography in place of manually calculating Radians and using the Haversine formula... only a few minor changes were required (view source):

0. Reference Microsoft.SqlServer.Types assembly


1. Replace the Distance() method completely
I haven't implemented miles here - just kilometres

The important bits being:

string lineString = String.Format("LINESTRING({0} {1}, {2} {3})"
    , compare.Longitude, compare.Latitude
    , this.Longitude, this.Latitude);
var g = Microsoft.SqlServer.Types.SqlGeography.STGeomFromText(
        new System.Data.SqlTypes.SqlChars(new System.Data.SqlTypes.SqlString(
            lineString
        ))
    , 4326); // IMPORTANT for distance calc
geoDistance = Convert.ToDouble(g.STLength().ToString()) / 1000; // to kilometers

2. Remove the conversion to Radians
The SqlGeography type uses Latitude and Longitude in degrees.


3. Minor text changes (miles to kilometres)


4. and finally, change the input from '25 miles' to '40 kilometres'
Here is a comparison of the results -- there appears to minor differences in the calculation results - I assume due to rounding...


Of course, WiredPrairie's solution would run anywhere (including Silverlight) whereas adding a dependency to Microsoft.SqlServer.Types means you can only use this code where the SqlServer Types redistributable has been installed. It's still interesting to see the additional features provided by SQL Server 2008 - I hope some spatial types (without the SQL dependency) make their way into a future core .NET Framework release (and Silverlight 3.0??).

Great Circles!

No, I'm not over-excited about a shape... Great Circles are a 'geographic concept' (via The Map Room).

A great circle is defined as any circle drawn on a globe (or other sphere) with a center that includes the center of the globe. Thus, a great circle divides the globe into two equal halves.

The post reminded me that some additional work on Geoquery 2008 is long overdue. Until then, it prompted me to blog my 'great circle drawing algorithm' as used in Geoquery 2008 to draw lines like this (a 'straight line' between two points, which looks curved on a Mercator projection but is clearly straight when viewed from above on a globe):

The above example uses SQL Server 2008's STBuffer() geography function - notice how the shaded area looks kidney-shaped on the 'flat' Mercator projection (the curve is more pronounced further north) but on the globe the sides are straight.

You can read the C# code to calculate a Great Circle arc or see it below. Given two System.Windows.Points on a Mercator grid (following the Virtual Earth 'scaling' pattern - notice the hardcoded zoomlevel 2 in the LatLongToPixel method calls) it will calculate a Great Circle between them and plot 10 points along that path (notice the hardcoded 0.1 in the for clause). It returns a System.Windows.Shapes.Polyline ready to draw!


It references Virtual Earth Tile System, and the algorithm itself is based on the work by Ed Williams on Aviation Formulary in particular calculating intermediate points on a great circle.

For those who are really paying attention, you may guess this isn't exactly the code used in Geoquery 2008... for a start 10 divisions was too coarse so Geoquery uses a variable number of points based on the length of the line and its proximity to the poles (where the 'curve' is more pronounced). Geoquery ALSO needs to stop/start drawing lines that cross Longitude 180° (so they wrap nicely on Mercator, and join when projected on a globe)... which is done by calculating the intersection of two great circles because SQL Server 2008 couldn't do it (but that's another story...)

P.S. there are some previous posts on the 'straight line' problem that might be an interesting read...

GeoData Visualization: Geoquery v. SQL Server 2008

Came across this post today - GeoData Visualization vs. Analysis - and thought it would be fun to try in Geoquery.

The dataset (au_sumg.e00) can be found at USGS and imported into SQL Server 2008 using Shape2Sql (by Morten). Among other statistics it contains data on known oil reserves (outside the USA, as of 2000).

The simplest query - SELECT * FROM au_sumg - looks like this in SQL Server 2008 Management Studio and Geoquery 2008 respectively:



Management Studio's Dundas graphic sure looks nice - but a little difficult to interpret in isolation. Geoquery's Shape tab looks worse without even the funky pastel colors.

However, the Map tab in Geoquery starts to make up for the green-ness by giving the shapes some context.


But that's not all - using the same 'scale' as this image and the [Fill] and [Thickness] 'special columns' that Geoquery supports we can write this:

select * ,
CASE
WHEN KWN_OIL < -6238 THEN '#FFF2ECB0'
WHEN KWN_OIL < 15769 THEN '#FFF3CE75'
WHEN KWN_OIL < 37775 THEN '#FFEA944B'
WHEN KWN_OIL < 59781 THEN '#FFB16B3A'
ELSE '#FF864E37'
END AS [Fill]
, 0 AS [Thickness]
FROM au_sumg

to produce this result:


Of course, viewing these results on the globe offers absolutely no advantage; but I'll include an obligatory screenshot anyway...

Geoquery 2008 v0.73 (for SQL Server 2008 RTM)

v0.72 of Geoquery2008 has been recompiled against Microsoft.NET Framework 3.5 SP1 and updated for the RTM of SQL Server 2008. It also has a new home, at

Geoquery2008.com!

It is now aware of the geography coordinate order swap and only requires the Client Profile portion of Microsoft.NET Framework 3.5 Service Pack 1 and the System CLR Types to run.

The next 'major' version - 0.8 - is still a work-in-progress - but will never be designed to compete with the 'built-in' Spatial Results Tab in Management Studio.

Whither Geoquery?

If you were reading this blog earlier in the year, you may have seen these posts on Geoquery for SQL Server 2008.

Since my last post (in March) a few things have happened to stall development:
1) Running. I was injured over summer (that's Dec-Feb down here) and couldn't train - all that extra free time went into developing Geoquery. Now I'm back in training.
2) Holidays. Spent a few weeks in Thailand and Bhutan.
3) Job. It's busier now.
4) I was using the demo version Visual Studio 2008. It's expired. The Express versions just don't cut it for "real" work.
5) I can't help thinking Microsoft is going to release something similar themselves. It was great fun building it up to v0.72, but now real hard work is required to get the functionality "to the next level". That's a lot of effort for "nothing", less fun and more sweat.

I saw today that Google Maps now has a Flash API (and demos), something Yahoo has had for a while. I might get back to playing with silverlightearth.com and deepzoompublisher.com in Silverlight 2.0.

There may yet be a way to tie them all together: converting silverlightearth.com to Silverlight 2.0 will result in a whole pile of C# and Xaml that would then plug right into Geoquery 2008. Will see how things develop... ;-)

Geoquery 2008 v0.72

Another minor release of Geoquery 2008 to finish off some of the 'known bugs' around line and polygon rendering... It can be downloaded here.

First up, Morten has helped fix a number of bugs in Geoquery, and this release is no exception. While trying out his OGC conformance test map I discovered some GEOMETRY idiosyncrasies I still hadn't fixed. If you read his post and download the schema/data, this sql works best in Geoquery.


Morten also pointed out the issue with holes... (notice how the GEOGRAPHY instance is rendered on both the Shape & Map tabs!)


And finally, after a month away from the code (some of it at the beach :) it was relatively quick to solve most of the problems around 'wrapping': with a bit more help from Ed Williams I was able to discard my previous efforts (using the SqlTypes) and solve the problem with some basic math - hoorah!




I expect it will be a bit of a wait for the next version of this code. I have a lot of ideas to incorporate, some of which require further architectural changes. I will also update the Examples to discuss the new features in CTP6. Until then - enjoy Geoquerying...

Back to "work" on Geoquery

After a great holiday at the beach, a rushed implementation of PayPal, some silverlight prototyping and fighting with IL at work, I've finally found some time to work on Geoquery again.

v0.8 is pretty ambitious, so it's likely to be an interim version that fixes a few more bugs, such as polygons that cross +/- 180 degrees Longitude.


and Geography displayed on the Shapes tab


and the dodgy example from my last post


Still a little more to do (what about when the polygon disappears above +/-85 degrees Latitude), but hopefully I'll be back working on street-level maps 'any time now'.



p.s. yeah, I went to the "2008 launch wave". What I don't get is when headlines that look like {data binding} became cool marketing?

Geoquery 2008 v0.71 (oops)

Geoquery 0.7 has been quickly updated to 0.71 to fix the bug from yesterday.

It was while I was describing the problem that I was reminded how confusing it is for users when something "fails silently".

Imagine if the a "File-Save" function didn't report an error when your latest masterpiece couldn't be saved... you'd close the application and expect the data to be available next time - but it would be gone... forever! I strongly discourage rampant use of try{}catch{} in C# because it often (inadvertently) leads to silent failures.

By drawing the polygons on the Map & Sphere, but drawing them incorrectly, Geoquery was giving the impression that it was working when it wasn't. OK, probably not as severe as losing a file, but confusing nevertheless. It seemed important enough to fix and update (even though just a day earlier I'd decided it was "good enough for now").

So now the polygons are drawn correctly, and when they can't be (across +/-180° longitude), it "fails loudly" :)

SELECT geography::STGeomFromText(
'POLYGON((-33.9 151.12,  0 -118.4, 33.8 -118.4,-33.9 151.12))', 4326
), 'Red' as color, '#33ff0000' as Fill

P.S. a couple of other issues were also fixed, thanks Morten!

Geoquery 2008 beta - polygons on sphere NOT!

The releast notes for the latest Geoquery beta release should state in BIG BOLD LETTERS that the POLYGON GEOGRAPHY is not correctly rendered. I knew this but rushed the release out regardless - but this MSDN Forums post reminded me how confusing this can be so I wanted to be sure people use LINESTRING for now.

To illustrate, the poster wants to know why these polygons don't BOTH STIntersect 'Vancouver'...

DECLARE @myPoint geography, @polySmall geography, @polyBIG geography
SET @myPoint = geography::Parse('POINT(49.274138 -123.098562)')
SET @polySmall = geography::Parse('POLYGON((47.0 -124.0, 47.0 -122.0, 50.0 -122.0, 50.0 -124.0, 47.0 -124.0))')
SET @polyBIG = geography::Parse('POLYGON((47.0 -155.0, 47.0  -85.0, 50.0  -85.0, 50.0 -155.0, 47.0 -155.0))')
SELECT @polySmall.STIntersects(@myPoint) AS Intersect_polySmall, @polyBIG.STIntersects(@myPoint) AS Intersect_polyBIG
SELECT @polyBIG , 'Red' as Color, '#44ff0000' as fill
SELECT @polySmall, 'Green' as Color, '#4400ff00' as fill
SELECT @myPoint, 'Blue' as Color, 2 as Thickness

because he's visualizing the shapes looking like this:

This is WRONG WRONG WRONG because Geoquery is not correctly rendering POLYGONs on the curvature of the Earth's surface... so you THINK they should intersect, but STIntersects() correctly returns false.

If we use LINESTRINGs to do the drawing (which Geoquery 2008 v0.7 DOES support) then it's rendered correctly and you can see that they don't overlap:

/* However converting to a LINESTRING which is rendered correctly, it's clear they don't intersect*/
DECLARE @myPoint geography, @polySmall geography, @polyBIG geography
SET @myPoint = geography::Parse('POINT(49.274138 -123.098562)')
SET @polySmall = geography::Parse('LINESTRING(47.0 -124.0, 47.0 -122.0, 50.0 -122.0, 50.0 -124.0, 47.0 -124.0)')
SET @polyBIG = geography::Parse('LINESTRING(47.0 -155.0, 47.0  -85.0, 50.0  -85.0, 50.0 -155.0, 47.0 -155.0)')
SELECT @polySmall, 'Green' as Color, '#4400ff00' as fill
SELECT @polyBIG, 'Red' as Color, '#44ff0000' as fill
SELECT @myPoint, 'Blue' as Color, 2 as Thickness

Even more interesting - the lines that should be parallel ARE, and the lines that shouldn't be (longitude) AREN'T, when drawn on a globe (as they should be).


p.s. if you're wondering why LINESTRING works and POLYGON doesn't; it was (fairly) trivial to handle lines that 'break' over the +/-180 longitude (~international date line), but less so to 'split' POLYGONs into parts to draw independently ... so I left it out (for now). Sorry 'bout that.

UPDATE 31-Jan: this has been fixed and Geoquery 2008 v0.71 is available for download

Geoquery 2008 v0.7 (beta) available for "testing"

The latest version of Geoquery, v0.7, is available for download. Unfortunately it was a bit of a rush to get it out, so I haven't updated all the instructions - hopefully you'll get the hang of it.

Obligatory screenshot:

STIntersects=true; STIntersection=null ??

Not really sure whether this is a 'bug' or me not just 'getting' the single-hemisphere restriction, but sometimes when STIntersects returns true, STIntersection returns null.

I first noticed it when 'drawing' LINESTRINGs across the 'international date line' (longitude=180) which is where the default map edges are -- the "line" from Sydney to Los Angeles necessarily crosses this boundary.

DECLARE @idl geography
DECLARE @idlSouth geography
DECLARE @SYDtoEquator geography
DECLARE @SYDtoLAX geography

SET @idl = geography::STGeomFromText('LINESTRING(85 180, -85 180)',4326);
SET @idlSouth = geography::STGeomFromText('LINESTRING(0 180, -85 180)',4326);
SET @SYDtoEquator = geography::STGeomFromText('LINESTRING(-33.9 151.12, 0 -118.4)', 4326);
SET @SYDtoLAX = geography::STGeomFromText('LINESTRING(-33.9 151.12, 33.8 -118.4)', 4326);

SELECT @SYDtoEquator.STIntersects(@idl) as [Intersects]
     , @SYDtoEquator.STIntersection(@idl) as [Intersection]
     , @SYDtoEquator, 'Sydney to Equator' as [Desc]
     , 'why is Intersection null?' as [Question]
SELECT @SYDtoLAX.STIntersects(@idl) as [Intersects]
     , @SYDtoLAX.STIntersection(@idl) as [Intersection]
     , @SYDtoLAX , 'Sydney to LAX' as [Desc]
/* shortening the 'idl' line 'fixes' it */
SELECT @SYDtoEquator.STIntersects(@idlSouth) as [Intersects]
     , @SYDtoEquator.STIntersection(@idlSouth ) as [Intersection]
     , @SYDtoEquator, 'Sydney to Equator' as [Desc]
SELECT @SYDtoLAX.STIntersects(@idlSouth) as [Intersects]
     , @SYDtoLAX.STIntersection(@idlSouth) as [Intersection]
     , @SYDtoLAX, 'Sydney to LAX' as [Desc]

This is the result "set"


And this is how the lines "should look"


As the line endpoint changes, so does the STIntersection result (although STIntersects returns true every time)

@SYDtoEquator = geography::STGeomFromText
    ('LINESTRING(-33.9 151.12, 10 -118.4)', 4326); -- NULL
@SYDtoEquator = geography::STGeomFromText
    ('LINESTRING(-33.9 151.12, 20 -118.4)', 4326); -- NULL
@SYDtoEquator = geography::STGeomFromText
    ('LINESTRING(-33.9 151.12, 30 -118.4)', 4326); -- GOOD
@SYDtoEquator = geography::STGeomFromText
    ('LINESTRING(-33.9 151.12, 40 -118.4)', 4326); -- GOOD
@SYDtoEquator = geography::STGeomFromText
    ('LINESTRING(-33.9 151.12, 50 -118.4)', 4326); -- NULL

So... any thoughts why some of these 'succeed' while others fail?

Spatial queries: not just for Earth

No idea what real applications this might have, but it was kinda fun:

/* http://www.fi.edu/pieces/schutte/landchart.html */
select         -- Apollo landings
  geography::STGeomFromText('POINT(0.67   23.49)', 4326)--1
, geography::STGeomFromText('POINT(-2.94 -23.45)', 4326)--2
, geography::STGeomFromText('POINT(-3.67 -17.46)', 4326)--14
, geography::STGeomFromText('POINT(26.11  3.66)', 4326)--15
, geography::STGeomFromText('POINT(-8.6  15.31)', 4326)--16
, geography::STGeomFromText('POINT(20.17 30.8)', 4326)--17
, 15.0 AS [Thickness], '#770000ff' AS [Color]

The map data in Geoquery is in a configurable Xml file - MapSources.xml so it's relatively simple to include different tile servers for Earth images (as well as the Moon & Mars)!

geography:STBuffer() IS 'straight'

It can be difficult to visualize "geographic" straight lines (intersections, overlaps, etc) on a 'flat' projection, so here's a set of STBuffers over the LAX-JFK line... it doesn't "look" straight until you see it on a globe:


Using the 'special columns' feature of Geoquery (and a different map...)

DECLARE @g GEOGRAPHY
SET @g = geography::STGeomFromText('LINESTRING(33 -118, 40 -73)', 4326); 
SELECT @g, 'Blue' AS [COLOR], 2 AS [Thickness]
UNION ALL SELECT @g.STBuffer(450000),  'Blue' as [COLOR], 1 AS [Thickness]
UNION ALL SELECT @g.STBuffer(1500000), 'Green' as [COLOR], 1 AS [Thickness]
UNION ALL SELECT @g.STBuffer(3000000), 'Purple' as [COLOR], 1 AS [Thickness]
UNION ALL SELECT @g.STBuffer(4000000), 'Orange' as [COLOR], 1 AS [Thickness]

Geoquery: when a straight line isn't straight...

I was already aware that straight lines on a sphere (aren't really straight), but had left the correct handling of such matters out of Geoquery 0.6.

Unfortunately, this made it look rather incomplete; simply drawing a straight line over a projected map...

...isn't quite right

There are quite a few references for how to calculate the 'correct line' - I liked Ed Williams Aviation Formulary - and finally got around to giving it a try.

Here is how the query (and result) should look; with a couple of extra STBuffers thrown in...
DECLARE @g geography;
-- LAX 33°57 -118°24
-- JFK 40°38 -73°47
SET @g = geography::STGeomFromText(
'LINESTRING(33 -118, 40 -73)', 4326);
SELECT @g, null,null,null, 'Blue' as [COLOR]
UNION ALL
SELECT @g.STBuffer(450000), @g.STBuffer(1500000), @g.STBuffer(3000000), @g.STBuffer(4000000), 'Orange' as [COLOR]


The interesting thing about this canvas is the 'straight line' was the only element that I needed to programmatically 'project' -- Katmai's spatial support in the GEOGRAPHY::STBuffer() result is already a POLYGON whose points can be applied directly to the map-projection (ie. it isn't just described as two semicircular shapes joined by two parallel straight lines).

Geoquery can also correctly show buffers on points [aren't] circular on a projected map; the query from that link is reproduced below

DECLARE @g geography;
-- Copenhagen 55°37 12°39
SET @g = geography::STGeomFromText('POINT(55.55 12.66)', 4326);
SELECT @g, 'Blue' as [COLOR]
UNION ALL
SELECT @g.STBuffer(4000000), 'Purple' as [COLOR]



NOTE: Geoquery 2008 0.7 (which contains the enhancements described above) is NOT YET AVAILABLE for download. When it is available, it will also include the complete set of GEOGRAPHY Examples in addition to the GEOMETRY ones already available.

Geoquery 2008 v0.6 (beta)

Geoquery 2008 now supports the GEOGRAPHY datatype as well - displaying the results on a Virtual Earth map rather than a white shape canvas.


v0.6 also adds an Export... function so you can save spatial query results graphically, for emailing/reports/etc.



Here are a couple of Exported resultsets from the MSDN Sample queries...

Geoquery: Geography MkII

Just displaying a series of points doesn't look like much



but it can be the foundation of some neat spatial reporting - this query (attempts to) show the density of votes by postcode.

select p.postcode, p.location,
CASE
WHEN p.Postcode < 2000 THEN '#44ff0000'
WHEN p.Postcode < 3000 THEN '#4400ff00'
WHEN p.Postcode < 4000 THEN '#440000ff'
WHEN p.Postcode < 5000 THEN '#44880088'
WHEN p.Postcode < 7000 THEN '#44FFFF00'
END AS [Color]
, CASE WHEN Totalvotes < 1000 THEN 1
WHEN Totalvotes < 2000 THEN 2
WHEN Totalvotes < 3000 THEN 3
WHEN Totalvotes < 4000 THEN 4
WHEN Totalvotes < 5000 THEN 5
ELSE 6 END AS [Thickness]
from Postcodes P
INNER JOIN PollingPlaces pp ON pp.Postcode = p.Postcode
INNER JOIN PollingPlaceVotes ppv ON ppv.PollingPlaceID = pp.PollingPlaceID



Note the status bar - that's 56,951 rows!! Try loading those points into Virtual Earth! :-)

Geoquery is learning maps

Turns out basic 'map' support was pretty easy (although it's not quite done and available for download just yet...)

This article came in handy - Virtual Earth Tile System - I'd previously used it as a reference for SilverlightEarth.com, although I had to convert it to jscript for that.

SQL Server 2008: Geoquery available for download

DOWNLOAD Geoquery 2008 beta (96Kb) here

SQL Server 2008 Geography: STExteriorRing, STInteriorRingN

MSDN describes STExteriorRing() with this query
DECLARE @g geometry;
SET @g = geometry::STGeomFromText('POLYGON((0 0, 3 0, 3 3, 0 3, 0 0),(2 2, 2 1, 1 1, 1 2, 2 2))', 0);
SELECT @g.STExteriorRing().ToString();
which returns
LINESTRING (0 0, 3 0, 3 3, 0 0, 0 0)

Here it is visually -- the Green line is the original geometry (one square inside another), the thicker Blue line shows that STExteriorRing and the Orange line shows STInteriorRing(1):
DECLARE @g geometry;
SET @g = geometry::STGeomFromText('POLYGON((0 0, 3 0, 3 3, 0 3, 0 0),(2 2, 2 1, 1 1, 1 2, 2 2))', 0);
SELECT @g.STExteriorRing(), 0.3 as thickness, 'Blue' as color
UNION ALL
SELECT @g.STInteriorRingN(1), 0.2 as thickness, 'Orange' as color
UNION ALL
SELECT @g, 0.1 as thickness, 'Green' as color


Interestingly, this query exposed a bug in Geoquery where the Exterior and Interior rings were being 'joined' (ie there was no "pen up" occuring as the shapes were generated) -- causing the dodgy looking lines on the world map in the last post. Here's the "fixed" map drawing...