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

Материал из GIS-Lab
Перейти к навигации Перейти к поиску
Эта страница является черновиком статьи.


Описание вариантов пересчета широты/долготы в плоские координаты проекции Меркатора на базе эллипсоида 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

Теория

Формулы пересчета:


Где:

  • long/lat - долгота/широта в радианах;
  • e - эксцентриситет эллипса,
  • a - большая полуось эллипса;
  • b - малая полуось эллипса;

Реализации

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

MS Excel

Консольные команды 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

WTF???

Код на 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)

Результат: 

4187591.89 7473789.46

Код на Python

import math
   
def LatLongToMerc(lon, lat):
    if lat>89.5:
        lat=89.5
    if lat<-89.5:
        lat=-89.5
        
    rLat = math.radians(lat)
    rLong = math.radians(lon)
    
    a=6378137.0
    b=6356752.3142
    f=(a-b)/a
    e=math.sqrt(2*f-f**2)
    x=a*rLong
    y=a*math.log(math.tan(math.pi/4+rLat/2)*((1-e*math.sin(rLat))/(1+e*math.sin(rLat)))**(e/2))
    return {'x':x,'y':y}

res = LatLongToMerc(37.617778,55.751667)
print res['x'], res['y']

Результат: 

4187591.89173 7473789.4619

Код на С


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

Теория

Реализации


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

Теория

Из-за идентичности полуосей сферойда, формулы пересчета координат для сферического Меркатора проще:


Где:

  • long/lat - долгота/широта в радианах;
  • a - полуось эллипса;


Реализации

Консольные команды 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 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

Результат: 

4187591.89      7509137.58

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

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

int LatLong2SpherMerc(double *x, double *y)
{
  projPJ pjLatlong = pj_init_plus("+proj=longlat +ellps=WGS84");
  projPJ pjSpherMerc = pj_init_plus("+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");

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

  pj_free(pjLatlong);
  pj_free(pjSpherMerc);

  return result;
}

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

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

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

Результат: 

4187591.89  7509137.58


Пересчет координат из проекции "Сферического Меркатора" в широту/долготу

Теория

Реализации


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

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

Источник — https://wiki.gis-lab.info/index.php?title=Пересчет_координат_из_Lat/Long_в_проекцию_Меркатора_и_обратно&oldid=4712