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 e₁ to e₀ (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 e₁ and e₀ 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.

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:

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.

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)

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

I have lately seen several blogposts confused about which datum and projection Microsoft’s Live Maps (Virtual Earth) and Google Maps use. As most people already know by now, they render the round earth onto a flat screen using a Mercator projection.

I think the confusion comes from that they actually use two spatial reference systems at the same time:

1. Geographic  Longitude/Latitude coordinatesystem based on the standard WGS84 datum.
2. Mercator projection using a datum based on WGS84, BUT modified to be spheric.

So when is what used?

The Javascript API’s use (1) as input when you want to add points, lines and polygons. That is, they expect you to input any geometry in geographical coordinates, and click events etc. will also return geometry in this spatial reference. This is the coordinate system most javascript developers will use. The API will automatically project it to the spheric mercator projection.

If you want to create image overlays, or roll your own tile server on top of the map, you will need to project your images into a spheric mercator projection. The JavaScript APIs are not able to do this for you.

Here’s a bit of facts about the two projections:

The valid range of (1) is: [-180,-85.05112877980659] to [180, 85.05112877980659].

The valid range of (2) is: [-20037508.3427892, -20037508.3427892] to [20037508.3427892, 20037508.3427892]

Well-known Text for (1):
GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]]

Well-known Text for (2):
PROJCS["Mercator Spheric", GEOGCS["WGS84basedSpheric_GCS", DATUM["WGS84basedSpheric_Datum", SPHEROID["WGS84based_Sphere", 6378137, 0], TOWGS84[0, 0, 0, 0, 0, 0, 0]], PRIMEM["Greenwich", 0, AUTHORITY["EPSG", "8901"]], UNIT["degree", 0.0174532925199433, AUTHORITY["EPSG", "9102"]], AXIS["E", EAST], AXIS["N", NORTH]], PROJECTION["Mercator"], PARAMETER["False_Easting", 0], PARAMETER["False_Northing", 0], PARAMETER["Central_Meridian", 0], PARAMETER["Latitude_of_origin", 0], UNIT["metre", 1, AUTHORITY["EPSG", "9001"]], AXIS["East", EAST], AXIS["North", NORTH]]

Proj.4 definition for (1):
+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs

Proj.4 definition for (2) (see here for an explanation of the weird ’nadgrids’ parameter):
+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs

So why this weird extent for the latitude? First of all, the poles in the Mercator projection extends towards infinity, so at some point they have to cut them off (and who cares about ice anyway - apart from that its melting). If you look at (2) these latitude/longitude values projected into the spheric mercator results in a perfect square, fitting perfectly with squared image tiles, that are simple to sub-divide over and over again, as you zoom in. I expect the reason for the spheric datum is for simplicity and perfomance when reprojecting points from longitude/latitude to screen coordinates. Charlie Savage also has a more mathematical approach to deriving these values.

For a quick introduction to projections, coordinate systems and datums see here.

If you want to know more about how these mapping api's work, keep an eye on Jayant's blog.

Update: We now have an offical EPSG code for the projection. See details here.

I've had a lot of requests for the .NET 2.0 based raytracer that I used in my Master Thesis to generate true orthophotos based on complex 3D surface models (not necessarily 2½D). To save me from having to dig it up all the time, I've now posted it here on my blog so you can download it from here.

Basically it will read in a text file with triangles (each triangle defined per line as comma-separated double values X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3), create a spatial index using axis-aligned boundingboxes, and perform ray/triangle intersections, returning any triangle intersecting the ray, including distance to triangle and angle of attack.

To save memory on larger DSM, all internal computation is done using floating-point insted of using double precision. This poses a problem for large values like UTM based coordinate systems, where rounding errors will occur when using floating point. The raytracer is internally applying an offset to the center of the data extents to accommodate for this, but it's worth remembering if your DSM covers a vast area. In most cases it's not a big issue though.

