In dem Artikel werde ich über die Ergebnisse meiner "Forschung" sprechen, einen Algorithmus zum Konvertieren einer beliebigen Rasterkarte in Kacheln erstellen, die für Anwendungen verständlich sind, und dabei Konzepte wie Ellipsoid, Datum, Koordinatensystem, Projektion kennenlernen.
Unsere Erde hat nicht die Form einer Kugel und nicht einmal die Form eines Ellipsoids. Diese komplexe Figur wird Geoid genannt. Tatsache ist, dass die Massen innerhalb der Erde nicht gleichmäßig verteilt sind, so dass die Erde an einigen Stellen leicht konkav ist, an anderen leicht konvex. Wenn wir das Territorium eines separaten Landes oder Kontinents einnehmen, kann seine Oberfläche mit ausreichender Genauigkeit mit der Oberfläche eines Ellipsoids ausgerichtet werden, wenn der Mittelpunkt dieses Ellipsoids entlang drei Koordinatenachsen relativ zum Massenschwerpunkt der Erde leicht verschoben ist. Ein solches Ellipsoid wird als Referenzellipsoid bezeichnet und eignet sich zur Beschreibung nur des lokalen Bereichs der Erde, für den es erstellt wurde. In großen Entfernungen von diesem Ort kann es zu einer sehr großen Diskrepanz mit der Erdoberfläche kommen. Ein Ellipsoid, dessen Mittelpunkt mit dem Massenmittelpunkt der Erde zusammenfällt, wird als gemeinsames terrestrisches Ellipsoid bezeichnet. Klar,dass das Referenzellipsoid seinen lokalen Teil der Erde besser beschreibt als das allgemeine terrestrische, aber das allgemeine terrestrische ist für die gesamte Erdoberfläche geeignet.
Zur Beschreibung des Ellipsoids reichen nur zwei unabhängige Werte aus: der äquatoriale Radius (normalerweise mit a bezeichnet) und der polare Radius (b). Anstelle des zweiten unabhängigen Werts wird jedoch üblicherweise die polare Kontraktion f = (ab) / a verwendet. Dies ist das erste, was wir in unserem Algorithmus als Objekt benötigen - ein Ellipsoid. Für verschiedene Teile der Erde in verschiedenen Jahren haben verschiedene Forscher viele Referenzellipsoide berechnet. Informationen darüber werden in Form von Daten angegeben: a (in Metern) und 1 / f (dimensionslos). Seltsamerweise gibt es für das gemeinsame terrestrische Ellipsoid auch viele verschiedene Varianten (verschiedene a, f), aber der Unterschied ist nicht sehr stark, sondern hauptsächlich auf die unterschiedlichen Methoden zur Bestimmung von a und f zurückzuführen.
struct Ellipsoid {
char *name;
double a; /* () */
double b; /* () */
double al; /* (a-b)/a */
double e2; /* (a^2-b^2)/a^2 */
};
struct Ellipsoid Ellipsoid_WGS84 = {
.name = "WGS84",
.a = 6378137.0,
.al = 1.0 / 298.257223563,
};
struct Ellipsoid Ellipsoid_Krasovsky = {
.name = "Krasovsky",
.a = 6378245.0,
.al = 1.0 / 298.3,
};
Das Beispiel zeigt zwei Ellipsoide: das gemeinsame terrestrische WGS84, das in der Satellitennavigation verwendet wird, und das Krasovsky-Referenzellipsoid, das für Europa und Asien gilt.
Betrachten Sie einen weiteren interessanten Punkt: Tatsache ist, dass die Form der Erde langsam ist, sich aber ändert, sodass das Ellipsoid, das heute die Oberfläche gut beschreibt, in hundert Jahren möglicherweise weit von der Realität entfernt ist. Dies hat seitdem wenig mit dem gemeinsamen terrestrischen Ellipsoid zu tun Abweichungen innerhalb desselben Fehlers, jedoch relevant für das Referenzellipsoid. Hier kommen wir zu einem anderen Konzept - dem Datum. Datum ist eine Reihe von Parametern des Ellipsoids (a, f), seiner Verschiebungen innerhalb der Erde (für das Referenzellipsoid), die zu einem bestimmten Zeitpunkt festgelegt sind. Genauer gesagt, das Datum muss nicht unbedingt ein Ellipsoid beschreiben, es kann Parameter einer komplexeren Figur sein, beispielsweise eines Quasigeoids.
Sicherlich hat sich bereits die Frage gestellt: Wie kann man von einem Ellipsoid oder Datum zu einem anderen wechseln? Dazu muss jedes Ellipsoid ein geodätisches Koordinatensystem haben: Breiten- und Längengrad (Phi, Lambda), der Übergang erfolgt durch Übersetzen von Koordinaten von einem Koordinatensystem in ein anderes.
Es gibt verschiedene Formeln zum Transformieren von Koordinaten. Sie können zuerst geodätische Koordinaten in einem Koordinatensystem in dreidimensionale Koordinaten X, Y, Z übersetzen, Verschiebungen und Rotationen mit ihnen durchführen und dann die resultierenden dreidimensionalen Koordinaten in geodätische Koordinaten in einem anderen Koordinatensystem konvertieren. Sie können es direkt tun. weil Da alle Formeln unendlich konvergierende Reihen sind, ist sie normalerweise auf wenige Mitglieder der Reihe beschränkt, um die erforderliche Genauigkeit zu erreichen. Als Beispiel werden wir die Helmert-Transformationen verwenden. Dies sind Transformationen, die einen Übergang zu dreidimensionalen Koordinaten verwenden. Sie bestehen aus den drei oben beschriebenen Stufen. Für Transformationen benötigen Sie zusätzlich zu zwei Ellipsoiden 7 weitere Parameter: drei Verschiebungen entlang drei Achsen, drei Drehwinkel und einen Skalierungsfaktor. Wie Sie vielleicht erraten haben, können alle Parameter aus Bezugspunkten extrahiert werden.Im Algorithmus werden wir jedoch kein solches Objekt als Bezugspunkt verwenden, sondern stattdessen ein Übergangsobjekt von einem Koordinatensystem in ein anderes einführen, das Folgendes enthält: Verweise auf zwei Ellipsoide und 7 Transformationsparameter. Dies wird das zweite Objekt unseres Algorithmus sein.
struct HelmertParam {
char *src, *dst;
struct Ellipsoid *esp;
struct Ellipsoid *edp;
struct {
double dx, dy, dz;
double wx, wy, wz;
double ms;
} p;
//
double a, da;
double e2, de2;
double de2__2, dxe2__2;
double n, n__2e2;
double wx_2e2__ro, wy_2e2__ro;
double wx_n__ro, wy_n__ro;
double wz__ro, ms_e2;
};
struct HelmertParam Helmert_SK42_WGS84 = {
.src = "SK42",
.dst = "WGS84",
.esp = &Ellipsoid_Krasovsky,
.edp = &Ellipsoid_WGS84,
// SK42->PZ90->WGS84 ( 51794-2001)
.p = {23.92, -141.27, -80.9, 0, -0.35, -0.82, -0.12e-6},
};
Das Beispiel zeigt die Parameter für die Konvertierung vom Pulkovo 1942-Koordinatensystem in das WGS84-Koordinatensystem. Die Transformationsparameter selbst sind ein separates Thema, für einige Koordinatensysteme sind sie offen, für andere werden sie empirisch ausgewählt, daher können sich ihre Werte in verschiedenen Quellen geringfügig unterscheiden.
Zusätzlich zu den Parametern wird auch eine Transformationsfunktion benötigt, die direkt sein kann und für die Transformation in die entgegengesetzte Richtung benötigen wir nur eine Transformation in die entgegengesetzte Richtung. Ich werde Tonnen von Mathe überspringen und meine optimierte Funktion geben.
void setupHelmert(struct HelmertParam *pp) {
pp->a = (pp->edp->a + pp->esp->a) / 2;
pp->da = pp->edp->a - pp->esp->a;
pp->e2 = (pp->edp->e2 + pp->esp->e2) / 2;
pp->de2 = pp->edp->e2 - pp->esp->e2;
pp->de2__2 = pp->de2 / 2;
pp->dxe2__2 = pp->de2__2 + pp->e2 * pp->da / pp->a;
pp->n = 1 - pp->e2;
pp->n__2e2 = pp->n / pp->e2 / 2;
pp->wx_2e2__ro = pp->p.wx * pp->e2 * 2 * rad(0,0,1);
pp->wy_2e2__ro = pp->p.wy * pp->e2 * 2 * rad(0,0,1);
pp->wx_n__ro = pp->p.wx * pp->n * rad(0,0,1);
pp->wy_n__ro = pp->p.wy * pp->n * rad(0,0,1);
pp->wz__ro = pp->p.wz * rad(0,0,1);
pp->ms_e2 = pp->p.ms * pp->e2;
}
void translateHelmertInv(struct HelmertParam *pp,
double lat, double lon, double h, double *latp, double *lonp) {
double sin_lat, cos_lat;
double sin_lon, cos_lon;
double q, n;
if (unlikely(!pp)) {
*latp = lat;
*lonp = lon;
return;
}
sin_lat = sin(lat);
cos_lat = cos(lat);
sin_lon = sin(lon);
cos_lon = cos(lon);
q = 1 / (1 - pp->e2 * sin_lat * sin_lat);
n = pp->a * sqrt(q);
*latp = lat
- ((n * (q * pp->de2__2 + pp->dxe2__2) * sin_lat + pp->p.dz) * cos_lat
- (pp->p.dx * cos_lon + pp->p.dy * sin_lon) * sin_lat
) / (n * q * pp->n + h)
+ (pp->wx_2e2__ro * sin_lon - pp->wy_2e2__ro * cos_lon)
* (cos_lat * cos_lat + pp->n__2e2)
+ pp->ms_e2 * sin_lat * cos_lat;
*lonp = lon
+ ((pp->p.dx * sin_lon - pp->p.dy * cos_lon) / (n + h)
- (pp->wx_n__ro * cos_lon + pp->wy_n__ro * sin_lon) * sin_lat
) / cos_lat
+ pp->wz__ro;
}
Woher kommt das alles? In einer verständlicheren Sprache finden sich Formeln im proj4-Projekt, aber seitdem Ich brauchte eine Optimierung für die Ausführungsgeschwindigkeit, dann wurden alle Berechnungen des Sinus der Summe der Winkel durch die Formeln transformiert, die Potenzierungen durch Verblendungen in Klammern optimiert und die Konstanten separat berechnet.
Um der ursprünglichen Aufgabe des Erstellens von Kacheln näher zu kommen, müssen Sie ein Koordinatensystem namens WebMercator in Betracht ziehen. Dieses Koordinatensystem wird in der OsmAnd-Anwendung und im Web verwendet, beispielsweise in Google Maps und in OpenStreetMap. WebMercator ist eine Mercator-Projektion, die auf einer Kugel basiert. Die Koordinaten in dieser Projektion sind die Koordinaten des Pixels X, Y und hängen von der Z-Skala ab. Bei einer Nullskala wird die gesamte Erdoberfläche (bis zu etwa 85 Breitengraden) auf einer 256 x 256 Pixel großen Kachel platziert. Die X, Y-Koordinaten ändern sich von links nach links von 0 auf 255 obere Ecke, für Skala 1 - bereits 4 Kacheln, X, Y - von 0 bis 511 und so weiter.
Die folgenden Funktionen werden zum Konvertieren von WebMercator nach WGS84 verwendet:
void XYZ_WGS84(unsigned x, unsigned y, int z, double *latp, double *lonp) {
double s = M_PI / ((1UL << 7) << z);
*lonp = s * x - M_PI;
*latp = asin(tanh(M_PI - s * y));
}
void WGS84_XYZ(int z, double lat, double lon, unsigned *xp, unsigned *yp) {
double s = ((1UL << 7) << z) / M_PI;
*xp = uint_round((lon + M_PI) * s);
*yp = uint_round((M_PI - atanh(sin(lat))) * s);
}
Und am Ende des ersten Teils des Artikels können wir bereits den Anfang unseres Algorithmus zum Erstellen einer Kachel skizzieren. Jede Kachel mit 256 x 256 Pixeln wird mit drei Werten adressiert: x, y, z, die Beziehung zu den Koordinaten X, Y, Z ist sehr einfach: x = (int) (X / 256); y = (int) (Y / 256); z = Z;
void renderTile(int z, unsigned long x, unsigned long y) {
int i, j;
double wlat, wlon;
double lat, lon;
for (i = 0; i < 255; ++i) {
for (j = 0; j < 255; ++j) {
XYZ_WGS84(x * 256 + j, y * 256 + i, z, &wlat, &wlon);
translateHelmertInv(&Helmert_SK42_WGS84, wlat, wlon, 0, &lat, &lon);
/* lat,lon - 42 */
}
}
}
Koordinaten in SK42 sind bereits transformierte Koordinaten in das Kartenkoordinatensystem. Jetzt muss ein Pixel auf der Karte anhand dieser Koordinaten gefunden und seine Farbe an den Koordinaten j, i in ein Kachelpixel eingegeben werden. Dies wird der nächste Artikel sein, in dem wir über geodätische Projektionen und affine Transformationen sprechen werden.