verbose-mode.de

Wir kennen das alle: Auf eine statische Karte, die aus OSM Tiles besteht sollen Punkte gezeichnet werden. Die Punkte liegen aber nur in Lat/Lon bezogen auf WGS84 vor. Passiert mir 3 mal täglich. Also, was tun?
Immer im Hinterkopf behalten das Latitude und Longitude Winkel auf einer Kugel beschreiben und keine kartesischen Koordinaten sind. Um Lat/Lon auf OSM-Karten (und Google Maps, Bing Maps, …) zu zeichnen, müssen sie daher erst in die Ebene projiziert werden; die meisten Kartendienst nutzen den Spherical Mercator.

Latitude / Longitude to Spherical Mercator

double* latLonToMercator( double lat, double lon ) const
{
	
	double* returnValue = new double[2];
	double wgs84Radius = 6378137; //Equatorial Radius, WGS84

	double shift  = M_PI * wgs84Radius;
	double x      = lon * shift / 180;
	double y      = log(tan((90 + lat) * M_PI / 360)) / (M_PI / 180);
	y = y * shift / 180;

	returnValue[0] = x;
	returnValue[1] = y;
	return returnValue;
}

Mit den so errechneten werden kann man endlich den Punkt auf eine OSM-Karte zeichnen – allerdings nicht auf einen statischen Kartenausschnitt. Eine recht simple Lösung ist die ebene Ähnlichkeitstransformation. Allerdings muss man dafür 2 identische Punkte in beiden Koordinatensystemen kennen – das funktioniert nicht, wenn die Karte aus nur einem Tile besteht.
Besteht die Karte aus mindestens 2 Tiles, kann man für die Pixelkoordinaten (0;0) und (0;256) (bzw (256;0)) die Lat/Lon Werte errechnen.

Obere linke Ecke in Lat/Lon einer Tile

// x -> lon (obere linke Ecke)
double upperLeftLon( double x, int zoom ) const
{
	return (x / pow(2,zoom) * 360) - 180;
}

// y -> lat (obere linke Ecke)
double upperLeftLat( double y, int zoom ) const
{
	return  atan(sinh(M_PI * ( 1-( 2*y / pow(2,zoom)  ))  )    )   * (180/M_PI);
}

Mit diesen 4 Startpunkten können dann alle weiteren Punkte in Pixelkoordinaten umgerechnet werden.

ebene Ähnlichkeitstransformation

double* ViewTrans::transformation( double urspX1, double urspY1, double urspX2, double urspY2, double zielX1, double zielY1, double zielX2, double zielY2, double transX, double transY ) const
{
/* 'urspX/Y' sind die Koordinaten im Ursprungssystem (also hier Lat/Lon); zielX/Y die dementsprechenden Punkte im Zielsystem (also Bildkoordinaten).
    transX / transY sind die Koordinaten des Punktes, welcher transformiert werden soll
*/
	double deltaUrspX = urspX2 - urspX1;
	double deltaUrspY = urspY2 - urspY1;
	double deltaZielX = zielX2 - zielX1;
	double deltaZielY = zielY2 - zielY1;
		double* returnValue = new double[2];
	double nenner = pow(deltaUrspX, 2) + pow(deltaUrspY,2);

	double a = (deltaUrspX*deltaZielX + deltaUrspY*deltaZielY) / nenner;
	double o = (deltaUrspX*deltaZielY + deltaUrspY+deltaZielX) / nenner;

	double massstab = sqrt( pow(a,2) + pow(o,2));

	double phi = acos(a/massstab); // drehwinkel

	double yNull = zielY1 - o*urspX1 - a * urspY1;
	double xNull =zielX1 - a*urspX1 + o* urspY1;
	


	double yNeu = floor(abs(yNull + o* transX + a*transY));
	double xNeu = floor(abs(xNull + a* transX - o*transY));
	returnValue[0] = xNeu;
	returnValue[1] = yNeu;

	return returnValue;
}

Endlich eine Lösung für dieses Alltagsproblem! (Okay, der C++ Code ist nicht wirklich schön …)

Wir kennen das alle: Auf eine statische Karte, die aus OSM Tiles besteht sollen Punkte gezeichnet werden. Die Punkte liegen aber nur in Lat/Lon bezogen auf WGS84 vor. Passiert mir...