The raytracer scales very well (O(log n)), and can handle hundreds of thousands of faces and still raytrace in a few milliseconds (see chapter 7.3.3 of the thesis for details). Performing the required indexing first might be a little slow, so if you want to optimize this, be my guest ;-).

Some uses this engine has been involved in is generating rough worldfiles for aerial images based on absolute camera orientation and a surface model or Simulating Navstar and Galileo GPS availability in dense urban areas. The latter was actually some very cool stuff developed by some students at DTU. You could click anywhere on a map, and it would instantly simulate the visibility of Navstar and Galileo GPS satellites and give you the PDOP values throughout the day, based on a detailed 3D city model. This was used to determine where the zones of GPS based roadpricing should be (you didn't want the zone edges to be close to a place where you had poor accuracy causing people to be charged incorrectly). The raytracer was used for doing the line-of-sight analysis between a virtual car and all the satellites.

Download: Raytracer.zip (86.29 KB)

People are often mixing the above as if they were one and the same, so here’s a recap of them. One of the things you often find people saying is that “my data is in the WGS84 coordinate system”. This doesn’t really make sense, but I will get back to this later.

This is a very confusing subject, and I might have gotten a few things wrong myself, so please add a comment and I’ll update it ASAP.

Coordinate systems

A coordinate system is simply put a way of describing a spatial property relative to a center. There is more than one way of doing this:

• The Geocentric coordinate system is based on a normal (X,Y,Z) coordinate system with the origin at the center of Earth. This is the system that GPS uses internally for doing it calculations, but since this is very unpractical to work with as a human being ( due to the lack of well-known concepts of east, north, up, down) it is rarely displayed to the user but converted to another coordinate system.
  
• The Spherical or Geographic coordinate system is probably the most well-known. It is based on angles relative to a prime meridian and Equator usually as Longitude and Latitude. Heights are usually given relative to either the mean sea level or the datum (I’ll get back to the datum later).
  
• The Cartesian coordinate system is defined as a “flat” coordinate system placed on the surface of Earth.  In some projections it’s not flat in the sense that it follows the earth’s curvature in one direction and has a known scale-error in the other direction relative to the distance of the origin. The most well-known coordinate system is the Universal Transverse Mercator (UTM), but surveyors define their own little local flat coordinate systems all the time. It is very easy to work with, fairly accurate over small distances making measurements such as length, angle and area very straightforward. Cartesian coordinate systems are strongly connected to projections that I will cover later.

Sidenote: The geocentric coordinate system is strictly speaking a cartesian coordinate system too, but this is the general terms I've seen used the most when talking about world coordinate systems.

Datums and ellipsoids

Some of the common properties of the above coordinate systems are that they are all relative to the center of Earth and except the Geocentric coordinate system, uses a height system relative to the surface of the earth.

This poses two immediate problems:

• Where is the center of the earth
• What is the shape of the earth?

By now most people should know that that the earth isn’t flat (although there are still some who doubts it). If we define the surface of Earth as being at the mean sea level (often referred to as the Geoid), we don’t get a spheroid or even an ellipsoid. Because of gravitational changes often caused by large masses such as mountain ranges etc, Earth is actually very irregular with variations of +/- 100 meters. Since this is not very practical to work with as a model of earth, we usually use an ellipsoid for approximation. The ellipsoid is defined by its semi-major axis, and either the flattening of the ellipoid or the semi-minor axis.

The center and orientation of the ellipsoid is what we call the datum. So the datum defines an ellipsoid and through the use of a set of points on the ground that we relate to points on the ellipsoid, we define the center of the Earth.  This poses another problem, because continental drift moves the points used to define the points around all the time. This is why the name of a datum usually have a year in it, often referring to the position of those points January 1st of that year (although that may vary).

There are a vast amount of datums, some used for measurements all over the world, and other local datums defined so they fit very well with a local area. Some common ones are: World Geodetic Datum 1984 (WGS84), European Datum 1950 (ED50) and North American Datum 1983 (NAD83).

The most well-known is WGS84 used by the GPS systems today. It is a good approximation of the entire world and with fix-points defined almost all over the world. When it was defined they forgot to include points in Europe though, so the Europeans now have their own ETRS89, which is usually referred to as the “realization of WGS84 in Europe”. The problem here was solely because of continental drift, so they defined some points relative to WGS84 in 1989, and keeps track of the changes. In most use-cases it is of no real importance and you can use one or the other.

I mentioned earlier that people often refer to having their data in WGS84, and you see now why this doesn’t make sense. All you know from that is that the data is defined using the WGS84 datum, but you don’t know which coordinate system it uses.

Read more on Datums and Spheroids.

Projections

The earth isn’t flat, and there is no simple way of putting it down on a flat paper map (or these days onto a computer screen), so people have come up with all sorts of ingenious solutions each with their pros and cons. Some preserves area, so all objects have a relative size to each other, others preserve angles (conformal) like the Mercator projection, some try to find a good intermediate mix with only little distortion on several parameters etc. Common to them all is that they transform the world onto a flat Cartesian coordinate system, and which one to choose depends on what you are trying to show.

A common statement that I hear in GIS is the following “My map doesn’t have a projection”, but this is simply not possible (unless you have a good old rotating globe). Often people are referring to data that is in longitude/latitude and displayed on a map without having specified any projection. What happens is that the system applies the simplest projection it can: Mapping Longitude directly to X and Latitude to Y. This results in an equirectangular projection, also called the “Plate Carree” projection. It results in very heavy distortion making areas look squashed close to the poles. You can almost say that the “opposite” of the Plate Carree is the Mercator projection which stretches areas close to the poles in the opposite direction, making them look very big. Mercator is the type of projection you see used on Live maps and Google maps, but as many often mistakenly thinks, they do NOT use WGS84 for the projected map, although WGS84 is used when you directly input longitude/latitude values using their API (read more on this here).

More on projected coordinate systems

Spatial reference

The spatial reference is a combination of all the above. It defines an ellipsoid, a datum using that ellipsoid, and either a geocentric, geographic or projection coordinate system. The projection also always has a geographic coordinate system associated with it. The European Petroleum Survey Group (EPSG) has a huge set of predefined spatial references, each given a unique ID. These ID’s are used throughout the industry and you can download an Access database with all them from their website, as well as some very good documents on projection (or see the Spatial References website).

So when you hear someone saying they have their data in WGS84, you can often assume that they have longitude/latitude data in WGS84 projected using Plate Carree.  The spatial reference ID of this is EPSG:4326.

Spatial References are often defined in a Well-known format defining all these parameters. The Spatial Reference EPSG:4326 can therefore also be written as:

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

As mentioned Live/Google maps use a Mercator projection, but although their datum is based on WGS84, they use a sphere instead of an ellipsoid. This means that they use the same center and orientation as WGS84, but without applying any flattening. The spatial reference string for their projection therefore becomes:

PROJCS["Mercator Spheric", GEOGCS["WGS84based_GCS", DATUM["WGS84based_Datum", SPHEROID["WGS84based_Sphere", 6378137, 0], TOWGS84[0, 0, 0, 0, 0, 0, 0]], PRIMEM["Greenwich", 0, AUTHORITY["EPSG", "8901"]], UNIT["degree", 0.0174532925199433, AUTHORITY["EPSG", "9102"]], AXIS["E", EAST], AXIS["N", NORTH]], PROJECTION["Mercator"], PARAMETER["False_Easting", 0], PARAMETER["False_Northing", 0], PARAMETER["Central_Meridian", 0], PARAMETER["Latitude_of_origin", 0], UNIT["metre", 1, AUTHORITY["EPSG", "9001"]], AXIS["East", EAST], AXIS["North", NORTH]]

This blog is written directly by the ArcGIS Server Development Team.  We are very excited about the ArcGIS Server 9.2 release.  We have so much to share with the user community.  Topics will range from detailed programming information ("What is a page life cycle in the Web ADF?") to GIS methodology ("What is the best methodology to produce map caches?") for serving spatial data with ArcGIS Server.  One of the key parts to this blog will be the use of an example site.  We intend to post useful entries for our users which are linked directly to examples that we have built that you can access live from the Internet.

Inspired by my petty shot at making SQL Server spatial, Ricardo Stuven has been working the new CLR in SQL Server 2005 and started adding almost all basic spatial functions, and implemented some aggregates, stored procedures and table-valued functions in SQL Server. He even started an Open Source project where you can check it out (yet another open source .NET GIS project we can add to the .NET tribe... ;-))

