<html>
<body>
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.<br><br>
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:<br><br>
Everything is using the WGS 84 datum.<br><br>
Osmarender (I surmise from osm@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.<br><br>
tiles@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.<br><br>
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@home and therefore using a flat
earth projection for each tile. If mapnik base layers are coming
from the mapnik library itself
(<a href="http://svn.berlios.de/svnroot/repos/mapnik/trunk" eudora="autourl">
http://svn.berlios.de/svnroot/repos/mapnik/trunk</a>) my guess from
source code inspection that it is also using a flat-earth
projection.<br><br>
Yahoo aerial photo imagery - no idea.<br><br>
Fuller analysis below.<br><br>
Mike<br>
Manila<br><br>
<br>
<b>GPS Devices<br><br>
</b>"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<br><br>
<b>OSM Database<br><br>
</b>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.<br><br>
<b>JOSM<br><br>
</b>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.<br><br>
EPSG:4326 is a common geographic projection that refers to WGS84
<a href="http://en.wikipedia.org/wiki//wiki/Latitude">latitude</a>
/<a href="http://en.wikipedia.org/wiki//wiki/Longitude">longitude</a>
coordinates in degrees with
<a href="http://en.wikipedia.org/wiki//wiki/Greenwich">Greenwich</a> as
the central
<a href="http://en.wikipedia.org/wiki//wiki/Meridian">meridian</a>.
Directly use latitude / longitude values as x/y.<br><br>
The key source code is here:<br><br>
<a href="http://www.eigenheimstrasse.de/svn/josm/src/org/openstreetmap/josm/data/projection/Mercator.java" eudora="autourl">
http://www.eigenheimstrasse.de/svn/josm/src/org/openstreetmap/josm/data/projection/Mercator.java<br>
</a>
<a href="http://www.eigenheimstrasse.de/svn/josm/src/org/openstreetmap/josm/data/projection/Epsg4326.java" eudora="autourl">
http://www.eigenheimstrasse.de/svn/josm/src/org/openstreetmap/josm/data/projection/Epsg4326.java<br>
<br>
</a><b>JOSM Landsat plugin<br><br>
</b>This pulls landsat data from
<a href="http://wms.jpl.nasa.gov/" eudora="autourl">
http://wms.jpl.nasa.gov</a> using a keyword srs=EPSG:4326<br><br>
<b>tiles@home<br><br>
</b>On cursory examination, the Perl source code implements exactly the
same Mercator algorythm used by JOSM.<br><br>
<b>osm@home<br><br>
</b>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:<br><br>
<font color="#008000"><!-- Derive the latitude of the middle of the
map --><br>
<x-tab> </x-tab>
<xsl:variable name='middleLatitude' select='($topRightLatitude +
$bottomLeftLatitude) div 2.0'/><br>
<x-tab> </x-tab>
<!--woohoo lets do trigonometry in xslt --><br>
<x-tab> </x-tab>
<!--convert latitude to radians --><br>
<x-tab> </x-tab>
<xsl:variable name='latr' select='$middleLatitude * 3.1415926 div
180.0' /><br>
<x-tab> </x-tab>
<!--taylor series: two terms is 1% error at lat<68 and 10% error
lat<83. we probably need polar projection by then --><br>
<x-tab> </x-tab>
<xsl:variable name='coslat' select='1 - ($latr * $latr) div 2 + ($latr
* $latr * $latr * $latr) div 24' /><br>
<x-tab> </x-tab>
<xsl:variable name='projection' select='1 div $coslat' /><br><br>
<x-tab> </x-tab>
<xsl:variable name='dataWidth'
select='(number($topRightLongitude)-number($bottomLeftLongitude))*10000*$scale'
/><br>
<x-tab> </x-tab>
<xsl:variable name='dataHeight'
select='(number($topRightLatitude)-number($bottomLeftLatitude))*10000*$scale*$projection'
/><br>
<x-tab> </x-tab>
<xsl:variable name='km' select='(0.0089928*$scale*10000*$projection)'
/><br>
<x-tab> </x-tab>
<xsl:variable name='documentWidth'><br>
<x-tab> </x-tab><x-tab>
</x-tab>
<xsl:choose><br>
<x-tab> </x-tab><x-tab>
</x-tab><x-tab>
</x-tab><xsl:when
test='$dataWidth > (number(/rules/@minimumMapWidth) *
$km)'><br>
<x-tab> </x-tab><x-tab>
</x-tab><x-tab>
</x-tab><x-tab>
</x-tab><xsl:value-of
select='$dataWidth'/><br>
<x-tab> </x-tab><x-tab>
</x-tab><x-tab>
</x-tab>
</xsl:when><br>
<x-tab> </x-tab><x-tab>
</x-tab><x-tab>
</x-tab>
<xsl:otherwise><xsl:value-of
select='number(/rules/@minimumMapWidth) *
$km'/></xsl:otherwise><br>
<x-tab> </x-tab><x-tab>
</x-tab>
</xsl:choose><br>
<x-tab> </x-tab>
</xsl:variable><br><br>
</font><b>Mapnik<br><br>
</b>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.<br><br>
<a href="http://svn.berlios.de/svnroot/repos/mapnik/trunk/src/proj_transform.cpp" eudora="autourl">
http://svn.berlios.de/svnroot/repos/mapnik/trunk/src/proj_transform.cpp<br>
</a>
<a href="http://svn.berlios.de/svnroot/repos/mapnik/trunk/src/projection.cpp" eudora="autourl">
http://svn.berlios.de/svnroot/repos/mapnik/trunk/src/projection.cpp<br>
</a>
<a href="http://svn.berlios.de/svnroot/repos/mapnik/trunk/src/scale_denominator.cpp" eudora="autourl">
http://svn.berlios.de/svnroot/repos/mapnik/trunk/src/scale_denominator.cpp</a>
:<br><br>
static const double meters_per_degree = 6378137 * 2 * pi/ 360;<br>
...<br>
double denom = map.scale() / 0.00028;<br><br>
<b>Google Earth<br><br>
</b>
<a href="http://cfis.savagexi.com/articles/2006/05/03/google-maps-deconstructed" eudora="autourl">
http://cfis.savagexi.com/articles/2006/05/03/google-maps-deconstructed<br>
<br>
</a>Uses the "Simple Mercator" projection, i.e. assumes the
earth is round.<br><br>
<b>Yahoo aerial photo imagery<br><br>
</b>?? </body>
</html>