Пересчет координат из 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 в.д

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.617778;
  double lat=55.751667;

  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.89  7473789.46

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

Код на С

Взято тут

#include <math.h>
#include <stdio.h>

#define D_R (M_PI / 180.0)
#define R_D (180.0 / M_PI)
#define R_MAJOR 6378137.0
#define R_MINOR 6356752.3142
#define RATIO (R_MINOR/R_MAJOR)
#define ECCENT (sqrt(1.0 - (RATIO * RATIO)))
#define COM (0.5 * ECCENT)
 
static double deg_rad (double ang) {
        return ang * D_R;
}
 
void LatLong2Merc(double lon, double lat, double* x, double* y) {
	*x =  R_MAJOR * deg_rad (lon);
	lat = fmin (89.5, fmax (lat, -89.5));
        double phi = deg_rad(lat);
        double sinphi = sin(phi);
        double con = ECCENT * sinphi;
        con = pow((1.0 - con) / (1.0 + con), COM);
        double ts = tan(0.5 * (M_PI * 0.5 - phi)) / con;
        *y = 0 - R_MAJOR * log(ts);
}

void main(){
  double x=0,y=0;
  LatLong2Merc(37.617778,55.751667,&x,&y);
  printf(""X: %10.2f Y: %10.2f\n"",x,y);
}

Результат:

4187591.89  7473789.46

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

Теория

Реализации

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

Теория

Из-за идентичности полуосей сферойда, формулы пересчета координат для сферического Меркатора проще:
$X=a*long$
$Y=a*\ln[\tan({\frac {\pi }{4}}+{\frac {lat}{2}})]$
Где:

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

Код на R

LatLongToSpherMerc <- function(lon, lat)
{
  rLat = lat*pi/180;
  rLong = lon*pi/180;
  a=6378137.0;
  X=a*rLong;
  Y=a*log(tan(pi/4+rLat/2));
  return(list('x'=X, 'y'=Y));
}
 
res<-LatLongToSpherMerc(37.617778,55.751667)
print(res['x'], digits=9)
print(res['y'], digits=9)

Результат:

4187591.89 7509137.58

Код на Python

import math
 
def LatLongSpherToMerc(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
    x=a*rLong
    y=a*math.log(math.tan(math.pi/4+rLat/2))
    return {'x':x,'y':y}
 
res = LatLongToSpherMerc(37.617778,55.751667)
print res['x'], res['y']

Результат:

4187591.89173 7509137.5811

Код на С

#include <math.h>
#include <stdio.h>

#define D_R (M_PI / 180.0)
#define R_MAJOR 6378137.0
 
static double deg_rad (double ang) {
        return ang * D_R;
}

void LatLong2SpherMerc(double lon, double lat, double* x, double* y) {
	lat = fmin (89.5, fmax (lat, -89.5));
	*x = R_MAJOR * deg_rad (lon);
	*y = R_MAJOR * log(tan(M_PI_4 + deg_rad(lat)/2 ));
}

void main(){
  double x=0,y=0;
  LatLong2SpherMerc(37.617778,55.751667,&x,&y);
  printf("X: %10.2f Y: %10.2f\n",x,y);
}

Результат:

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

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

Содержание

Введение

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

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

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

Теория

Реализации

MS Excel

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

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

Код на R

Код на Python

Код на С

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

Теория

Реализации

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

Теория

Реализации

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

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

Код на R

Код на Python

Код на С

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

Теория

Реализации

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

Навигация

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

Введение

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

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

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

Теория

Реализации

MS Excel

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

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

Код на R

Код на Python

Код на С

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

Теория

Реализации

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

Теория

Реализации

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

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

Код на R

Код на Python

Код на С

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

Теория

Реализации

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

Навигация

Поиск