Na parte anterior, paramos no fato de que é necessário obter a cor de um pixel do mapa raster original por coordenadas geodésicas (latitude,
longitude), já traduzidas para o SC deste mapa.
As coordenadas geodésicas são definidas na superfície do elipsóide e as coordenadas dos pixels estão no plano, ou seja, você precisa de uma maneira de ir de uma superfície convexa para uma plana. O problema é que uma superfície convexa (imagine que algum tipo de desenho ou grade de coordenadas seja aplicada a ela) não pode ser transferida para uma superfície plana sem qualquer distorção. Pode estar distorcido: forma (ângulos), área, dimensões lineares. Existe, é claro, a possibilidade, e não a única, de se transferir para uma superfície plana sem distorcer apenas uma coisa, mas o resto será inevitavelmente distorcido.

Obviamente, em escalas menores (todo o planeta, continente), as distorções são mais pronunciadas do que nas maiores (região, cidade), e nas maiores (plano de uma pequena área), não estamos falando sobre elas, porque a superfície do elipsóide nessas escalas não é mais distinguível do plano.
Aqui chegamos ao conceito de projeção de mapa. Não vou dar exemplos com imagens de projeção de uma esfera em um cilindro ou cone com posterior desenvolvimento, para nós agora é importante generalizar.
A projeção de um mapa é uma forma matematicamente definida de mapear a superfície de um elipsóide em um plano. Simplificando, essas são fórmulas matemáticas para converter coordenadas geodésicas (latitude, longitude) em coordenadas cartesianas planas - exatamente o que precisamos.
Uma grande variedade de projeções cartográficas foi inventada, todas elas encontram aplicação para seus próprios fins: tamanhos iguais (em que a área é preservada) para mapas políticos, mapas de solo, conformados (em que a forma é preservada) - para navegação, equidistantes (mantendo a escala na direção escolhida) - para computador processamento e armazenamento de matrizes de geodados. Existem também projeções que combinam as características mencionadas acima em certas proporções, existem projeções completamente esotéricas. Exemplos de projeções podem ser encontrados na Wikipedia na página Lista de projeções de mapas.
Para cada projeção, ou fórmulas exatas são derivadas, ou na forma de uma soma de séries convergentes infinitas e, neste último caso, pode até haver várias opções. As fórmulas de projeção convertem latitude e longitude em coordenadas cartesianas, normalmente o metro é usado como unidade de medida em tais coordenadas. Essa grade retangular de 1 metro às vezes pode ser traçada em um mapa (na forma de uma grade de quilômetros), mas na maioria dos casos não é traçada. Mas agora sabemos que ele ainda existe de forma implícita. A escala do mapa, que é indicada no mapa, é calculada apenas em relação ao tamanho desta grade. Deve ser claramente entendido que 1 metro em tal grade de coordenadas corresponde a 1 metro no solo não em todo o mapa, mas apenas em alguns pontos, ao longo de uma certa linha, ou ao longo de linhas em uma determinada direção,em outros pontos ou direções aparecem distorções e 1 metro no solo não corresponde a 1 metro da grade de coordenadas.
Uma pequena digressão. A função da primeira parte do artigo WGS84_XYZ é precisamente a transformação de coordenadas do WGS84 SC em coordenadas retangulares, mas expressas não em metros, mas em pixels. Como você pode ver, a fórmula lá é extremamente simples, se a projeção de Mercator não fosse usada em uma esfera, mas em um elipsóide, a fórmula seria mais complicada. Por isso, para facilitar a vida dos navegadores, foi escolhida uma esfera, posteriormente a projeção do WebMercator se enraizou com sua esfera, pela qual é frequentemente criticada.
Para minhas próprias necessidades, usei um documento chamado "Projeções de mapas usadas pelo US Geological Survey" em formato pdf, que pode ser encontrado na Internet. O documento fornece instruções detalhadas para cada projeção para tornar mais fácil escrever uma função de transformação em uma linguagem de programação.
Vamos continuar escrevendo nosso algoritmo. Vamos implementar uma das projeções populares chamada de Transverse Mercator e uma de suas variantes chamada de projeção de 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,
};
Uma característica da projeção transversal de Mercator é que ela é conforme, ou seja, a forma dos objetos no mapa e os ângulos são preservados, assim como o fato de que a escala é preservada ao longo de um meridiano central. A projeção é adequada para todo o globo, mas as distorções aumentam com a distância do meridiano central, então toda a superfície da Terra é cortada em faixas estreitas ao longo dos meridianos, que são chamados de zonas, para cada um dos quais seu próprio meridiano central é usado. Exemplos da implementação de tais projeções são a projeção de Gauss-Kruger e UTM, em que o método de distribuição de coordenadas sobre zonas também é definido, ou seja, definido e com o mesmo nome SC.

