[Talk-de] Tilekoordinaten

Dimitri Junker OSM at dimitri-junker.de
Sa Mai 8 21:09:12 UTC 2010


Hallo,

bei dem Versuch Worldfiles zu erzeugen kam die Frage nach der richtigen 
Berechnung der OSM-Tilekoordinaten auf. Lt.
<http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#tile_numbers_to_lon
..2Flat>

lat_rad = arctan(sinh(Pi * (1 - 2 * ytile / n)))
mit n = 2 ^ zoom

In taho wird fast die gleiche Formel verwendet, aber statt Pi
ProjectF(85.0511), dies ergibt: 3.1415868112787813, also fast genau Pi, aber 
eben nicht genau. Bei Zoomlevel 12 ergibt sich da z.B. ein Unterschied von 
etwa 8m für Aachen. Ich hänge mal die entsprechenden Funktionen aus taho.pl 
und der c++ Version an. Letztere habe ich aus der perl Version abgeleitet.

Wer denkt sich soetwas kompliziertes aus wenn es falsch ist? Also was ist 
richtig?

Dimitri

C++ Version aus taho.exe


//##############################################################################
//# transform the x,y tilenumbers to the S,W,N,E edge ccordinates
//##############################################################################
void CPixmap::Project(int x, int y, int zoom, CGeoRect &gr)
{
	double unit = 1. / (1<<zoom);
	double relY1 = y * unit;
	double relY2 = relY1 + unit;


	double limitY = ProjectF(85.0511);
	limitY=3.14159265;
	double rangeY = 2 * limitY;

	relY1 = limitY - rangeY * relY1;
	relY2 = limitY - rangeY * relY2;

	gr.m_n = ProjectMercToLat(relY1);
	gr.m_s = ProjectMercToLat(relY2);

	unit = 360. / (1<<zoom);
	gr.m_w = -180 + x * unit;
	gr.m_e=gr.m_w+unit;
}


//##############################################################################
//# helper sub for Project()
//##############################################################################
double CPixmap::ProjectF(double winkel)
{
  double lat = 3.14159265/180 * winkel;
  return(log(tan(lat) + 1/cos(lat)));
}
double CPixmap::ProjectMercToLat(double mercY)
{
  return( 180/3.14159265* atan(sinh(mercY)));
}





Perl Version von taho.pl:




sub Project{
  my ($X,$Y, $Zoom) = @_;

  my $Unit = 1 / (2 ** $Zoom);
  my $relY1 = $Y * $Unit;
  my $relY2 = $relY1 + $Unit;


  my $LimitY = ProjectF(85.0511);
  my $RangeY = 2 * $LimitY;

  $relY1 = $LimitY - $RangeY * $relY1;
  $relY2 = $LimitY - $RangeY * $relY2;

  my $Lat1 = ProjectMercToLat($relY1);
  my $Lat2 = ProjectMercToLat($relY2);

  $Unit = 360 / (2 ** $Zoom);
  my $Long1 = -180 + $X * $Unit;

  return(($Lat2, $Long1, $Lat1, $Long1 + $Unit)); # S,W,N,E
}

##############################################################################
# helper sub for Project()
##############################################################################
sub ProjectMercToLat($){
  my $MercY = shift();
  return( 180/pi* atan(sinh($MercY)));
}
##############################################################################
# helper sub for Project()
##############################################################################
sub ProjectF($){
  my $Lat = pi/180 * shift();
  my $Y = log(tan($Lat) + sec($Lat));
  return($Y);
}




Mehr Informationen über die Mailingliste Talk-de