В статье расскажу о результатах своих исследований, составим алгоритм преобразования произвольной растровой карты в тайлы, понятные для приложений и попутно познакомимся с такими понятиями как эллипсоид, датум, система координат, проекция.
Наша Земля имеет не форму шара, и даже не форму эллипсоида, эта сложная фигура называется геоид. Дело в том, что массы внутри Земли распределены не равномерно, поэтому в одних местах Земля немного вогнутая, в других немного выпуклая. Если взять территорию отдельной страны или материка, то ее поверхность с достаточной точностью можно совместить с поверхностью некоторого эллипсоида, если центр этого эллипсоида немного сдвинуть по трем осям координат относительно центра масс Земли. Такой эллипсоид называется референц-эллипсоидом, он пригоден для описания только локального участка Земли, для которого он был создан. На больших расстояниях от этого участка, он может иметь очень большое расхождение с поверхностью Земли. Эллипсоид, центр которого совпадает с центром масс Земли, называется общеземным эллипсоидом. Понятно, что референц-эллипсоид лучше описывает свой локальный участок Земли чем общеземной, но зато общеземной пригоден для всей поверхности Земли.
Для описания эллипсоида достаточно только двух независимых значений: экваториального радиуса (обычно обозначается a) и полярного радиуса (b), но вместо второго независимого значения обычно пользуются полярным сжатием f=(a-b)/a. Это первое, что нам понадобится в нашем алгоритме как объект эллипсоид. Для разных участков Земли в разные годы разными исследователями было вычислено множество референц-эллипсоидов, информация о них приводится в виде данных: a (в метрах) и 1/f (безразмерная). Как это ни странно, для общеземного эллипсоида также существует множество отличающихся вариантов (разные a,f), но отличие не очень сильное, связано оно в основном с различием в методиках определения a и 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,};
В примере приведены два эллипсодида: общеземной WGS84, используемой в спутниковой навигации, и референц-эллипсоид Красовского, применимый для территории Европы и Азии.
Рассмотрим еще один интересный момент, дело в том, что форма Земли медлено, но меняется, поэтому эллипсоид, который сегодня хорошо описывает поверхнось, через сотню лет может быть далек от реальности. Это мало касается общеземного эллипсоида, т.к. отклонения в пределах его же погрешности, но актуально для референц-эллипсоида. Тут мы подошли еще к одному понятию датум. Датум это совокупность параметров эллипсоида (a,f), его смещения внутри Земли (для референц-эллипсоида), зафиксированные в определенный момент времени. Если говорить более точно, то датум может описывать не обязательно эллипсоид, это могут быть параметры более сложной фигуры, например, квазигеоида.
Наверняка уже возник вопрос: как переходить от одного эллипсоида или датума к другому? Для этого на каждом эллипсоиде должна быть система геодезических координат: широта и долгота (фи, лямбда), переход осуществляется переводом координат из одной системы координат в другую.
Для преобразования координат существуют различные формулы. Можно сначала геодезичесике координаты в одной системе координат переводить в трехмерные координаты X,Y,Z, с ними выполнять сдвиги и повороты и затем полученные трехмерные координаты переводить в геодезические в другой системе координат. Можно это делать и напрямую. Т.к. все формулы это бесконечные сходящиеся ряды, то обычно ограничиваются несколькими членами рядов для достижения требуемой точности. В качестве примера воспользуемся преобразованиями Гельмерта (Helmert), это преобразования с использование перехода в трехмерные координаты, состоят из трех этапов описанных выше. Для преобразований кроме двух эллипсоидов понадобятся еще 7 параметров: три сдвига по трем осям, три угла поворота и масштабный коэффициент. Как можно догадаться, все параметры можно извлечь из датумов. Но в алгоритме мы не будем пользоваться таким объектом как датум, а вместо этого введем объект перехода из одной системы координат в другую, который будет содержать: ссылки на два эллипсоида и 7 параметров преобразования. Это будет вторым объектом нашего алгоритма.
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},};
В примере приведены параметры для преобразования из системы координат Пулково 1942 в систему координат WGS84. Сами параметры преобразования это отдельная тема, для некоторых систем координат они открыты, для других подобраны опытным путем, поэтому в разных источниках их значения могут незначительно отличаться.
Кроме параметров необходима и функция преобразования, она может быть прямая и для преобразования в обратном направлении, нам понадобится только преобразование в обратном направлении. Пропущу тонны матана и приведу свою оптимизированную функцию.
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;}
Откуда все это берется? На более понятном языке формулы можно найти в проекте proj4, но т.к. мне была необходима оптимизация по скорости выполнения, то всякие вычисления синуса суммы углов были преобразованы по формулам, возведения в степень оптимизированы венесениями за скобки, а константы посчитаны отдельно.
Теперь, чтобы приблизиться к выполнению изначальной задачи создать тайлы, необходимо рассмотреть систему координат под названием WebMercator. Эта система координат используется в приложении OsmAnd и в web, например в картах от Google и в OpenStreetMap. WebMercator это проекция Меркатора, построенная на сфере. Координаты в этой проекции это координаты пикселя X,Y и они зависят от масштаба Z, для нулевого масштаба вся земная поверхность (примерно до 85 градуса широты) помещается на одном тайле 256x256 пикселей, координаты X,Y меняются от 0 до 255, начиная с левого верхнего угла, для масштаба 1 уже 4 тайла, X,Y от 0 до 511 и так далее.
Для преобразования из WebMercator в 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);}
И под конец первой части статьи мы уже сможем набросать начало нашего алгоритма создания тайла. Каждый тайл 256x256 пикселей адресуется тремя значениями: x,y,z, соотношение с координатами X,Y,Z очень простое: 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 */ } }}
Координаты в СК42 это уже преобразованные координаты в систему координат карты, теперь осталось найти на карте пиксель по этим координатам и занести его цвет в пиксель тайла по координатам j,i. Об этом будет следующая статья, в которой мы поговорим о геодезических проекциях и афинных преобразованиях.
