Dans l'article, je vais vous parler des résultats de ma "recherche", composer un algorithme pour convertir une carte raster arbitraire en tuiles compréhensibles pour les applications et, en cours de route, me familiariser avec des concepts tels qu'un ellipsoïde, une donnée, un système de coordonnées, une projection.
Notre Terre n'a pas la forme d'une boule, ni même la forme d'un ellipsoïde, cette figure complexe s'appelle un géoïde. Le fait est que les masses à l'intérieur de la Terre ne sont pas uniformément réparties, donc à certains endroits, la Terre est légèrement concave, dans d'autres, elle est légèrement convexe. Si nous prenons le territoire d'un pays ou d'un continent séparé, sa surface avec une précision suffisante peut être alignée avec la surface d'un ellipsoïde, si le centre de cet ellipsoïde est légèrement décalé le long de trois axes de coordonnées par rapport au centre de masse de la Terre. Un tel ellipsoïde est appelé un ellipsoïde de référence, il convient pour décrire uniquement la zone locale de la Terre pour laquelle il a été créé. À de grandes distances de ce site, il peut avoir un très grand écart avec la surface de la Terre. Un ellipsoïde dont le centre coïncide avec le centre de masse de la Terre est appelé ellipsoïde terrestre commun. Clair,que l'ellipsoïde de référence décrit mieux sa zone locale de la Terre que le terrestre général, mais le terrestre général convient à toute la surface de la Terre.
Pour décrire l'ellipsoïde, seules deux valeurs indépendantes suffisent: le rayon équatorial (généralement désigné par a) et le rayon polaire (b), mais au lieu de la deuxième valeur indépendante, la contraction polaire f = (ab) / a est généralement utilisée. C'est la première chose dont nous avons besoin dans notre algorithme en tant qu'objet - un ellipsoïde. Pour différentes parties de la Terre à différentes années, différents chercheurs ont calculé de nombreux ellipsoïdes de référence, les informations les concernant sont données sous forme de données: a (en mètres) et 1 / f (sans dimension). Curieusement, pour l'ellipsoïde terrestre commun, il existe également de nombreuses variantes différentes (différentes a, f), mais la différence n'est pas très forte, elle est principalement due à la différence dans les méthodes de détermination de a et f.
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,
};
L'exemple montre deux ellipsoïdes: le WGS84 terrestre commun, utilisé en navigation par satellite, et l'ellipsoïde de référence Krasovsky, applicable au territoire de l'Europe et de l'Asie.
Prenons un autre point intéressant, le fait est que la forme de la Terre est lente, mais changeante, de sorte que l'ellipsoïde, qui décrit aujourd'hui bien la surface, dans cent ans peut être loin de la réalité. Cela n'a pas grand-chose à voir avec l'ellipsoïde terrestre commun, car écarts dans sa même erreur, mais pertinents pour l'ellipsoïde de référence. Nous arrivons ici à un autre concept - une donnée. Le datum est un ensemble de paramètres d'un ellipsoïde (a, f), ses déplacements à l'intérieur de la Terre (pour un ellipsoïde de référence), fixés à un certain moment dans le temps. Plus précisément, la donnée peut ne pas nécessairement décrire un ellipsoïde, il peut s'agir de paramètres d'une figure plus complexe, par exemple une quasi-géoïde.
Sûrement la question s'est déjà posée: comment passer d'un ellipsoïde ou d'une donnée à une autre? Pour ce faire, chaque ellipsoïde doit avoir un système de coordonnées géodésiques: latitude et longitude (phi, lambda), la transition est réalisée en traduisant les coordonnées d'un système de coordonnées à un autre.
Il existe différentes formules pour transformer les coordonnées. Vous pouvez d'abord traduire les coordonnées géodésiques d'un système de coordonnées en coordonnées tridimensionnelles X, Y, Z, effectuer des décalages et des rotations avec elles, puis convertir les coordonnées tridimensionnelles résultantes en coordonnées géodésiques dans un autre système de coordonnées. Vous pouvez le faire directement. Parce que Puisque toutes les formules sont des séries convergentes infinies, elles sont généralement limitées à quelques membres de la série pour obtenir la précision requise. A titre d'exemple, nous utiliserons les transformées de Helmert, ce sont des transformations utilisant une transition vers des coordonnées tridimensionnelles, elles sont constituées des trois étapes décrites ci-dessus. Pour les transformations, en plus de deux ellipsoïdes, 7 paramètres supplémentaires sont nécessaires: trois décalages le long de trois axes, trois angles de rotation et un facteur d'échelle. Comme vous pouvez le deviner, tous les paramètres peuvent être extraits des datums.Mais dans l'algorithme, nous n'utiliserons pas un tel objet comme une donnée, mais à la place nous allons introduire un objet de transition d'un système de coordonnées à un autre, qui contiendra: des références à deux ellipsoïdes et 7 paramètres de transformation. Ce sera le deuxième objet de notre algorithme.
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},
};
L'exemple montre les paramètres de conversion du système de coordonnées Pulkovo 1942 en système de coordonnées WGS84. Les paramètres de transformation eux-mêmes sont un sujet distinct, pour certains systèmes de coordonnées, ils sont ouverts, pour d'autres, ils sont sélectionnés de manière empirique, leurs valeurs peuvent donc légèrement différer selon les sources.
En plus des paramètres, une fonction de transformation est également nécessaire, elle peut être directe et pour une transformation dans la direction opposée, nous n'avons besoin que d'une transformation dans la direction opposée. Je vais sauter des tonnes de matan et donner ma fonction optimisée.
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;
}
D'où vient tout cela? Dans un langage plus compréhensible, des formules peuvent être trouvées dans le projet proj4, mais depuis J'avais besoin d'optimisation pour la vitesse d'exécution, puis tous les calculs du sinus de la somme des angles ont été transformés par les formules, les exponentiations ont été optimisées par des placages entre parenthèses et les constantes ont été calculées séparément.
Maintenant, pour nous rapprocher de la tâche initiale de création de tuiles, nous devons envisager un système de coordonnées appelé WebMercator. Ce système de coordonnées est utilisé dans l'application OsmAnd et sur le Web, par exemple dans Google maps et dans OpenStreetMap. WebMercator est une projection Mercator construite sur une sphère. Les coordonnées de cette projection sont les coordonnées du pixel X, Y et dépendent de l'échelle Z, pour une échelle zéro, toute la surface de la Terre (jusqu'à environ 85 degrés de latitude) est placée sur une tuile de 256x256 pixels, les coordonnées X, Y changent de 0 à 255, en commençant par la gauche coin supérieur, pour l'échelle 1 - déjà 4 carreaux, X, Y - de 0 à 511 et ainsi de suite.
Les fonctions suivantes sont utilisées pour convertir de WebMercator en WGS84:
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);
}
Et à la fin de la première partie de l'article, nous pouvons déjà esquisser le début de notre algorithme pour créer une tuile. Chaque carreau de 256x256 pixels est adressé par trois valeurs: x, y, z, la relation avec les coordonnées X, Y, Z est très simple: x = (int) (X / 256); y = (entier) (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 */
}
}
}
Les coordonnées dans SK42 sont déjà des coordonnées converties dans le système de coordonnées de la carte, il reste maintenant à trouver un pixel sur la carte par ces coordonnées et à entrer sa couleur dans un pixel de tuile aux coordonnées j, i. Ce sera le prochain article, dans lequel nous parlerons des projections géodésiques et des transformations affines.