…or why spherical coordinate systems are not flat!

One of the most common ways the round world is displayed on a map is using the simplest projection we have:

x = longitude
y = latitude

The name of this projection is “Plate Carree”, and is widely used because it is so simple. However we often seem to forget that we are talking about a projection. Therefore the spatial reference for this projection is very often (mis)referenced as a spherical coordinate system like the following for EPSG:4326:

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],
PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]

GEOGCS denotes “Geographical Coordinate System”, meaning a spherical coordinate system using angles for latitude and longitudes, not X’s and Y’s. However, the moment we project this map onto a flat screen or piece of paper, we projected the coordinate system (using the simple formula above), so how can the spatial reference of our map be the above one?

I found that in ArcMap you can manually define a projected coordinate system that renders the same map, but correctly identifies the coordinate system as being projected. Below is the Well-known-text representation of this projection string (note how PROJCS denotes projected coordinate system, and how EPSG:4326 is defined inside it):

PROJCS["World_Plate_Carree_Degrees",GEOGCS["GCS_WGS_1984",
DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],
PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],
PROJECTION["Plate_Carree"],PARAMETER["False_Easting",0],
PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",0],
UNIT["Degree",111319.49079327357264771338267056]]

As far as I know we don’t have a spatial reference ID for the above (EPSG:54001 comes close but uses a different linear unit), so we often incorrectly start using EPSG:4326 to describe a projection. This has been done for so long, that this is now common (mal)practice. Even OGC’s WMS specification allows you request flat maps in a spherical coordinate system, even though this doesn’t really make any sense.

So why is this a problem? Well until recently it hasn’t really been a problem, mostly because it was rare that the spherical coordinate were correctly handled, and most applications assumed that the world was flat. However, Microsoft SQL Server 2008 can correctly handle spherical coordinates, and suddenly this becomes a major problem.

Below is a query on all the counties that follow the red line defined by two points. The blue polygons is the result returned. The map claims that its spatial reference is EPSG:4326. Because of the curvature Earth, a straight line between the two endpoints looks like a curve on the projected map, so the results returned seems like they don’t match with the line, but they are in fact the counties on the line between the two endpoints if you think in a spherical coordinate system.

The real problem here is not that the line is not really a curve (I drew the line like that because I want the features along that latitude). The problem is that the line didn’t know that it was projected, because the map it was displayed on was incorrectly set to EPSG:4326. However had the application known that this was not a spherical coordinate system, but a Plate Carree projection, it would have known that when converting from the flat coordinate system to the spherical coordinate system, it should ensure that my line follows the latitude. But because the input line already claims to be in spherical units, the application can not know that it needs to do anything to it (aside the fact that the application of course knows that this came from a flat screen map, but the business logic is of course separate from the UI).

This brings me to another even worse issue: I’d like to make the claim that almost all of the data you have that claims to be in EPSG:4326 has never been EPSG:4326 ! (point data excluded).

Example: Most of the northern US border roughly follows the 49’th latitude. Let’s for sake of argument define the whole northern border as two points, from coast to coast, 49 degrees north. If I then think that my data is in EPSG:4326 (or any other geographic coordinate system), the border will NOT be following the 49th latitude, but go along a great circle that cuts into Canada.

So let us thank Microsoft for creating a database that handles spherical coordinates correctly, and for giving us major headaches when trying to handle these things correctly in a clean clear way :-). Of course there wouldn’t be any headaches if we never mixed spherical and flat coordinate systems in the first place.

I attended the Microsoft PDC in Los Angeles Convention Center this week, and here's a recap of various statements, quotes and sessions that I found interesting.

WCF Preconference notes

It looks like a class, it feels like a class, it smells like a class, it works like a class – but it’s a service!

This was a full day seminar on Windows Communication Foundation. It was very cool to see how ordinary looking classes could be turned into WCF services that can be exposed through a multitude of protocols using just two class/method attributes. One of the interesting statements were how .NET had replaced COM 7 years after its introduction, and now 8 years later, WCF is replacing .NET (well at least the part that replaced COM). Or to put it in the presenters word: ".NET as you know it is dead!"

WCF is host agnostic. No difference between hosting on IIS, in-proc, self-hosting, Windows Activation Service etc.

Possible transport protocols: HTTP, HTTPS, TCP, P2P, IPC, MSMQ.
Possible message formats: Plain text, Binary, MTOM.
Security options: None, Transport security, Message security, Authentication and authorizing callers.

