[OSM-talk] OSM et al map projections

[Mike Collinson]

Can anyone reviewing the following tell me whether it has any glaring 
errors and can fill in any of the missing data?  It may also be 
useful for the current discussion on projections.

I want to be able to overlay like with like in my own OSM-related 
data gathering and display software and have had a fun weekend 
working out what datums and map projections are being use where; my 
conclusions so far:

Everything is using the WGS 84 datum.

Osmarender (I surmise from osm at home) assumes the earth is flat, it 
maps lat/lon directly to pixels using a scaling factor from the 
middle of the image to get the x and y scales roughly equal.  JOSM in 
its default(?) configuration uses EPSG:4326 projection, which seems 
to be a formal way of saying the earth is flat.  The JOSM landsat 
plugin also uses EPSG:4326.

tiles at home, JOSM's "Mercator" projection and Google Earth uses what 
I've seen termed as "Simple Mercator".  This assumes the earth is 
perfectly round when projecting onto a flat surface, i.e. it does not 
take into account the WGS 84 spheroid.

Slippy map:  I've no real idea what the difference between the mapnik 
db, mapnik, mapnik WMS-C and osmarender base layers are.  I'd assume 
Osmarender is coming from osm at home and therefore using a flat earth 
projection for each tile.  If mapnik base layers are coming from the 
mapnik library itself 
(http://svn.berlios.de/svnroot/repos/mapnik/trunk) my guess from 
source code inspection that it is also using a flat-earth projection.

Yahoo aerial photo imagery - no idea.

Fuller analysis below.

Mike
Manila


GPS Devices

"The only coordinate system "natural" to GPS is the cartesian 
geocentric (X,Y,Z), and it is used internally by ALL GPS receivers. 
For example it is easy to convince a sirf-based device to provide 
this information as output. (X,Y,Z) is usually converted by the NMEA 
machinery into (long,lat,h;WGS84) and the OSM database stores only 
(long,lat) by dropping 'h' and implicitly assuming WGS84 datum. If 
you are taking the (long,lat) pair with another datum, the error can 
be hundreds of meters." Johnny Doe

OSM Database

The OSM database simple stores lat/lons so does not formally define 
any datum. However, since it is based primarilly on GPS data, it is 
implicitly WGS 84 datum-based.

JOSM

The JOSM Editor grabs WGS 84 datum-based OSM nodes and/or or your WGS 
84 datum-based GPS track data presents them to you in one of two 
projections onto the flat screen of your display: EPSG:4326 and Mercator.

EPSG:4326 is a common geographic projection that refers to WGS84 
<http://en.wikipedia.org/wiki//wiki/Latitude>latitude/<http://en.wikipedia.org/wiki//wiki/Longitude>longitude 
coordinates in degrees with 
<http://en.wikipedia.org/wiki//wiki/Greenwich>Greenwich as the 
central 
<http://en.wikipedia.org/wiki//wiki/Meridian>meridian.  Directly use 
latitude / longitude values as x/y.

The key source code is here:

http://www.eigenheimstrasse.de/svn/josm/src/org/openstreetmap/josm/data/projection/Mercator.java
http://www.eigenheimstrasse.de/svn/josm/src/org/openstreetmap/josm/data/projection/Epsg4326.java

JOSM Landsat plugin

This pulls landsat data from http://wms.jpl.nasa.gov using a keyword 
srs=EPSG:4326

tiles at home

On cursory examination, the Perl source code implements exactly the 
same Mercator algorythm used by JOSM.

osm at home

Direct lat/lon to pixel mapping with scaling to get the x and y 
scales to match. The key source code appears to be this:

<!-- Derive the latitude of the middle of the map -->
         <xsl:variable name='middleLatitude' 
select='($topRightLatitude + $bottomLeftLatitude) div 2.0'/>
         <!--woohoo lets do trigonometry in xslt -->
         <!--convert latitude to radians -->
         <xsl:variable name='latr' select='$middleLatitude * 
3.1415926 div 180.0' />
         <!--taylor series: two terms is 1% error at lat<68 and 10% 
error lat<83. we probably need polar projection by then -->
         <xsl:variable name='coslat' select='1 - ($latr * $latr) div 
2 + ($latr * $latr * $latr * $latr) div 24' />
         <xsl:variable name='projection' select='1 div $coslat' />

         <xsl:variable name='dataWidth' 
select='(number($topRightLongitude)-number($bottomLeftLongitude))*10000*$scale' 
/>
         <xsl:variable name='dataHeight' 
select='(number($topRightLatitude)-number($bottomLeftLatitude))*10000*$scale*$projection' 
/>
         <xsl:variable name='km' 
select='(0.0089928*$scale*10000*$projection)' />
         <xsl:variable name='documentWidth'>
                 <xsl:choose>
                         <xsl:when test='$dataWidth > 
(number(/rules/@minimumMapWidth) * $km)'>
                                 <xsl:value-of select='$dataWidth'/>
                         </xsl:when>
                         <xsl:otherwise><xsl:value-of 
select='number(/rules/@minimumMapWidth) * $km'/></xsl:otherwise>
                 </xsl:choose>
         </xsl:variable>

Mapnik

I couldn't track down what is happening, the code is not exactly 
written for linear traceibility.  The key files seem to be the 
following, if so I'd stab a guess that it is using an "earth is flat" 
direct-from-lat/lon projection but using the major axis of the WGS 84 
spheroid for degrees to metres conversion.

http://svn.berlios.de/svnroot/repos/mapnik/trunk/src/proj_transform.cpp
http://svn.berlios.de/svnroot/repos/mapnik/trunk/src/projection.cpp
http://svn.berlios.de/svnroot/repos/mapnik/trunk/src/scale_denominator.cpp :

static const double meters_per_degree = 6378137 * 2 * pi/ 360;
...
double denom = map.scale() / 0.00028;

Google Earth

http://cfis.savagexi.com/articles/2006/05/03/google-maps-deconstructed

Uses the "Simple Mercator" projection, i.e. assumes the earth is round.

Yahoo aerial photo imagery

??  
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.openstreetmap.org/pipermail/talk/attachments/20070121/fab5c3bd/attachment.html>