SharpGIS

#GIS from a .NET developer's perspective

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.

Visual Studio 2008 JavaScript Intellisense - Take Two

...Not for namespaces either.

We all know that global variables in JavaScript API’s are bad (although they are usually OK in application-specific contexts). So therefore the solution is to create a namespace and put all your objects inside that. That way you only end up with one global, the name of the “root” namespace. You should choose a root name that probably isn’t used by anyone else, like your company name or similar.

With the Microsoft AJAX Clientside Library, Microsoft gave us registerNamespace for this thing, like for instance:

Type.registerNamespace('myCompany.Utility.Math');
myCompany.Utility.Math.Calculator = function() { }
myCompany.Utility.Math.Calculator.prototype = {
 add : function(a,b) { return a+b; }
}
myCompany.Utility.Math.Calculator.registerClass('myCompany.Utility.Math.Calculator');

Basically what Type.registerNamespace does is creating a set of objects:

myCompany = {};
myCompany.Utility = {};
myCompany.Utility.Math = {};

There's is just one problem to this. Visual Studio 2008’s intellisense is not evaluating the result of Type.registerNamespace, and thus, never realizes that those three objects exist. The result is that you don’t get any intellisense on your namespaced objects. Bummer! The workaround is either to create the namespaces yourself (like above) before calling Type.registerNamespace, or create aliases for the classes (which would re-introduce the globals we were trying to get rid of in the first place).

Of course I considered this to be a bug, and as I mentioned at the end of my previous post there was a problem with using namespaces and JavaScript Intellisense in Visual Studio 2008, that I had submitted a bug report for.

I was surprised how quickly they investigated and escalated the issue, but my biggest surprise was the result, which is pretty much this:

We are not going to support the pattern we told you to use until version 20xx.

OK, maybe a bit over-dramatized but still... It's funny because they are actually able to give me full intellisense for the namespaced objects in the MS AJAX Sys.* namespace (but they probably hardcoded that in there ;-)

With this and the fact that javascript files included by webcontrols are not intellisensed either, I’ve come to the conclusion that the hyped and long awaited JavaScript Intellisense in Orcas is pretty much useless. What do you think?