Read his comments on his approach, including query examples and a list of features which he posted in an update on the SharpMap forum today.

...well at least according to this "white" paper from Refractions. Quote:

Open Source GIS software can be categorized into two largely independent development tribes. [...]

• The ‘C’ tribe, consisting of developers working on UMN Mapserver, GRASS, GDAL/OGR, OSSIM, Proj4, GEOS, PostGIS, QGIS, MapGuide OS and OpenEV. The ‘C’ tribe also includes users of scripting languages that bind easily to C libraries, such as Python, Perl and PHP.
• The ‘Java’ tribe, consisting of developers working on GeoTools, uDig, GeoServer, JTS, JUMP, and DeeGree.

Don't get me wrong - I think the above libraries are GREAT (I use many of them daily) but I think this is a strange way of categorizing software. First of all what about the .NET "tribe" and all the other tribes out there? I can think of several great .NET/Mono based open source GIS application and libraries (MonoGIS, Appomattox, NetTopologySuite, GeoTools.NET, SharpMap and many more). Secondly these libraries aren't even fixed to a specific language, merely a framework where you decide what .NET/Mono compatible language you decide to use when linking to these libraries. Actually you can even mix languages within the same application.

It's funny to see that people still can't see .NET as a language used for Open Source applications - probably because the big bad wolf (Microsoft) invented todays fastest growing language.