WCF is the first platform that provides Interoperability, Productivity/Quality AND Extensibility all at once.

Favorite sessions

I must admit that many of the sessions were a little dissappointing. They were often marked as advanced or expert level, but rarely lived up to that. There were also way too many Twitter demos. But having said that, here are the ones I attended that didn't dissappoint. You can watch them all online by clicking the links next to them or see the whole list of sessions here.

• TL16: "The Future of C#" by Anders Hejlsberg. My Danish hero and the brain behind C# talks about where C# is headed with upcoming v4 and 5 of C#. View
• TL49: "Microsoft .NET Framework: Overview and Applications for Babies" by Scott Hanselman. Using a a wide varity of .NET core technologies, Scott takes his BabySmash application and ports it to WPF, Silverlight, Mobile and Surface, in his usual entertaining way. Great demonstration of the power of .NET. View
• PC06: "Deep Dive: Building an Optimized, Graphics-Intensive Application in Microsoft Silverlight" by Seema Ramchandani. Seema goes through all the gory internal details of how Silverlight works and uses this knowledge to present a great set of tips and tricks to making the performance of your Silverlight applications scream. A must-see if you are doing serious work with Silverlight. View
• TL26: "Parallel Programming for Managed Developers with the Next Version of Microsoft Visual Studio" by Daniel Moth. Processors are not really getting much faster anymore, but instead we get more and more cores to work with, but this also requires us to start changing the way we make software. The Parallel Framework that comes in .NET 4 makes this task really easy. All I can say is that this framework rocks! View
• BB24: "SQL Server 2008: Deep Dive into Spatial Data" by Isaac Kunen. Isaac's deep dive on the spatial support in SQL Server (thanks for the plug Isaac :-). View
• PC29: "Microsoft Silverlight 2: Control Model" by Karen Corby. Good information on how to build reusable, skinnable controls for Silverlight. View
• PC32: "ASP.NET AJAX Futures" by Bertrand Le Roy. Upcoming features for ASP.NET AJAX. View
• TL46: "Microsoft Visual C# IDE: Tips and Tricks" by Dustin Campbell. LOTS of great shortcuts and features in Visual Studio that you never knew was there and you wonder how you could ever live without. MUST SEE if you want to be more efficient when coding C# in Visual Studio. View

My favorite comment was Juval Lowy's reaction when he the first time heard that Microsoft was working on the successor to VB6:

Giving VB developers access to a multithreaded environment (VB.NET) is like giving razorblades to babies.

It’s not that C++ developers are better off with C#, but they are more used to seeing blood.

Microsoft is already planning a new PDC next year November 17-20, 2009 (unfortunately same place).

According to a Microsoft announcement the Express Edition of SQL Server 2008 won't be available until the end of August.

This however turns out not to be entirely true. Itching to try out the new version, I ended up downloading the evaluation version. It turns out that during the installation you get to choose between installing the evaluation and the express versions.

Unfortunately I can't get the management studio to install though (would have loved to try out the new query visualizer). I don't even see it in the install options, and executing the install directly from the DVD throws an error. At least my tool seems to continue to work with the RTM version.

It looks like Ed Katibah (Spatial Program Manager of SQL Server) finally started blogging.

His first blogpost covers some interesting changes in the latest February CTP.

There are some new methods added to the list:

• Geometry::Filter
• Geography::Filter
• Geography::EnvelopeCenter
• Geography::EnvelopeAngle
• Geography::STDistance is no longer limited to having one of the inputs of type Point.
• Geography::Reduce

Unfortunately I can't try these new features today, because a bug prevents you from using SQL Server 2008 on a Leap Day and the doc on MSDN doesn't cover these new methods yet, so there's not telling how these methods really work.

Ed's and Isaac Kunens blogs are a must-read if you plan on doing some serious SQL Server Spatial work.

The Simple Features Specification v1.1.0 has a conformance test based on "Joe's Blue Lake vicinity map", which roughly looks like this:

To create the map in Microsoft SQL Server 2008, first create the tables:

-- Lakes
CREATE TABLE lakes (fid INTEGER NOT NULL PRIMARY KEY,name VARCHAR(64),shore Geometry);
-- Road Segments
CREATE TABLE road_segments (fid INTEGER NOT NULL PRIMARY KEY,name VARCHAR(64),aliases VARCHAR(64),num_lanes INTEGER,centerline Geometry);
-- Divided Routes
CREATE TABLE divided_routes (fid INTEGER NOT NULL PRIMARY KEY,name VARCHAR(64),num_lanes INTEGER,centerlines Geometry);
-- Forests
CREATE TABLE forests (fid INTEGER NOT NULL PRIMARY KEY,name VARCHAR(64),boundary Geometry);
-- Bridges
CREATE TABLE bridges (fid INTEGER NOT NULL PRIMARY KEY,name VARCHAR(64),position Geometry);
-- Streams
CREATE TABLE streams (fid INTEGER NOT NULL PRIMARY KEY,name VARCHAR(64),centerline Geometry);
-- Buildings
CREATE TABLE buildings (fid INTEGER NOT NULL PRIMARY KEY,address VARCHAR(64),position Geometry,footprint Geometry);
-- Ponds
CREATE TABLE ponds (fid INTEGER NOT NULL PRIMARY KEY,name VARCHAR(64),type VARCHAR(64),shores Geometry);
-- Named Places
CREATE TABLE named_places (fid INTEGER NOT NULL PRIMARY KEY,name VARCHAR(64),boundary Geometry);
-- Map Neatline
CREATE TABLE map_neatlines (fid INTEGER NOT NULL PRIMARY KEY,neatline Geometry);

Next step is to insert the rows:

-- Lakes
INSERT INTO lakes VALUES (101, 'BLUE LAKE',Geometry::STPolyFromText('POLYGON((52 18,66 23,73 9,48 6,52 18),(59 18,67 18,67 13,59 13,59 18))', 101));
-- Road segments
INSERT INTO road_segments VALUES(102, 'Route 5', NULL, 2,Geometry::STLineFromText('LINESTRING( 0 18, 10 21, 16 23, 28 26, 44 31 )' ,101));
INSERT INTO road_segments VALUES(103, 'Route 5', 'Main Street', 4,Geometry::STLineFromText('LINESTRING( 44 31, 56 34, 70 38 )' ,101));
INSERT INTO road_segments VALUES(104, 'Route 5', NULL, 2,Geometry::STLineFromText('LINESTRING( 70 38, 72 48 )' ,101));
INSERT INTO road_segments VALUES(105, 'Main Street', NULL, 4,Geometry::STLineFromText('LINESTRING( 70 38, 84 42 )' ,101));
INSERT INTO road_segments VALUES(106, 'Dirt Road by Green Forest', NULL, 1,Geometry::STLineFromText('LINESTRING( 28 26, 28 0 )',101));
-- DividedRoutes
INSERT INTO divided_routes VALUES(119, 'Route 75', 4,Geometry::STMLineFromText('MULTILINESTRING((10 48,10 21,10 0),(16 0,16 23,16 48))', 101));
-- Forests
INSERT INTO forests VALUES(109, 'Green Forest',Geometry::STMPolyFromText('MULTIPOLYGON(((28 26,28 0,84 0,84 42,28 26),(52 18,66 23,73 9,48 6,52 18)),((59 18,67 18,67 13,59 13,59 18)))', 101));
-- Bridges
INSERT INTO bridges VALUES(110, 'Cam Bridge', Geometry::STPointFromText('POINT( 44 31 )', 101));
-- Streams
INSERT INTO streams VALUES(111, 'Cam Stream',Geometry::STLineFromText('LINESTRING( 38 48, 44 41, 41 36, 44 31, 52 18 )', 101));
INSERT INTO streams VALUES(112, NULL,Geometry::STLineFromText('LINESTRING( 76 0, 78 4, 73 9 )', 101));
-- Buildings
INSERT INTO buildings VALUES(113, '123 Main Street',Geometry::STPointFromText('POINT( 52 30 )', 101),Geometry::STPolyFromText('POLYGON( ( 50 31, 54 31, 54 29, 50 29, 50 31) )', 101));
INSERT INTO buildings VALUES(114, '215 Main Street',Geometry::STPointFromText('POINT( 64 33 )', 101),Geometry::STPolyFromText('POLYGON( ( 66 34, 62 34, 62 32, 66 32, 66 34) )', 101));
-- Ponds
INSERT INTO ponds VALUES(120, NULL, 'Stock Pond',Geometry::STMPolyFromText('MULTIPOLYGON( ( ( 24 44, 22 42, 24 40, 24 44) ),( ( 26 44, 26 40, 28 42, 26 44) ) )', 101));
-- Named Places
INSERT INTO named_places VALUES(117, 'Ashton',Geometry::STPolyFromText('POLYGON( ( 62 48, 84 48, 84 30, 56 30, 56 34, 62 48) )', 101));
INSERT INTO named_places VALUES(118, 'Goose Island',Geometry::STPolyFromText('POLYGON( ( 67 13, 67 18, 59 18, 59 13, 67 13) )', 101));
-- Map Neatlines
INSERT INTO map_neatlines VALUES(115,Geometry::STPolyFromText('POLYGON( ( 0 0, 0 48, 84 48, 84 0, 0 0 ) )', 101)); 

