En la parte anterior nos detuvimos en el hecho de que es necesario obtener el color de un píxel del mapa ráster original por coordenadas geodésicas (latitud,
longitud), ya traducidas al SC de este mapa.
Las coordenadas geodésicas se establecen en la superficie del elipsoide, y las coordenadas de píxeles están en el plano, es decir, necesita una forma de pasar de una superficie convexa a una plana. El problema es que una superficie convexa (imagine que se le aplica algún tipo de dibujo o una cuadrícula de coordenadas) no se puede transferir a una superficie plana sin ninguna distorsión. Puede estar distorsionado: forma (ángulos), área, dimensiones lineales. Por supuesto, existe la posibilidad, y no la única, de pasar a una superficie plana sin distorsionar solo una cosa, pero el resto inevitablemente se distorsionará.
Obviamente, en escalas menores (todo el planeta, continente), las distorsiones son más pronunciadas que en las mayores (región, ciudad), y en las mayores (plano de un área pequeña), no estamos hablando de ellas en absoluto, porque la superficie del elipsoide en tales escalas ya no se puede distinguir del plano.
Aquí llegamos al concepto de proyección cartográfica. No daré ejemplos con imágenes de proyectar una esfera sobre un cilindro o un cono con un desarrollo posterior, para nosotros ahora es importante generalizar.
Una proyección de mapa es una forma matemáticamente definida de mapear la superficie de un elipsoide en un plano. En pocas palabras, estas son fórmulas matemáticas para convertir coordenadas geodésicas (latitud, longitud) en planas cartesianas, justo lo que necesitamos.
Se ha inventado una gran variedad de proyecciones cartográficas, todas encuentran aplicación para sus propios fines: de igual tamaño (en el que se conserva el área) para mapas políticos, mapas de suelo, conforme (en el que se conserva la forma) - para navegación, equidistante (manteniendo la escala en la dirección elegida) - para computadora procesamiento y almacenamiento de matrices de geodatos. También hay proyecciones que combinan los rasgos antes mencionados en determinadas proporciones, hay proyecciones completamente esotéricas. Se pueden encontrar ejemplos de proyecciones en Wikipedia en la página Lista de proyecciones de mapas.
Para cada proyección, se derivan fórmulas exactas o en forma de suma de series convergentes infinitas, y en este último caso puede haber incluso varias opciones. Las fórmulas de proyección convierten la latitud y la longitud a coordenadas cartesianas, normalmente el metro se utiliza como unidad de medida en dichas coordenadas. Esta cuadrícula rectangular de 1 metro a veces se puede trazar en un mapa (en forma de cuadrícula de kilómetros), pero en la mayoría de los casos no se traza. Pero ahora sabemos que todavía existe de forma implícita. La escala del mapa, que se indica en el mapa, se calcula simplemente en relación con el tamaño de esta cuadrícula. Debe entenderse claramente que 1 metro en dicha cuadrícula de coordenadas corresponde a 1 metro en el suelo, no en todo el mapa, sino solo en algunos puntos, a lo largo de una determinada línea o a lo largo de las líneas en una determinada dirección,en otros puntos o direcciones aparecen distorsiones y 1 metro en el suelo no corresponde a 1 metro de la cuadrícula de coordenadas.
Una pequeña digresión. La función de la primera parte del artículo WGS84_XYZ es precisamente la transformación de coordenadas del WGS84 SC en coordenadas rectangulares, pero expresadas no en metros, sino en píxeles. Como puede ver, la fórmula es extremadamente simple, si la proyección de Mercator no se usara en una esfera, sino en un elipsoide, entonces la fórmula sería más complicada. Por eso, para facilitar la vida de los navegadores, se optó por una esfera, posteriormente se ha arraigado la proyección WebMercator con su esfera, por lo que muchas veces es criticada.
Para mis necesidades, utilicé un documento llamado "Proyecciones cartográficas utilizadas por el Servicio Geológico de los Estados Unidos" en formato pdf, que se puede encontrar en Internet. El documento proporciona instrucciones detalladas para cada proyección para facilitar la escritura de una función de transformación en un lenguaje de programación.
Sigamos escribiendo nuestro algoritmo. Implementemos una de las proyecciones populares llamada Transverse Mercator y una de sus variantes llamada proyección Gauss-Kruger.
struct TransverseMercatorParam {
struct Ellipsoid *ep;
double k; /* */
double lat0; /* ( ) */
double lon0; /* ( ) */
double n0; /* */
double e0; /* */
double zw; /* ( ) */
double zs; /* */
//
double e2__a2k2, ie2__a2k2, m0, mp, imp;
double f0, f2, f4, f6;
double m1, m2, m4, m6;
double q, q1, q2, q3, q4, q6, q7, q8;
double q11, q12, q13, q14, q15, q16, q17, q18;
// - 2
double apk, n, b, c, d;
double b1, b2, b3, b4;
};
struct TransverseMercatorParam TM_GaussKruger = {
.ep = &Ellipsoid_Krasovsky,
.zw = rad(6,0,0),
.lon0 = -rad(3,0,0),
.e0 = 5e5,
.zs = 1e6,
};
Una característica de la proyección transversal de Mercator es que es conforme, es decir, se conservan la forma de los objetos en el mapa y los ángulos, así como el hecho de que la escala se conserva a lo largo de un meridiano central. La proyección es adecuada para todo el globo, pero las distorsiones aumentan con la distancia del meridiano central, por lo que toda la superficie de la tierra se corta en franjas estrechas a lo largo de los meridianos, que se denominan zonas, para cada una de las cuales se utiliza su propio meridiano central. Ejemplos de la implementación de tales proyecciones son la proyección de Gauss-Kruger y UTM, en las que también se define el método de distribución de coordenadas sobre zonas, es decir. definido y el mismo nombre SC.
Y, de hecho, el código de las funciones de inicialización y transformación de coordenadas. En la función de inicialización, las constantes se calculan una sola vez, que serán reutilizadas por la función de conversión.
void setupTransverseMercator(struct TransverseMercatorParam *pp) {
double sin_lat, cos_lat, cos2_lat;
double q, n, rk, ak;
if (!pp->k)
pp->k = 1.0;
sin_lat = sin(pp->lat0);
cos_lat = cos(pp->lat0);
cos2_lat = cos_lat * cos_lat;
q = pp->ep->e2 / (1 - pp->ep->e2);
// n = (a-b)/(a+b)
n = (pp->ep->a - pp->ep->b) / (pp->ep->a + pp->ep->b);
rk = (pp->ep->a + pp->ep->b) * pp->k / 2;
ak = pp->ep->a * pp->k;
pp->e2__a2k2 = pp->ep->e2 / (ak * ak);
pp->ie2__a2k2 = (1 - pp->ep->e2) / (ak * ak);
pp->f6 = 1097.0/4 * n*n*n*n;
pp->f4 = (151.0/3 - 3291.0/8 * n) * n*n*n;
pp->f2 = (21.0/2 + (-151.0/3 + 5045.0/32 * n) * n) * n*n;
pp->f0 = (3.0 + (-21.0/4 + (31.0/4 - 657.0/64 * n) * n) * n) * n;
pp->m6 = rk * 315.0/4 * n*n*n*n;
pp->m4 = rk * (-70.0/3 - 945.0/8 * n) * n*n*n;
pp->m2 = rk * (15.0/2 + (70.0/3 + 1515.0/32 * n) * n) * n*n;
pp->m1 = rk * (-3.0 + (-15.0/4 + (-4.0 - 255.0/64 * n) * n) * n) * n;
// polar distance
pp->mp = rk * (1.0 + (1.0/4 + 1.0/64 * n*n) * n*n);
pp->imp = 1 / pp->mp;
pp->m0 = pp->n0 - pp->mp * pp->lat0 - sin_lat * cos_lat *
(pp->m1 + (pp->m2 + (pp->m4 + pp->m6 * cos2_lat) * cos2_lat) * cos2_lat);
pp->q = q;
pp->q1 = 1.0/6 * q*q;
pp->q2 = 3.0/8 * q;
pp->q3 = 5.0/6 * q;
pp->q4 = 1.0/6 - 11.0/24 * q;
pp->q6 = 1.0/6 * q;
pp->q7 = 3.0/5 * q;
pp->q8 = 1.0/5 - 29.0/60 * q;
pp->q11 = - 5.0/12 * q;
pp->q12 = -5.0/24 + 3.0/8 * q;
pp->q13 = - 1.0/240 * q*q;
pp->q14 = 149.0/360 * q;
pp->q15 = 61.0/720 - 63.0/180 * q;
pp->q16 = - 1.0/40 * q*q;
pp->q17 = - 1.0/60 * q;
pp->q18 = 1.0/24 + 1.0/15 * q;
// - 2
double e2 = pp->ep->e2;
pp->apk = ak * (1 + n*n / 4 + n*n*n*n / 64) / (1 + n);
pp->n = n;
pp->b = (5 - e2) * e2 * e2 / 6;
pp->c = (104 - 45 * e2) * e2 * e2 * e2 / 120;
pp->d = 1237.0/1260 * e2 * e2 * e2 * e2;
pp->b1 = (1.0/2 + (-2.0/3 + (5.0/16 + 41.0/180 * n) * n) * n) * n;
pp->b2 = (13.0/48 + (-3.0/5 + 557.0/1440 * n) * n) * n*n;
pp->b3 = (61.0/240 - 103.0/140 * n) * n*n*n;
pp->b3 = 49561.0/161280 * n*n*n*n;
}
void translateTransverseMercator(struct TransverseMercatorParam *pp, int zone,
double lat, double lon, double *ep, double *np) {
double lon2, v, m;
double k4, k6, h3, h5;
double sin_lat = sin(lat);
double cos_lat = cos(lat);
double cos2_lat = cos_lat * cos_lat;
lon -= zone * pp->zw + pp->lon0;
while (unlikely(lon <= -M_PI))
lon += 2*M_PI;
lon2 = lon * lon;
//
v = 1 / sqrt(pp->e2__a2k2 * cos2_lat + pp->ie2__a2k2);
m = ((pp->m6 * cos2_lat + pp->m4) * cos2_lat + pp->m2) * cos2_lat + pp->m1;
k4 = ((pp->q1 * cos2_lat + pp->q2) * cos2_lat + 1.0/4 ) * cos2_lat - 1.0/24;
k6 = ((pp->q3 * cos2_lat + pp->q4) * cos2_lat - 1.0/12) * cos2_lat + 1.0/720;
h3 = (( pp->q6) * cos2_lat + 1.0/3 ) * cos2_lat - 1.0/6;
h5 = ((pp->q7 * cos2_lat + pp->q8) * cos2_lat - 1.0/6 ) * cos2_lat + 1.0/120;
// ( )
*np = pp->m0 + pp->mp * lat
+ (m + v * ((k6 * lon2 + k4) * lon2 + 0.5) * lon2) * cos_lat * sin_lat;
*ep = pp->e0 + pp->zs * zone
+ ( v * ((h5 * lon2 + h3) * lon2 + 1.0) * lon ) * cos_lat;
}
En la salida de la función de transformación, tendremos coordenadas: el desplazamiento este y norte (e, n) son coordenadas cartesianas rectangulares en metros.
Ya estamos muy cerca de completar nuestro algoritmo. Solo tenemos que encontrar las coordenadas del píxel (x, y) en el archivo de mapa. Porque las coordenadas de los píxeles también son cartesianas, luego podemos encontrarlas por transformación afín (e, n) a (x, y). Volveremos a encontrar los parámetros de la transformación más afín un poco más adelante.
struct AffineParam {
double c00, c01, c02;
double c10, c11, c12;
};
void translateAffine(struct AffineParam *app, double e, double n,
double *xp, double *yp) {
*xp = app->c00 + app->c01 * e + app->c02 * n;
*yp = app->c10 + app->c11 * e + app->c12 * n;
}
Y finalmente, el algoritmo completo de creación de mosaicos:
void renderTile(ImagePtr tile, int z, unsigned long x, unsigned long y) {
int i, j;
double wlat, wlon;
ImagePtr srcimg;
double lat, lon;
double e, n;
double x, y;
for (i = 0; i < 256; ++i) {
for (j = 0; j < 256; ++j) {
XYZ_WGS84(x * 256 + j, y * 256 + i, z, &wlat, &wlon);
translateHelmertInv(&Helmert_SK42_WGS84, wlat, wlon, 0, &lat, &lon);
findSrcImg(&srcimg, lat, lon);
translateTransverseMercator(&TM_GaussKruger, srcimg->zone, lat, lon, &e, &n);
translateAffine(&srcimg->affine, e, n, &x, &y);
setPixelColor(tile, j, i, interpolatePixelColor(srcimg, x, y));
}
}
}
El resultado del trabajo para z = 12, y = 1377, x = 2391:
En el algoritmo, la función de encontrar la imagen original del mapa srcimg a partir de las coordenadas geodésicas lat, lon especificadas en el SC del mapa quedó sin describir. Creo que no habrá problemas con él y el número de la zona srcimg-> zone, pero nos detendremos en encontrar los parámetros de la transformación afín srcimg-> affine con más detalle.
En algún lugar de Internet, hace mucho tiempo, se encontró una función de este tipo para encontrar los parámetros de una transformación afín, la cito incluso con comentarios originales:
struct TiePoint {
struct TiePoint *next;
double lon, lat;
double e, n;
double x, y;
};
void setupAffine(struct AffineParam *app, struct TiePoint *plist) {
/*
* :
* x = c00 + c01 * e + c02 * n
* y = c10 + c11 * e + c12 * n
*/
struct TiePoint *pp; /* */
double a11, a12, a13; /* */
double a21, a22, a23; /* 3*3 */
double a31, a32, a33; /* */
double b1, b2, b3; /* */
int m; /* z: m=0 -> z=x, m=1 -> z=y */
double z; /* x, y */
double t; /* */
/* 2- 3 ,
. */
/* */
a11 = a12 = a13 = a22 = a23 = a33 = 0;
for (pp = plist; pp; pp = pp->next) {
a11 += 1;
a12 += pp->e;
a13 += pp->n;
a22 += pp->e * pp->e;
a23 += pp->e * pp->n;
a33 += pp->n * pp->n;
}
/* ( ) */
a21 = a12;
a31 = a13;
a12 /= a11;
a13 /= a11;
a22 -= a21 * a12;
a32 = a23 - a31 * a12;
a23 = a32 / a22;
a33 -= a31 * a13 + a32 * a23;
/* , X Y */
for (m = 0; m < 2; m++) { /* m=0 -> X, m=1 -> Y */
/* */
b1 = b2 = b3 = 0;
for (pp = plist; pp; pp = pp->next) {
z = !m ? pp->x : pp->y;
b1 += z;
b2 += pp->e * z;
b3 += pp->n * z;
}
/* */
b1 /= a11;
b2 = (b2 - a21 * b1) / a22;
b3 = (b3 - a31 * b1 - a32 * b2) / a33;
/* */
t = b2 - a23 * b3;
*(!m ? &app->c00 : &app->c10) = b1 - a12 * t - a13 * b3;
*(!m ? &app->c01 : &app->c11) = t;
*(!m ? &app->c02 : &app->c12) = b3;
}
}
En la entrada, se deben enviar al menos tres puntos de anclaje, en la salida obtenemos parámetros listos para usar. Los puntos de anclaje son puntos para los que se conocen tanto las coordenadas de píxeles (x, y) como las coordenadas de desplazamiento este y norte (e, n). Los puntos de intersección de la cuadrícula de kilómetros en el mapa original se pueden usar como tales puntos. ¿Qué pasa si no hay una cuadrícula de kilómetros en el mapa? Luego puede establecer los pares (x, y) y (lon, lat), ya que tales puntos toman los puntos de intersección de los paralelos y meridianos, siempre están en el mapa. Solo queda convertir (lon, lat) a (e, n), esto se hace mediante la función de transformación para la proyección, en nuestro caso es translateTransverseMercator ().
Como puede ver, los puntos de anclaje son necesarios para decirle al algoritmo qué parte de la superficie terrestre describe el archivo con la imagen del mapa. Dado que ambos sistemas de coordenadas eran cartesianos, no importa cuántos puntos de anclaje establezcamos y no importa qué tan lejos estén el uno del otro, la discrepancia en todo el plano del mapa estará dentro del error al determinar los puntos de anclaje. La mayoría de los errores se deben a que se utiliza la proyección, el dato o el elipsoide incorrectos (con parámetros no especificados exactamente), como resultado, las coordenadas (e, n) en la salida no se obtienen en el sistema de coordenadas cartesiano, sino en un sistema de coordenadas ligeramente curvado con respecto al cartesiano. Visualmente, esto se puede visualizar como una "hoja arrugada". Naturalmente, aumentar el número de puntos de anclaje no resuelve este problema. El problema se puede resolver ajustando los parámetros de proyección, datum y elipsoide,en este caso, una gran cantidad de puntos de anclaje le permitirá suavizar la “hoja” con mayor precisión y no perder las áreas no lisas.
Y al final del artículo te diré cómo usar mosaicos prefabricados en OpenLayers y en qué forma prepararlos para el programa OsmAnd.
Para OpenLayers, solo necesita colocar los mosaicos en la web y nombrarlos para que pueda resaltar (z, y, x) en el nombre del archivo o directorio, por ejemplo:
/tiles/topo/12_1377_2391.jpg
o, incluso mejor:
/ tiles / topo /12/1377/2391.jpg
Entonces se pueden usar así:
map = new OpenLayers.Map (...);
map.addLayer(new OpenLayers.Layer.XYZ("Topo Maps", "/tiles/topo/${z}_${y}_${x}.jpg", {
isBaseLayer: false, visibility: false,
}));
Para el programa OsmAnd, es fácil determinar el formato de cualquier archivo existente con un conjunto de mosaicos. Te contaré de inmediato los resultados. Los mosaicos se empaquetan en un archivo de base de datos sqlite, que se crea así:
CREATE TABLE info AS SELECT 99 AS minzoom, 0 AS maxzoom;
CREATE TABLE tiles (x int, y int, z int, s int, image blob, PRIMARY KEY (x,y,z,s));
CREATE INDEX IND on tiles (x,y,z,s);
UPDATE info SET minzoom=$minzoom, maxzoom=$maxzoom
La columna s siempre se llena con cero, por lo que, no entendí, la imagen se ingresa en la forma binaria original, el formato (jpg, png, gif) se pierde, pero OsmAnd lo determina por su contenido. Los mosaicos en diferentes formatos se pueden empaquetar en una base de datos. En lugar de $ minzoom y $ maxzoom, es necesario sustituir los límites de escala de los mosaicos ingresados en la base. El archivo de base de datos completo, por ejemplo, Topo.sqlitedb, se transfiere a un teléfono inteligente o tableta en el directorio osmand / tiles. Reinicie OsmAnd, en Menú -> "Configurar mapa" -> "Capa superior" aparecerá una nueva opción "Topo": este es un mapa con nuestros nuevos mosaicos.