[Программирование, OpenStreetMap, C, Геоинформационные сервисы] Создание тайлов из растровых карт
Автор
Сообщение
news_bot ®
Стаж: 6 лет 9 месяцев
Сообщений: 27286
Как-то я озадачился вопросом создания карт, пригодных для использования в OsmAnd и OpenLayers. О ГИС я тогда вообще не имел ни малейшего понятия, поэтому разбирался со всем с нуля.
В статье расскажу о результатах своих «исследований», составим алгоритм преобразования произвольной растровой карты в тайлы, понятные для приложений и попутно познакомимся с такими понятиями как эллипсоид, датум, система координат, проекция.
Наша Земля имеет не форму шара, и даже не форму эллипсоида, эта сложная фигура называется геоид. Дело в том, что массы внутри Земли распределены не равномерно, поэтому в одних местах Земля немного вогнутая, в других немного выпуклая. Если взять территорию отдельной страны или материка, то ее поверхность с достаточной точностью можно совместить с поверхностью некоторого эллипсоида, если центр этого эллипсоида немного сдвинуть по трем осям координат относительно центра масс Земли. Такой эллипсоид называется референц-эллипсоидом, он пригоден для описания только локального участка Земли, для которого он был создан. На больших расстояниях от этого участка, он может иметь очень большое расхождение с поверхностью Земли. Эллипсоид, центр которого совпадает с центром масс Земли, называется общеземным эллипсоидом. Понятно, что референц-эллипсоид лучше описывает свой локальный участок Земли чем общеземной, но зато общеземной пригоден для всей поверхности Земли.
Для описания эллипсоида достаточно только двух независимых значений: экваториального радиуса (обычно обозначается 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. Об этом будет следующая статья, в которой мы поговорим о геодезических проекциях и афинных преобразованиях.
===========
Источник:
habr.com
===========
Похожие новости:
- [Разработка под Android, Хакатоны, Kotlin, Искусственный интеллект] Челлендж по разговорному ИИ на хакатоне Junction: создай чатбота или голосовой навык и выиграй 10 000 евро
- [Разработка под iOS, Управление персоналом, Карьера в IT-индустрии] Сценарий идеального технического собеседования
- [Настройка Linux, Программирование, Разработка под Linux] Знакомимся с семафорами в Linux (перевод)
- [Гаджеты, Компьютерное железо, Искусственный интеллект, Настольные компьютеры] Nvidia представила новую версию одноплатного ПК Jetson Nano всего за $59
- [Алгоритмы, Машинное обучение, Искусственный интеллект] Ученые разработали метод обучения ИИ с меньшим числом параметров, который превзошел GPT-3
- [Космонавтика, Научно-популярное, Промышленное программирование, IT-компании] Определённо не Windows 95: какие операционные системы поддерживают работу в космосе? (перевод)
- [Accessibility, Дизайн, Здоровье] Доброшрифт: так пишутся добрые дела
- [Управление проектами, Управление продуктом, IT-компании] Прививка инноваций через инициативные проекты
- [Исследования и прогнозы в IT, Микросервисы, Программирование, Распределённые системы] Обработка распределенных транзакций в микросервисной архитектуре (перевод)
- [Go, Конференции] Технические доклады Lamoda на GolangLive 2020
Теги для поиска: #_programmirovanie (Программирование), #_openstreetmap, #_c, #_geoinformatsionnye_servisy (Геоинформационные сервисы), #_gis (гис), #_datum (датум), #_sistemy_koordinat (системы координат), #_proektsii (проекции), #_osmand, #_openlayers, #_openstreetmap, #_programmirovanie (
Программирование
), #_openstreetmap, #_c, #_geoinformatsionnye_servisy (
Геоинформационные сервисы
)
Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете голосовать в опросах
Вы не можете прикреплять файлы к сообщениям
Вы не можете скачивать файлы
Текущее время: 22-Ноя 23:21
Часовой пояс: UTC + 5
| Автор | Сообщение |
|---|---|
|
news_bot ®
Стаж: 6 лет 9 месяцев |
|
|
Как-то я озадачился вопросом создания карт, пригодных для использования в OsmAnd и OpenLayers. О ГИС я тогда вообще не имел ни малейшего понятия, поэтому разбирался со всем с нуля. В статье расскажу о результатах своих «исследований», составим алгоритм преобразования произвольной растровой карты в тайлы, понятные для приложений и попутно познакомимся с такими понятиями как эллипсоид, датум, система координат, проекция. Наша Земля имеет не форму шара, и даже не форму эллипсоида, эта сложная фигура называется геоид. Дело в том, что массы внутри Земли распределены не равномерно, поэтому в одних местах Земля немного вогнутая, в других немного выпуклая. Если взять территорию отдельной страны или материка, то ее поверхность с достаточной точностью можно совместить с поверхностью некоторого эллипсоида, если центр этого эллипсоида немного сдвинуть по трем осям координат относительно центра масс Земли. Такой эллипсоид называется референц-эллипсоидом, он пригоден для описания только локального участка Земли, для которого он был создан. На больших расстояниях от этого участка, он может иметь очень большое расхождение с поверхностью Земли. Эллипсоид, центр которого совпадает с центром масс Земли, называется общеземным эллипсоидом. Понятно, что референц-эллипсоид лучше описывает свой локальный участок Земли чем общеземной, но зато общеземной пригоден для всей поверхности Земли. Для описания эллипсоида достаточно только двух независимых значений: экваториального радиуса (обычно обозначается 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. Об этом будет следующая статья, в которой мы поговорим о геодезических проекциях и афинных преобразованиях. =========== Источник: habr.com =========== Похожие новости:
Программирование ), #_openstreetmap, #_c, #_geoinformatsionnye_servisy ( Геоинформационные сервисы ) |
|
Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете голосовать в опросах
Вы не можете прикреплять файлы к сообщениям
Вы не можете скачивать файлы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете голосовать в опросах
Вы не можете прикреплять файлы к сообщениям
Вы не можете скачивать файлы
Текущее время: 22-Ноя 23:21
Часовой пояс: UTC + 5