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