I was just reading a few articles on XAML - Microsoft new UI language for rendering vector content and user interfaces - when it hit me that this probably can get you much further in web mapping interaction than SVG can. XAML is part of the "Windows Presentation Foundation" a part of .NET Framework 3.0

My first attempt to generate XAML was done by creating a new renderer for SharpMap, and below you can see my first results in XamlPad. The data here has been created based on a SharpMap map by my quick'n'dirty XamlRenderer for SharpMap.

All it takes with the new renderer is initializing the new renderer with your map, and request a new map:

SharpMap.Xaml.Renderer XamlRenderer = new SharpMap.Xaml.Renderer(map);

string xaml = XamlRenderer.GetMap();

If you zip the generated XAML file you even get a 20% reduction over a PNG image.

Next step is to add interaction like zooming, querying etc…

I just read this days newsletter from OGC, containing some news on the release of a simple profile of the GML v3.1.1 standard. From the standard one can read the following in the introduction:

This profile defines a restricted but useful subset of XML-Schema and GML to lower the “implementation bar” of time and resources required for an organization to commit for developing software that supports GML. [...] this profile facilitates the ability to use WFS for interoperable feature data exchange with much less software development investment.

The GML standard is a big and complex standard, and I really think this is a great move by the OGC. This is a 30 page long specification (well 117 if you include the annex), and it's nothing compared to the 267 (601) pages full GML spec.

Hopefully we'll see more GML support in the future. At least this finally got me started with implementing GML/WFS support into SharpMap.

Well, while I'm at it, if you want to know more on the GML specification, I can recommend the book "Geography Markup Language, Foundation for the Geo-Web" by Ron Lake from Wiley.

Charlie Savage has a great post introducing coordinate systems in his blog. I too have had many questions on coordinate systems, projections and datums, and Charlie have made a great introduction to the world of "spatial reference systems". He promises to post some more on the topic, so keep track of his blog if you find these topics confusing (I know I did when I first was told about datums and coordinate systems).