Пересчет координат из Lat/Long в проекцию Меркатора и обратно

Эта страница является черновиком статьи.

Описание вариантов пересчета широты/долготы в плоские координаты проекции Меркатора на базе эллипсоида WGS84 и на базе сфероида

Введение

Картографический веб-сервис Google Maps для отображения карт использует проекцию Меркатора на сфере. Для начала рассмотрим более общий случай проекции Меркатора для эллипсоида WGS84, а потом будет не сложно перейти к сфере. !Переписать!

Описание проекции Меркатора/WGS84

Равноугольная цилиндрическая проекция Меркатора — одна из основных картографических проекций. Разработана Герардом Меркатором для применения в его «Атласе» в 1596г.
В наше время достаточно часто используется вариант этой проекции WGS 84/World Mercator (EPSG:3395) на базе общеземного эллипсоида WGS84 (EPSG:7030). Данный эллипсоид имеет размер большой полуоси (экваториальный радиус) равный 6378137,0 метров и меньшей полуоси (полярный радиус) — ~6356752,314245 метра. Величина обратного уплощения (1/f) равна 298,257223563.
WKT представление это проекции:

PROJCS["WGS 84 / World Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.01745329251994328,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    AUTHORITY["EPSG","3395"],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH]]

Proj.4 представление:

+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs

Описание проекции "Cферического Меркатора"

В последнее время, в связи со стремительным развитием картографических веб сервисов, большое распространение получил другой вариант проекции Меркатора - на базе не эллипсоида, а сферы. Этот выбор обусловлен более простыми расчетами, которые могут быть быстро выполнены тонкими клиентами этих сервисов (например js скрипт в браузере). Часто, эту проекцию называют "сферическим Меркатором" (EPSG:3857).
Эта проекция имеет множество названий и кодов: OSGEO:41001, Google Mercator EPSG:900913, ESRI Web Mercator ESRI:102100 и ESRI:102113, Popular Visualisation CRS / Mercator EPSG:3785. Но наиболее правильно использовать код EPSG:3857, так как этим кодом были заменены все остальные, как в пространстве кодов EPSG, так и в ESRI.
Отличительное свойство этой проекции - равные полуоси сфероида, размером 6378137,0 метров.
WKT представление проекции:

PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["Popular Visualisation CRS",
        DATUM["Popular_Visualisation_Datum",
            SPHEROID["Popular Visualisation Sphere",6378137,0,
                AUTHORITY["EPSG","7059"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6055"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.01745329251994328,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4055"]],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    AUTHORITY["EPSG","3785"],
    AXIS["X",EAST],
    AXIS["Y",NORTH]]

Proj.4 представление:

+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

Пересчет координат из широты/долготы в проекцию Меркатора/WGS84

Теория

Формулы пересчета:
$X=a*long$
$Y=a*\ln[\tan({\frac {\pi }{4}}+{\frac {lat}{2}})*({\frac {1-e*\sin lat}{1+e*\sin lat}})^{e/2}]$
Где:

• long/lat - долгота/широта в радианах;
• e - эксцентриситет эллипса, $e={\sqrt {1-({\frac {b}{a}})^{2}}}$
• a - большая полуось эллипса;
• b - малая полуось эллипса;

Реализации

Для примера будем использовать координаты г. Москва: 55.751667 с.ш., 37.617778 в.д

Консольные команды Proj.4

Для пересчета используется команда cs2cs:

echo 37.617778 55.751667 |cs2cs +proj=latlong +ellps=WGS84 +to +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs

Результат:

4187591.89   7473789.46

Код на С\С++ для Proj.4

#include <iostream>
#include <iomanip>
#include <proj_api.h>

int LatLong2Merc(double *x, double *y)
{
  projPJ pjLatlong = pj_init_plus("+proj=longlat +ellps=WGS84");
  projPJ pjMerc = pj_init_plus("+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs");

  int result=pj_transform(pjLatlong, pjMerc, 1, 1, x, y, NULL );

  pj_free(pjLatlong);
  pj_free(pjMerc);

  return result;
}

int main(int argc, char **argv)
{
  double lon=37.617778f;
  double lat=55.751667f;

  std::cout<<"Long/Lat:"<<std::setprecision(9)<<lon<<"  "<<lat<<std::endl;
    
  lon*=DEG_TO_RAD;
  lat*=DEG_TO_RAD;

  LatLong2Merc(&lon, &lat);
  std::cout<<"X/Y: "<<lon<<"  "<<lat<<std::endl;
}

Результат:

4187591.98  7473789.47

Код на R

LatLongToMerc <- function(lon, lat)
{
  rLat = lat*pi/180;
  rLong = lon*pi/180;
  a=6378137.0;
  b=6356752.3142;
  f=(a-b)/a;
  e=sqrt(2*f-f^2);
  X=a*rLong;
  Y=a*log(tan(pi/4+rLat/2)*((1-e*sin(rLat))/(1+e*sin(rLat)))^(e/2));
  return(list('x'=X, 'y'=Y));
}

res<-LatLongToMerc(37.617778,55.751667)
print(res['x'], digits=9)
print(res['y'], digits=9)

Ссылки по теме

Mercator projection(eng wiki)
Варианты пересчета на вики OSM
Web Mercator: Non-Conformal, Non-Mercator
Tiles à la Google Maps: Coordinates, Tile Bounds and Projection
Spatial references