No artigo, vou falar sobre os resultados da minha "pesquisa", compor um algoritmo para converter um mapa raster arbitrário em blocos que são compreensíveis para aplicações e, ao longo do caminho, me familiarizarei com conceitos como elipsóide, datum, sistema de coordenadas, projeção.
Nossa Terra não tem a forma de uma bola, nem mesmo de um elipsóide, esta figura complexa é chamada de geóide. O fato é que as massas dentro da Terra não estão uniformemente distribuídas, então em alguns lugares a Terra é ligeiramente côncava, em outros é ligeiramente convexa. Se tomarmos o território de um país ou continente separado, então sua superfície com precisão suficiente pode ser alinhada com a superfície de algum elipsóide, se o centro deste elipsóide for ligeiramente deslocado ao longo de três eixos coordenados em relação ao centro de massa da Terra. Esse elipsóide é chamado de elipsóide de referência e é adequado para descrever apenas a área local da Terra para a qual foi criado. A grandes distâncias deste local, pode haver uma discrepância muito grande com a superfície da Terra. Um elipsóide cujo centro coincide com o centro de massa da Terra é denominado elipsóide terrestre comum. Claro,que o elipsóide de referência descreve melhor sua área local da Terra do que o terrestre geral, mas o terrestre geral é adequado para toda a superfície da Terra.
Para descrever o elipsóide, apenas dois valores independentes são suficientes: o raio equatorial (geralmente denotado por a) e o raio polar (b), mas em vez do segundo valor independente, a contração polar f = (ab) / a é normalmente usada. Esta é a primeira coisa que precisamos em nosso algoritmo como um objeto - um elipsóide. Para diferentes partes da Terra em diferentes anos, diferentes pesquisadores calcularam muitos elipsóides de referência, as informações sobre eles são fornecidas na forma de dados: a (em metros) e 1 / f (adimensionais). Curiosamente, para um elipsóide terrestre comum, também existem muitas variantes diferentes (diferentes a, f), mas a diferença não é muito forte, deve-se principalmente à diferença nos métodos para determinar a e 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,
};
O exemplo mostra dois elipsóides: o comum terrestre WGS84, usado na navegação por satélite, e o elipsóide de referência de Krasovsky, aplicável ao território da Europa e Ásia.
Considere outro ponto interessante, o fato é que a forma da Terra é lenta, mas está mudando, então o elipsóide, que hoje descreve bem a superfície, daqui a cem anos pode estar longe da realidade. Isso tem pouco a ver com o elipsóide terrestre comum, uma vez que desvios dentro de seu mesmo erro, mas relevantes para o elipsóide de referência. Aqui chegamos a outro conceito - datum. Datum é um conjunto de parâmetros de um elipsóide (a, f), seus deslocamentos no interior da Terra (para um elipsóide de referência), fixados em um determinado ponto no tempo. Mais precisamente, o datum pode não necessariamente descrever um elipsóide, podem ser parâmetros de uma figura mais complexa, por exemplo, um quase-geóide.
Certamente a pergunta já se levantou: como passar de um elipsóide ou datum para outro? Para fazer isso, cada elipsóide deve ter um sistema de coordenadas geodésicas: latitude e longitude (phi, lambda), a transição é realizada pela tradução de coordenadas de um sistema de coordenadas para outro.
Existem várias fórmulas para transformar as coordenadas. Você pode primeiro traduzir as coordenadas geodésicas em um sistema de coordenadas em coordenadas tridimensionais X, Y, Z, realizar mudanças e rotações com elas e, em seguida, converter as coordenadas tridimensionais resultantes em coordenadas geodésicas em outro sistema de coordenadas. Você pode fazer isso diretamente. Porque todas as fórmulas são séries convergentes infinitas, então geralmente são limitadas a alguns membros da série para alcançar a precisão necessária. Como exemplo, usaremos as transformadas de Helmert, são transformações que usam uma transição para coordenadas tridimensionais, consistem nos três estágios descritos acima. Para transformações, além de dois elipsóides, você precisa de mais 7 parâmetros: três deslocamentos ao longo de três eixos, três ângulos de rotação e um fator de escala. Como você pode imaginar, todos os parâmetros podem ser extraídos dos datums.Porém, no algoritmo, não usaremos tal objeto como datum, mas, em vez disso, introduziremos um objeto de transição de um sistema de coordenadas para outro, que conterá: referências a dois elipsóides e 7 parâmetros de transformação. Este será o segundo objeto de nosso 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},
};
O exemplo mostra os parâmetros para a conversão do sistema de coordenadas Pulkovo 1942 para o sistema de coordenadas WGS84. Os próprios parâmetros de transformação são um tópico separado, para alguns sistemas de coordenadas eles estão abertos, para outros eles são selecionados empiricamente, portanto, seus valores podem diferir ligeiramente em fontes diferentes.
Além dos parâmetros, também é necessária uma função de transformação, que pode ser direta e para transformação na direção oposta, precisamos apenas de uma transformação na direção oposta. Vou pular toneladas de matemática e dar minha função otimizada.
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 onde vem tudo isso? Em uma linguagem mais compreensível, as fórmulas podem ser encontradas no projeto proj4, mas desde Eu precisava de otimização para a velocidade de execução, então quaisquer cálculos do seno da soma dos ângulos eram transformados de acordo com as fórmulas, as exponenciações eram otimizadas por folheados entre parênteses e as constantes eram calculadas separadamente.
Agora, para chegar mais perto de completar a tarefa original de criar blocos, precisamos considerar um sistema de coordenadas chamado WebMercator. Este sistema de coordenadas é usado no aplicativo OsmAnd e na web, por exemplo, no Google maps e no OpenStreetMap. WebMercator é uma projeção de Mercator construída sobre uma esfera. As coordenadas nesta projeção são as coordenadas do pixel X, Y e dependem da escala Z, para uma escala zero toda a superfície da Terra (até cerca de 85 graus de latitude) é colocada em um bloco de 256x256 pixels, as coordenadas X, Y mudam de 0 a 255, começando da esquerda canto superior, para a escala 1 - já 4 peças, X, Y - de 0 a 511 e assim por diante.
As seguintes funções são usadas para converter de WebMercator para 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);
}
E no final da primeira parte do artigo, já podemos esboçar o início de nosso algoritmo para a criação de um ladrilho. Cada bloco de 256x256 pixels é endereçado por três valores: x, y, z, a relação com as coordenadas X, Y, Z é muito simples: 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 */
}
}
}
As coordenadas no SK42 já são coordenadas transformadas no sistema de coordenadas do mapa, agora resta encontrar um pixel no mapa por essas coordenadas e inserir sua cor em um pixel de ladrilho nas coordenadas j, i. Este será o próximo artigo, no qual falaremos sobre projeções geodésicas e transformações afins.