В предыдущей части мы остановились на том, что необходимо получить цвет пикселя из исходной растровой карты по геодезическим координатам широта,
долгота), уже переведенным в СК этой карты.
Геодезические координаты заданы на поверхности эллипсоида, а координаты пикселя на плоскости, т.е. необходим способ перехода от выпуклой поверхности к плоской. Проблема состоит в том, что выпуклую поверхность (представим что на ней нанесен какой-то рисунок или координатная сетка), невозможно перенести на плоскую без каких-либо искажений. Искажаться может: форма (углы), площадь, линейные размеры. Есть, конечно, возможность, причем не единственная, переноса на плоскую поверхность не искажая только что-то одно, но неизбежно будет искажаться остальное.
Очевидно, что на более мелких масштабах (вся планета, континет), искажения проявляются сильнее, чем на более крупных (область, город), а на самых крупных (план небольшого участка) о них вообще речь не идет, т.к. поверхность эллипсоида на таких масштабах уже не отличима от плоскости.
Здесь мы подошли к понятию картографической проекции. Не буду приводить примеры с картинками проецирования сферы на цилиндр или на конус с последующей разверткой, для нас сейчас важно обобщить.
Картографическая проекция это заданный математически способ отображения поверхности эллипсоида на плоскость. Проще говоря, это математические формулы для перевода геодезических координат (широта, долгота) в плоские декартовы как раз то, что нам нужно.
Картографических проекций изобретено огромное множество, все они находят применение для своих целей: равновеликие (в кторых сохраняется площадь) для политических карт, карт грунтов, равноугольные (в которых сохраняется форма) для навигации, равнопромежуточные (сохранение масштаба по выбранному направлению) для компьютерной обработки и хранения массивов геоданных. Есть также проекции, которые сочетают в себе в определенных пропорциях вышеперечисленные особенности, есть совсем эзотерические проекции. Примеры проекций можно найти в Википедии на странице Список картографических проекций.
Для каждой проекции выведены либо точные формулы, либо в виде суммы бесконечных сходящихся рядов, причем в последнем случае может существовать даже несколько вариантов. Формулы для проекций переводят широту и долготу в декартовы координаты, как правило единицей измерения в таких координатах используется метр. Такая метровая прямоугольная координатная сетка иногда может быть нанесена и на карту (в виде километровой сетки), но в большинстве случаев ее не наносят. Но мы теперь знаем, что в неявном виде она все-таки есть. Масштаб карты, который указан на карте, какраз и вычислен относительно размера этой сетки. Следует четко понимать, что 1 метр на такой координатной сетке соответствует 1 метру на местности не на всей карте, а лишь в некоторых точках, по определенной линии, либо по линиям в определенном направлении, в остальных же точках или направлениях проявляются искажения и 1 метр на местности не соответствет 1 метру координатной сетки.
Небольшое отступление. Функция из первой части статьи WGS84_XYZ как раз и является преобразованием координат из СК WGS84 в прямоугольные координаты, но выраженные не в метрах, а в пикселях. Как видно, формула там предельно простая, если бы использовалась проекция Меркатора не на сфере, а на эллипсоиде, то формула была бы посложнее. Именно поэтому, чтобы облегчить жизнь браузерам, была выбрана сфера, в последствии проекция WebMercator так и прижилась со своей сферой, за что ее часто ругают.
Для своих нужд я использовал документ под названием Map projections used by the U.S. Geological Survey в формате pdf, который можно найти в Интернет. В документе для каждой проекции приведены подробные инструкции, по которым легко написать функцию преобразования на языке программирования.
Продолжим писать наш алгоритм. Реализуем одну из популярных проекций под названием поперечная проекция Меркатора (Transverse Mercator) и один из его вариантов под названием проекция Гаусса-Крюгера (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,};
Особенностью поперечной проекции Меркатора является то, что она конформна, т.е сохраняется форма объектов на карте и углы, а также то, что сохраняется масштаб вдоль одного центрального меридиана. Проекция пригодна для всего земного шара, но при отдалении от центрального меридиана растут искажения, поэтому всю земную поверхность нарезают на узкие полоски по меридианам, которые называются зонами, для каждой из которых используется свой центральный меридиан. Примерами реализации таких проекций являются проекция Гаусса-Крюгера и UTM, в которых также определен способ распределения координат по зонам, т.е. определена и одноименная СК.
И, собственно, код функций инициализации и преобразования координат. В функции инициализации одноразово вычисляются константы, которые будут многократно использоваться функцией преобразования.
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;}
На выходе функции преобразования будем иметь координаты: восточное и северное смещение (e,n) это прямоугольные декартовы координаты в метрах.
Мы уже очень близки к завершению нашего алгоритма. Нам осталось только найти координаты пикселя (x,y) в файле карты. Т.к. координаты пикселей тоже декартовы, то мы можем найти их афинным преобразованием (e,n) к (x,y). К нахождению параметров самого афинного преобразования мы вернемся чуть позже.
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;}
И, наконец, полный алгоритм создания тайла:
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)); } }}
Результат работы для z=12, y=1377, x=2391:
В алгоритме осталась не описанной функция нахождения исходного изображения карты srcimg по заданным в СК карты геодезическим координатам lat, lon. Думаю, с ней и номером зоны srcimg->zone проблем не возникнет, а вот на нахождении параметров афинного преобразования srcimg->affine остановимся более подробно.
Где-то на просторах Интернета очень давно была найдена вот такая функция для нахождения параметров афинного преобразования, привожу ее даже с оригинальными комментариями:
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; }}
На вход необходимо подать не менее трех точек привязки, на выходе получаем готовые параметры. Точки привязки это точки, для которых одновременно известны координаты пикселя (x,y) и координаты восточного и северного смещения (e,n). В качестве таких точек можно использовать точки пересечения километровой сетки на исходной карте. А что если на карте нет километровой сетки? Тогда можно задать пары (x,y) и (lon,lat), в качестве таких точек взять точки пересечения параллелей и меридианов, они на карте есть всегда. Осталось только преобразовать (lon,lat) в (e,n), это делается функцией преобразования для проекции, в нашем случае это translateTransverseMercator().
Как видно, точки привязки необходимы, чтобы сообщить алгоритму, какой участок земной поверхности описывает файл с изображением карты. Поскольку обе координатные системы были декартовыми, то сколько бы точек привязки мы не задавали и как далеко бы они не находились друг от друга, расхождение по всей плоскости карты будет в пределах погрешности определения точек привязки. Большинство ошибок состоит в том, что применяется не правильная (с не точно заданными параметрами) проекция, датум или эллипсоид, в результате координаты (e,n) на выходе получаются не в декартовой системе координат, а в немного искривленной относительно декартовой. Визуально это можно представить как мятую простыню. Естественно, увеличение количества точек привязки эту проблему не решает. Проблему может решить тюнинг параметров проекции, датума и эллипсоида, в данном случае большое количество точек привязки позволит разгладить простыню точнее и не пропустить неразглаженных участков.
И в завершение статьи расскажу, как готовые тайлы использовать в OpenLayers и в каком виде их приготовить для программы OsmAnd.
Для OpenLayers тайлы необходимо просто выложить в web и назвать так, чтобы в названии файла или каталога можно было выделить (z,y,x), например так:
/tiles/topo/12_1377_2391.jpg
или, еще лучше так:
/tiles/topo/12/1377/2391.jpg
Тогда их использовать можно таким образом:
map = new OpenLayers.Map (...);map.addLayer(new OpenLayers.Layer.XYZ("Topo Maps", "/tiles/topo/${z}_${y}_${x}.jpg", { isBaseLayer: false, visibility: false, }));
Для программы OsmAnd легко определить формат по каким-либо уже имеющимся файлам с набором тайлов. Расскажу сразу про результаты. Тайлы упаковываются в файл базы данных sqlite, который создается таким образом:
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
Колонка s всегда заполняется нулем, для чего она, я не разбирался, в image заносится картинка в исходном бинарном виде, формат (jpg, png, gif) теряется, но OsmAnd определяет его сам по содержимому. В одной базе могут быть упакованы тайлы в разных форматах. Вместо $minzoom и $maxzoom необходимо подставить пределы масштаба занесенных в базу тайлов. Заполненный файл базы данных, например, Topo.sqlitedb переносим на смартфон или планшет в каталог osmand/tiles. Перезапускаем OsmAnd, в Меню -> Configure Map -> Верхний слой появится новая опция Topo это карта с нашими новыми тайлами.
