Creating OGC conformance test map in SQL Server 2008

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)

Straight lines on a sphere

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.

The community works

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.

400km buffer around Copenhagen

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.

Shapefile to SqlServer 2008 to WPF

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.

Read more and download SqlSpatialTools

BattleShip - Now for SqlServer 2008

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!

Linear maps

The Map Room blogs about San Francisco's new public transport maps that will employ more "straight lines", also covered more in depth here. I find these sort of map abstractions really interesting, so I thought I would chip in with a couple of similar "straight-lines" public transport maps from Copenhagen. It must take a good degree of ingenuity to weigh accuracy with simplicity.

The first one shows the most used map for the network using straight lines as much as possible. You can still clearly see the "five-finger plan" that had been used for city planning for many years by expanding out in "fingers" and allow for greener less dense areas in between. Click the image for a larger view:

The second shows a much more dense map (usually used on a larger poster) which also shows some of the major bus routes. This map is much more geographically correct with a slight grey outline of the land behind it. It still uses straight lines to simplify the map (click for larger view):

There's an even simpler version, where the map is far from geographically correct, but instead shows the network and connections as a set of mostly horizontal lines. This map is mostly used inside the train, so at this point you (hopefully) already know where you are headed and don't need any geographic resemblence. Therefore this "map" only shows which stations the train lines stops at/skips and where they connect to each other,and it fits nicely over the door :-) (unfortunately this was the best picture I could find):

Comments working again

Thanks to Al for pointing out that my comments section was failing with JavaScript errors. Apparently the blowery http compressor that this site is using was messing up the embedded ASP.NET JavaScripts so they return empty script files. Fortunately a workaround exists for this.

So if you meant to comment about one of the screw-ups in my recent blogposts, here's your chance.

Huge wind turbines in Virtual Earth

Since today is Blog Action Day and the topic is the environment, I thought I would share this amazing birds eye view of one of the 2mW windmills in the offshore windfarm "Middelgrunden" just off the coast of Copenhagen that was published last week. I've knew these windmills were BIG, but looking from the peer it's hard to tell how big they really are. The boat next to it really gives a sense of its size. I can't even begin to imagine how big the 5mW mills they make these days are... 


There are 20 of these mills, all together producing 3% of Copenhagens total electricity consumption. Around 20% of Denmarks electricity comes from windmills.

Warning... bad joke coming up: Now you know why I moved to Sunny California - The weather really "blows" in Denmark :-)

 

New Virtual Earth release

I just visited http://maps.live.com and it looks like they finally put the new release out there. A revamped design and Birds Eye View imagery inside the 3D viewer. You can now also create animations using your collections.

Nyhavn, Copenhagen - A place I missed a lot this summer

Another new feature is the new "Create model" button in the collection editor when you are in 3D mode. Using it will prompt you to download "3dvia" enabling you to model your own buildings and insert them into VE. It looks somewhat similar to SketchUp used for the same thing in Google Earth. A neat feature is that you get the map as a background to assist you in your modelling, and comes with a bunch of default textures as well:

Accurate distance calculations in UTM projections

Recently a friend of mine asked me what projection would be the best to create point-buffer circles. All map projections adds distortion so you rarely end up with a circle. If you’re lucky, you’ll have an ellipse but often they get even more complex than this. In many cases a “normal” circle is sufficient for accuracy, but still it got me thinking how I would do this 100% accurate. There are many approaches to this, fx. creating a new projection for each point, do the math in spherical coordinate systems etc, but none of them are completely trivial.

The UTM projection is one of the most commonly used projections. It’s fairly accurate to measure from point to point within small distances and close to the meridian, so this post is based on that. Most of the math here can also be applies to Mercator, where the scale reductions are applied to Y instead of X.

The Mercator projection is created by projecting the globe onto a cylinder whose center follows the Earth’s rotation axis. The Transverse Mercator is basically a "tilted" cylinder parallel to equator. Along the line that the cylinder "touches" the surface the distortion is 0, and increasing away from it in the east/west direction. There is no distortion in the North/South direction. The Universal Transverse Mercator is slightly smaller than the globe, so it will intersect the surface two places. This also ensures that the distortion between and around these two lines are minimal. Along the meridian (the center of the two lines) the distortion for the UTM projection is 0.9996, meaning that this is the amount you have to reduce an infinitely small line in the east/west direction (For Mercator this is normally 1 along Equator). As you move away from the meridian, the scale factor increases up toward infinity. Normally UTM projections are used close to the meridian, and every time you get too far east/west, you would switch to use a new UTM projection. That’s why The Earth is divided into 60 UTM ‘zones’, one for each 6 degrees. Often for practical reasons you might want to use a UTM zone even though it’s outside its defined usage area. Here you have to be specifically careful when measuring distances. For example, Denmark is covered by zone 32 and 33, but it’s common to only use zone 32 which means that distance calculations at the east end can be very inaccurate without applying a scale reduction.

Warning: The following contains a lot of math. You can skip to the bottom and download a C# class that does all the calculations for you.

The scale factor SR(x) for any given point at distance 'x' from the meridian is given by:

Since this is the scale reduction at a point, we need to sum up all of these values along a line to find the correct distance, thus the east/west distance from e1 to e0 (expressed in distance from the meridian) becomes an integral expression:

To this you need to add the north/south distance using the well-known Euclidean distance formula.

For atypical UTM using WGS84 ellipsoid and a 500000m false easting it becomes:

If you haven’t refreshed your integral math lately, or just want to write this in pure code, one can also write the expression as:

Where m is the scalefactor at the meridian, r is the earth radius and e1 and e0 have the false easting subtracted.

Often we can do with an approximate value by using the scale reduction between the start and end point. This should only be used on distances <1km in easting because it gets very inaccurate for greater distances.

Simpson’s Rule

One can also use Simpson’s rule for calculating the integral solution. My tests showed it’s a very accurate approximation for this type of integral usually giving millimeter accuracy. It’s defined as:

I performed a few performance tests for 10 million calculations using all four methods for finding a distance in a UTM projection. Below is the processing time for each method:

Euclidean 4.38 s
Average / Linear 4.49 s
Simpsons 6.94 s
Integral 9.26 s

Bottom line of this: If you want performance and distance is small use the averaging method, if you want accuracy, go with the Integral method, or if you want both, use the Simpsons approach.

Distance calculation - example

start point : { E=400,000 ; N=6,100,000 }

end point : { E=600,000 ; N=6,250,000 }

  • Distance using Simpsons rule: 249,942.5569m (error: 0.0003m)
  • Distance with approximate scale reduction: 249936.0046m (error: 6.54m)
  • Distance without scale reduction: 250,000m (error=57.46m)
Where is scale reduction=1?

Solving SR(x) =1 gives us the distance from the meridian where the scale reduction is exactly 1:

Area calculation To calculate the area of a polygon, apply the scale reduction to the line segments and then calculate the area as you normally would.Creating perfect circles

To create a circle with a fixed radius taking scale reduction into account will give you an irregular shape in the projection. You can use the approximate method (using center as scale factor) which will give you an ellipse since the scale factor for x is constant in this case. The formula for approximate circle becomes r2=(SR*x)2+y2 where SR is the scale reduction at the center of the circle.

For "accurate" circles SR becomes dependent on X. To find this solution you need to isolate e1 from the above equations to find the scale corrected distance. My math is too rusty to figure that one out yet, but if you know, feel free to post it in the comments.

You can download a C# class that performs all four types of distance calculations here.