And voilá... you can now for instance run all the conformance tests, but the map is actually also a good base-map for learning the various spatial operations.

If you use my SQL Spatial Query Tool, here's a query you can use to visualize the map:

SELECT neatline, 'Background' as name, 'map_neatlines' as layer, 'White' as FillColor, 'Transparent' as LineColor, 1 as LineThickness FROM map_neatlines UNION ALL
SELECT shores, name, 'ponds' as layer, 'LightBlue' as FillColor, 'Blue' as LineColor, 1 as LineThickness FROM ponds UNION ALL
SELECT boundary, name, 'forests' as layer, 'Green' as FillColor, '#ff005500' as LineColor, 4 as LineThickness FROM Forests UNION ALL
SELECT shore, name, 'lakes' as layer, 'Blue' as FillColor, 'Transparent' as LineColor, 2 as LineThickness FROM lakes UNION ALL
SELECT boundary, name, 'named_places' as layer, '#00ffff99' as FillColor, 'Brown' as LineColor, 2 as LineThickness FROM named_places UNION ALL
SELECT centerline, name, 'streams_outline' as layer, 'Blue' as FillColor, 'Blue' as LineColor, 5 as LineThickness FROM streams UNION ALL
SELECT centerline, name, 'streams_fill' as layer, 'LightBlue' as FillColor, 'LightBlue' as LineColor, 3 as LineThickness FROM streams UNION ALL
SELECT centerline, name, 'road_segments' as layer, 'LightGreen' as FillColor, 'Black' as LineColor, num_lanes*2 as LineThickness FROM road_segments UNION ALL
SELECT centerlines, name, 'divided_routes' as layer, 'LightGreen' as FillColor, '#aaff0000' as LineColor, num_lanes*2 as LineThickness FROM divided_routes UNION ALL
SELECT footprint, address, 'buildings_footprint' as layer, 'Black' as FillColor, 'Transparent' as LineColor, 1 as LineThickness FROM buildings UNION ALL
SELECT position.STBuffer(0.5), address, 'buildings_position' as layer, 'Red' as FillColor, 'Transparent' as LineColor, 1 as LineThickness FROM buildings UNION ALL
SELECT position.STBuffer(1), name, 'bridges' as layer, '#550000ff' as FillColor, 'Black' as LineColor, 1 as LineThickness FROM bridges

Which will generate a map like this: 

Download SQL queries (1.71 kb)

There was a question in the Sql Server Spatial/Katmai forums about why a straight line on a sphere didn't "look straight".

Katmai both have planar and spherical geometry types. The spherical type makes it easier to handle things like straight lines, distances and areas over large distances, and depending on which type you choose, you can get very different results.

Consider traveling on a straight line from Vancouver to Frankfurt. They are both located roughly at latitude 50°N. So the line would follow latitude 50°N right? Well that's true for a planar coordinate system:

So according to The Flat Earth Society (who actually believes the Earth is flat) the green line above shows the shortest distance between the two cities. But what we are actually seeing is an artifact of working with a planar coordinate system depicting a world that is not planar.

Using a spheric coordinate system, we get a different result. When looking at this line in a projected coordinate system, like the Mercator projection I've used here, the line doesn't look straight (left image below). However if we look at this in a 3D view (right image), it becomes apparent that the "curved line" is actually the shortest line between the two cities.

The difference between using the two geometry types gets really apparent when you try to find the shortest distance from New York to the line in the two coordinate systems. In the planar case, Quebec is approximately the point where the line is the closest to New York. In the spherical case, the line never gets closer to New York than from Vancouver where it started out.

Another way that I like to think of why following a latitude is not the shortest route, is imagining two guys standing a few feet from the North Pole with the pole between them. Which way is the shortest path between them: Walking straight over the pole or following the latitude at 89.999°N around it? Now imagine how this line would look on a map.

