En el artículo, le contaré los resultados de mi "investigación", compondré un algoritmo para convertir un mapa ráster arbitrario en mosaicos que sean comprensibles para las aplicaciones y, en el camino, me familiarizaré con conceptos tales como elipsoide, datum, sistema de coordenadas, proyección.
Nuestra Tierra no tiene la forma de una bola, y ni siquiera la forma de un elipsoide, esta figura compleja se llama geoide. El hecho es que las masas dentro de la Tierra no están distribuidas de manera uniforme, por lo que en algunos lugares la Tierra es ligeramente cóncava, en otros es ligeramente convexa. Si tomamos el territorio de un país o continente separado, entonces su superficie con suficiente precisión puede alinearse con la superficie de algún elipsoide, si el centro de este elipsoide se desplaza ligeramente a lo largo de tres ejes de coordenadas con respecto al centro de masa de la Tierra. Tal elipsoide se llama elipsoide de referencia, es adecuado para describir solo el área local de la Tierra para la que fue creado. A grandes distancias de este sitio, puede haber una gran discrepancia con la superficie de la Tierra. Un elipsoide cuyo centro coincide con el centro de masa de la Tierra se llama elipsoide terrestre común. Claro,que el elipsoide de referencia describe mejor su área local de la Tierra que el terrestre general, pero el terrestre general es adecuado para toda la superficie de la Tierra.
Para describir el elipsoide, solo dos valores independientes son suficientes: el radio ecuatorial (generalmente denotado por a) y el radio polar (b), pero en lugar del segundo valor independiente, generalmente se usa la contracción polar f = (ab) / a. Esto es lo primero que necesitamos en nuestro algoritmo como objeto: un elipsoide. Para diferentes partes de la Tierra en diferentes años, diferentes investigadores han calculado muchos elipsoides de referencia, la información sobre ellos se da en forma de datos: a (en metros) y 1 / f (adimensional). Curiosamente, para el elipsoide terrestre común, también hay muchas variantes diferentes (diferentes a, f), pero la diferencia no es muy fuerte, se debe principalmente a la diferencia en los métodos para determinar a y 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,
};
El ejemplo muestra dos elipsoides: el común terrestre WGS84, utilizado en la navegación por satélite, y el elipsoide de referencia de Krasovsky, aplicable al territorio de Europa y Asia.
Considere otro punto interesante, el caso es que la forma de la Tierra es lenta, pero cambiante, por lo que el elipsoide, que hoy describe bien la superficie, en cien años puede estar lejos de la realidad. Esto tiene poco que ver con el elipsoide terrestre común, ya que desviaciones dentro de su mismo error, pero relevantes para el elipsoide de referencia. Aquí llegamos a otro concepto: un dato. Datum es un conjunto de parámetros de un elipsoide (a, f), sus desplazamientos dentro de la Tierra (para un elipsoide de referencia), fijados en un momento determinado. Más precisamente, el datum puede no necesariamente describir un elipsoide, pueden ser parámetros de una figura más compleja, por ejemplo, un cuasigeoide.
Seguro que ya ha surgido la pregunta: ¿cómo pasar de un elipsoide o datum a otro? Para ello, cada elipsoide debe tener un sistema de coordenadas geodésicas: latitud y longitud (phi, lambda), la transición se realiza traduciendo coordenadas de un sistema de coordenadas a otro.
Existen varias fórmulas para transformar coordenadas. Primero puede traducir las coordenadas geodésicas en un sistema de coordenadas en coordenadas tridimensionales X, Y, Z, realizar cambios y giros con ellas, y luego traducir las coordenadas tridimensionales resultantes en coordenadas geodésicas en otro sistema de coordenadas. Puedes hacerlo directamente. Porque todas las fórmulas son series convergentes infinitas, por lo que normalmente se limitan a unos pocos miembros de la serie para lograr la precisión requerida. Como ejemplo, usaremos las transformadas de Helmert, estas son transformaciones que utilizan una transición a coordenadas tridimensionales, consisten en las tres etapas descritas anteriormente. Para las transformaciones, además de dos elipsoides, necesita 7 parámetros más: tres cambios a lo largo de tres ejes, tres ángulos de rotación y un factor de escala. Como puede adivinar, todos los parámetros se pueden extraer de los datums.Pero en el algoritmo no usaremos tal objeto como datum, sino que introduciremos un objeto de transición de un sistema de coordenadas a otro, que contendrá: referencias a dos elipsoides y 7 parámetros de transformación. Este será el segundo objeto de nuestro algoritmo.
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},
};
El ejemplo muestra los parámetros para convertir del sistema de coordenadas Pulkovo 1942 al sistema de coordenadas WGS84. Los parámetros de transformación en sí mismos son un tema aparte, para algunos sistemas de coordenadas están abiertos, para otros se seleccionan empíricamente, por lo tanto, sus valores pueden diferir ligeramente en diferentes fuentes.
Además de los parámetros, también se necesita una función de transformación, puede ser directa y para la transformación en la dirección opuesta, solo necesitamos una transformación en la dirección opuesta. Me saltaré toneladas de matan y daré mi función optimizada.
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;
}
¿De dónde viene todo esto? En un lenguaje más comprensible, las fórmulas se pueden encontrar en el proyecto proj4, pero como Necesitaba optimización para la velocidad de ejecución, luego los cálculos del seno de la suma de los ángulos se transformaron mediante las fórmulas, las exponenciaciones se optimizaron mediante revestimientos entre paréntesis y las constantes se calcularon por separado.
Ahora, para acercarnos a completar la tarea original de crear mosaicos, debemos considerar un sistema de coordenadas llamado WebMercator. Este sistema de coordenadas se usa en la aplicación OsmAnd y en la web, por ejemplo, en mapas de Google y en OpenStreetMap. WebMercator es una proyección de Mercator construida sobre una esfera. Las coordenadas en esta proyección son las coordenadas del píxel X, Y y dependen de la escala Z, para una escala cero, toda la superficie de la tierra (hasta aproximadamente 85 grados de latitud) se coloca en un mosaico de 256x256 píxeles, las coordenadas X, Y cambian de 0 a 255, comenzando desde la izquierda. esquina superior, para escala 1 - ya 4 mosaicos, X, Y - de 0 a 511 y así sucesivamente.
Las siguientes funciones se utilizan para convertir de WebMercator a 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);
}
Y al final de la primera parte del artículo, ya podemos esbozar el comienzo de nuestro algoritmo para crear un mosaico. Cada mosaico de 256x256 píxeles se direcciona mediante tres valores: x, y, z, la relación con las coordenadas X, Y, Z es muy simple: 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 */
}
}
}
Las coordenadas en SK42 ya son coordenadas convertidas en el sistema de coordenadas del mapa, ahora queda por encontrar un píxel en el mapa por estas coordenadas e ingresar su color en un píxel de mosaico en las coordenadas j, i. Este será el próximo artículo, en el que hablaremos de proyecciones geodésicas y transformaciones afines.