E, de fato, o código das funções de inicialização e transformação de coordenadas. A função de inicialização avalia constantes uma vez que serão reutilizadas pela função de conversão.
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;
}
Na saída da função de transformação, teremos as coordenadas: os deslocamentos leste e norte (e, n) são coordenadas cartesianas retangulares em metros.

Já estamos muito perto de completar nosso algoritmo. Só temos que encontrar as coordenadas do pixel (x, y) no arquivo de mapa. Porque coordenadas de pixel também são cartesianas, então podemos encontrá-los por transformação afim (e, n) para (x, y). Voltaremos a encontrar os parâmetros da transformação mais afim um pouco mais tarde.
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;
}
E, finalmente, o algoritmo completo de criação de blocos:
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));
}
}
}
O resultado do trabalho para z = 12, y = 1377, x = 2391:

No algoritmo, a função de encontrar a imagem original do mapa srcimg a partir das coordenadas geodésicas lat, lon especificadas no SC do mapa permaneceu não descrita. Acho que não haverá problemas com isso e com o número da zona srcimg-> zona, mas nos deteremos em encontrar os parâmetros da transformação afim srcimg-> afim com mais detalhes.
Em algum lugar da Internet, muito tempo atrás, essa função foi encontrada para encontrar os parâmetros de uma transformação afim, eu cito até com comentários originais:
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;
}
}
Na entrada, pelo menos três pontos de ancoragem devem ser enviados, na saída temos parâmetros prontos. Os pontos de ancoragem são pontos para os quais as coordenadas de pixel (x, y) e as coordenadas de deslocamento leste e norte (e, n) são conhecidas. Os pontos de interseção da grade de quilômetros no mapa original podem ser usados como tais pontos. E se não houver uma grade de quilômetros no mapa? Em seguida, você pode definir os pares (x, y) e (lon, lat), já que esses pontos tomam os pontos de interseção dos paralelos e meridianos, eles estão sempre no mapa. Resta apenas converter (lon, lat) em (e, n), isso é feito pela função de transformação para a projeção, no nosso caso é translateTransverseMercator ().
Como você pode ver, os pontos de ancoragem são necessários para informar ao algoritmo qual parte da superfície da Terra o arquivo com a imagem do mapa descreve. Como os dois sistemas de coordenadas eram cartesianos, não importa quantos pontos de ancoragem definimos e não importa a que distância eles estejam um do outro, a discrepância em todo o plano do mapa estará dentro do erro na determinação dos pontos de ancoragem. A maioria dos erros é que a projeção, datum ou elipsóide errado (com parâmetros não especificados precisamente) é usado, como resultado, as coordenadas (e, n) na saída não estão no sistema de coordenadas cartesianas, mas em um ligeiramente curvado em relação ao cartesiano. Visualmente, isso pode ser visualizado como uma “folha amassada”. Naturalmente, aumentar o número de pontos de ancoragem não resolve esse problema. O problema pode ser resolvido ajustando os parâmetros da projeção, datum e elipsóide,neste caso, um grande número de pontos de ancoragem permitirá que você alise a “folha” com mais precisão e não perca as áreas não suavizadas.
E, no final do artigo, direi como usar tiles prontos em OpenLayers e de que forma prepará-los para o programa OsmAnd.
Para OpenLayers, os blocos simplesmente precisam ser postados na web e nomeados de forma que (z, y, x) possam ser destacados no nome do arquivo ou diretório, por exemplo:
/tiles/topo/12_1377_2391.jpg
ou, ainda melhor:
/ tiles / topo /12/1377/2391.jpg
Então eles podem ser usados assim:
map = new OpenLayers.Map (...);
map.addLayer(new OpenLayers.Layer.XYZ("Topo Maps", "/tiles/topo/${z}_${y}_${x}.jpg", {
isBaseLayer: false, visibility: false,
}));
Para o programa OsmAnd, é fácil determinar o formato de qualquer arquivo existente com um conjunto de blocos. Conto imediatamente sobre os resultados. Os blocos são compactados em um arquivo de banco de dados sqlite criado assim:
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
A coluna s sempre é preenchida com zero, para o qual, eu não entendi, a imagem é inserida na forma binária original, o formato (jpg, png, gif) é perdido, mas OsmAnd determina pelo seu conteúdo. Ladrilhos em diferentes formatos podem ser compactados em um banco de dados. Em vez de $ minzoom e $ maxzoom, é necessário substituir os limites de escala dos tiles inseridos na base. O arquivo de banco de dados completo, por exemplo, Topo.sqlitedb, é transferido para um smartphone ou tablet no diretório osmand / tiles. Reinicie o OsmAnd, em Menu -> "Configurar Mapa" -> "Camada Superior" uma nova opção "Topo" aparecerá - este é um mapa com nossos novos blocos.