We can extend this further and have the two men go further and further away from each other.  If one man gets a little further south than the other, the line wouldn't pass the pole, but go by just next to the it - it still wouldn't follow the latitude. Since the North Pole is just a point we humans "made up", we can abstract it and think of a pole in the ground anywhere on Earth. The bottom line is that the North Pole and the latitudes are all imaginary, and doesn't really exist, and definitely doesn't help describing the shortest distance.

Only at Equator would it be a good idea to follow the latitude, which is because this is the only latitude that is a great circle, and it so happens to be that any straight line between two points will describe parts of a great circle. Longitudes however are all great circles, and the line between the two men describe parts of a longitude. The reason we think that the horizontal line between Vancouver and Germany is the shortest path, is simply because we have to project the round Earth onto a flat piece of paper, and this causes all sorts of things to get distorted.

It looks like Microsoft listened to the community and plans on switching latitude and longitude ordering for WKT and WKB:
http://blogs.msdn.com/isaac/archive/2007/12/27/latitude-longitude-ordering.aspx

This change will only be for WKT and WKB. That does make some sense, since the ordering is related to the SRID which is not included in the WKT and WKB, thus making it impossible to decipher without any prior knowledge. It would still be nice if the two geometry models shared a more common interface for anything that is remotely similar, but this is a good start.

I’m actually starting to like working with the geography type (to begin with I hated the idea of having two different although similar geometry models). It’s fun to do buffers on points and see that they don’t look circular on a projected map, because they are, well… projected. Showing the same circles on a 3D globe reveals that they are in fact perfectly circular. It makes distance queries so much easier than for instance using this approach.

One thing that still bugs me about the geography datatype though… why is the linear unit type that is used for calculating area, length, buffers etc. declared as part of the SRID (which in this case is always geographic) ? This requires one to always first have to look up the linear unit in the spatial_ref_sys table for each returned row (since homogeneous SRIDs aren’t enforced on the table) and divide by the factor returned to ensure that we know the unit we get back. In my mind a linear unit should never be related to a geographic coordinate system. After all, that’s why we only define linear units in projected and geocentric coordinate systems where they makes sense.

It would be nice if the methods either took a unit type as input (although that would break with the common interface I just wished for) or simply just always outputs the results in SI units (then it would just be a matter of multiplying with a constant well-known conversion factor later if you want to display it in for instance miles, thus separating code from presentation). Luckily most of Katmai’s SR definitions use meters as default, but you never know if someone chose to change the conversion factors in the system table so you will always have to check first.

I've been playing with the spatial support of SQL Server 2008 today, but quickly hit two obstacles:

1. How do I get data into the database.
2. How do I visualize a query (apart from the Well-Known Text output).

Surprisingly there aren't really any good (free) tools out there yet to upload for instance shapefiles to Katmai. I pulled out some old code I had lying around and threw a small tool together that can upload a shapefile to the database either as geometry or geography. Katmai is pretty picky when it comes to the geography datatype, so make sure your shapefiles contain valid geometry and has vertices in the valid domain, or use the geometry format. You can download the tool below.

Next obstacle was visualizing the data. Playing with the spatial queries and not beeing able to see the results on a map quickly got boring. I haven't had a chance to play much with WPF, so I thought this was a good occation to play with it. Based on this MSDN article I managed to fairly quickly make a little tool that takes a SQL query as input and outputs the results in a WPF control. At the moment it only works with polylines and polygons - Points are simply ignored (tip: use STBuffer() on points to see them). It will look for the first column that contains a SqlGeometry or SqlGeography datatype and put that on the map. The remaining columns are shown in the box to the right when you hover over the shapes with the mouse.

The two tools are not perfect in any way (and there's no bug tracking support here :-)), but they provide a quick way to play with the spatial queries. Both tools requires Windows Authentification to access the sql server and runs on .NET 3.5. Furthermore the query tool uses the Microsoft.SqlServer.Types.dll that comes with installing Katmai (to parse the results), so for now you'll need to run it from a computer that has Katmai installed.

Now here's an interesting use of the spatial support in Microsofts SQL Server 2008: Battleship.

Using points, linestrings and hit-testing, Craig developed a small battleship game you can play by executing some stored procedures:

EXEX NewGame 'Craig'
EXEC Cheat 13,4
EXEC Shoot 13,4,1,1 -- Miss
EXEC Shoot 13,4,3,1 -- Hit
EXEC Shoot 13,4,3,2 -- Sunk

I'm not sure this was what Microsoft had in mind, but definitely nice work Craig!