【地图投影】墨卡托投影（三）

本次主要说说用不同语言来实现墨卡托的正算和反算，即经纬度和平面坐标的相互转换。

由于编写仓促，文中有不明白的地方，过几天我会依次增加注释。

 

正球墨卡托

 

   JavaScript

 

function y2lat(a) { return 180/Math.PI * (2 * Math.atan(Math.exp(a*Math.PI/180)) - Math.PI/2); }
function lat2y(a) { return 180/Math.PI * Math.log(Math.tan(Math.PI/4+a*(Math.PI/180)/2)); }

C

#include <math.h>
#define deg2rad(d) ((d*M_PI)/180)
#define rad2deg(d) ((d*180)/M_PI)
#define earth_radius 6378137
 
/* The following functions take or return there results in degrees */
 
double y2lat_d(double y) { return rad2deg(2 * atan(exp(  deg2rad(y) ) ) - M_PI/2); }
double x2lon_d(double x) { return x; }
double lat2y_d(double lat) { return rad2deg(log(tan(M_PI/4+ deg2rad(lat)/2))); }
double lon2x_d(double lon) { return lon; }
 
/* The following functions take or return there results in something close to meters, along the equator */
 
double y2lat_m(double y) { return rad2deg(2 * atan(exp( (y / earth_radius ) ) - M_PI/2)); }
double x2lon_m(double x) { return rad2deg(x / earth_radius); }
double lat2y_m(double lat) { return earth_radius * log(tan(M_PI/4+ deg2rad(lat)/2)); }
double lon2x_m(double lon) { return deg2rad(lon) * earth_radius; }

PostGIS / SQL

INSERT INTO spatial_ref_sys (srid, auth_name, auth_srid, srtext, proj4text) VALUES 
 (900913,'EPSG',900913,'PROJCS["WGS84 / Simple Mercator",GEOGCS["WGS 84",
 DATUM["WGS_1984",SPHEROID["WGS_1984", 6378137.0, 298.257223563]],PRIMEM["Greenwich", 0.0],
 UNIT["degree", 0.017453292519943295],AXIS["Longitude", EAST],AXIS["Latitude", NORTH]],
 PROJECTION["Mercator_1SP_Google"],PARAMETER["latitude_of_origin", 0.0],
 PARAMETER["central_meridian", 0.0],PARAMETER["scale_factor", 1.0],PARAMETER["false_easting", 0.0],
 PARAMETER["false_northing", 0.0],UNIT["m", 1.0],AXIS["x", EAST],
 AXIS["y", NORTH],AUTHORITY["EPSG","900913"]]',
 '+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs');

Java

import java.lang.Math;
 
public class SphericalMercator {
  public static double y2lat(double aY) {
    return Math.toDegrees(2* Math.atan(Math.exp(Math.toRadians(aY))) - Math.PI/2);
  }
 
  public static double lat2y(double aLat) {
    return Math.toDegrees(Math.log(Math.tan(Math.PI/4+Math.toRadians(aLat)/2)));
  }
}

PHP

function lon2x($lon) { return deg2rad($lon) * 6378137.0; }
function lat2y($lat) { return log(tan(M_PI_4 + deg2rad($lat) / 2.0)) * 6378137.0; }
function x2lon($x) { return rad2deg($x / 6378137.0); }
function y2lat($y) { return rad2deg(2.0 * atan(exp($y / 6378137.0)) - M_PI_2); }

Python

import math
def y2lat(a):
  return 180.0/math.pi*(2.0*math.atan(math.exp(a*math.pi/180.0))-math.pi/2.0)
def lat2y(a):
  return 180.0/math.pi*math.log(math.tan(math.pi/4.0+a*(math.pi/180.0)/2.0))

C#

public double YToLatitude(double y)
		{
			return 180.0/System.Math.PI * 
				(2 * 
				 System.Math.Atan(
					System.Math.Exp(y*System.Math.PI/180)) - System.Math.PI/2);
		}
		public double LatitudeToY (double latitude)
		{
			return 180.0/System.Math.PI * 
				System.Math.Log(
					System.Math.Tan(
						System.Math.PI/4.0+latitude*(System.Math.PI/180.0)/2));
		}

椭球墨卡托

 

JavaScript

function deg_rad(ang) {
    return ang * (Math.PI/180.0)
}
function merc_x(lon) {
    var r_major = 6378137.000;
    return r_major * deg_rad(lon);
}
function merc_y(lat) {
    if (lat > 89.5)
        lat = 89.5;
    if (lat < -89.5)
        lat = -89.5;
    var r_major = 6378137.000;
    var r_minor = 6356752.3142;
    var temp = r_minor / r_major;
    var es = 1.0 - (temp * temp);
    var eccent = Math.sqrt(es);
    var phi = deg_rad(lat);
    var sinphi = Math.sin(phi);
    var con = eccent * sinphi;
    var com = .5 * eccent;
    con = Math.pow((1.0-con)/(1.0+con), com);
    var ts = Math.tan(.5 * (Math.PI*0.5 - phi))/con;
    var y = 0 - r_major * Math.log(ts);
    return y;
}
function merc(x,y) {
    return [merc_x(x),merc_y(y)];
}

 

var Conv=({
	r_major:6378137.0,//Equatorial Radius, WGS84
	r_minor:6356752.314245179,//defined as constant
	f:298.257223563,//1/f=(a-b)/a , a=r_major, b=r_minor
	deg2rad:function(d)
	{
		var r=d*(Math.PI/180.0);
		return r;
	},
	rad2deg:function(r)
	{
		var d=r/(Math.PI/180.0);
		return d;
	},
	ll2m:function(lon,lat) //lat lon to mercator
	{
		//lat, lon in rad
		var x=this.r_major * this.deg2rad(lon);
 
		if (lat > 89.5) lat = 89.5;
		if (lat < -89.5) lat = -89.5;
 
 
		var temp = this.r_minor / this.r_major;
		var es = 1.0 - (temp * temp);
		var eccent = Math.sqrt(es);
 
		var phi = this.deg2rad(lat);
 
		var sinphi = Math.sin(phi);
 
		var con = eccent * sinphi;
		var com = .5 * eccent;
		var con2 = Math.pow((1.0-con)/(1.0+con), com);
		var ts = Math.tan(.5 * (Math.PI*0.5 - phi))/con2;
		var y = 0 - this.r_major * Math.log(ts);
		var ret={'x':x,'y':y};
		return ret;
	},
	m2ll:function(x,y) //mercator to lat lon
	{
		var lon=this.rad2deg((x/this.r_major));
 
		var temp = this.r_minor / this.r_major;
		var e = Math.sqrt(1.0 - (temp * temp));
		var lat=this.rad2deg(this.pj_phi2( Math.exp( 0-(y/this.r_major)), e));
 
		var ret={'lon':lon,'lat':lat};
		return ret;
	},
	pj_phi2:function(ts, e) 
	{
		var N_ITER=15;
		var HALFPI=Math.PI/2;
 
 
		var TOL=0.0000000001;
		var eccnth, Phi, con, dphi;
		var i;
		var eccnth = .5 * e;
		Phi = HALFPI - 2. * Math.atan (ts);
		i = N_ITER;
		do 
		{
			con = e * Math.sin (Phi);
			dphi = HALFPI - 2. * Math.atan (ts * Math.pow((1. - con) / (1. + con), eccnth)) - Phi;
			Phi += dphi;
 
		} 
		while ( Math.abs(dphi)>TOL && --i);
		return Phi;
	}
});
//usage
var mercator=Conv.ll2m(47.6035525, 9.770602);//output mercator.x, mercator.y
var latlon=Conv.m2ll(5299424.36041, 1085840.05328);//output latlon.lat, latlon.lon

C

#include <math.h>
 
/*
 * Mercator transformation
 * accounts for the fact that the earth is not a sphere, but a spheroid
 */
#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;
}
 
double merc_x (double lon) {
        return R_MAJOR * deg_rad (lon);
}
 
double merc_y (double lat) {
        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;
        return 0 - R_MAJOR * log(ts);
}
 
static double rad_deg (double ang) {
        return ang * R_D;
}
 
double merc_lon (double x) {
        return rad_deg(x) / R_MAJOR;
}
 
double merc_lat (double y) {
        double ts = exp ( -y / R_MAJOR);
        double phi = M_PI_2 - 2 * atan(ts);
        double dphi = 1.0;
        int i;
        for (i = 0; fabs(dphi) > 0.000000001 && i < 15; i++) {
                double con = ECCENT * sin (phi);
                dphi = M_PI_2 - 2 * atan (ts * pow((1.0 - con) / (1.0 + con), COM)) - phi;
                phi += dphi;
        }
        return rad_deg (phi);
}

// Add this line before including math.h:
#define _USE_MATH_DEFINES
// Additions for MS Windows compilation:
#ifndef M_PI
	#define M_PI acos(-1.0)
#endif
#ifndef M_PI_2
	#define M_PI_2 1.57079632679489661922
#endif
inline double fmin(double x, double y) { return(x < y ? x : y); }
inline double fmax(double x, double y) { return(x > y ? x : y); }

C#

using System;
 
public static class MercatorProjection
{
    private static readonly double R_MAJOR = 6378137.0;
    private static readonly double R_MINOR = 6356752.3142;
    private static readonly double RATIO = R_MINOR / R_MAJOR;
    private static readonly double ECCENT = Math.Sqrt(1.0 - (RATIO * RATIO));
    private static readonly double COM = 0.5 * ECCENT;
 
    private static readonly double DEG2RAD = Math.PI / 180.0;
    private static readonly double RAD2Deg = 180.0 / Math.PI;
    private static readonly double PI_2 = Math.PI / 2.0;
 
    public static double[] toPixel(double lon, double lat)
    {
        return new double[] { lonToX(lon), latToY(lat) };
    }
 
    public static double[] toGeoCoord(double x, double y)
    {
        return new double[] { xToLon(x), yToLat(y) };
    }
 
    public static double lonToX(double lon)
    {
        return R_MAJOR * DegToRad(lon);
    }
 
    public static double latToY(double lat)
    {
        lat = Math.Min(89.5, Math.Max(lat, -89.5));
        double phi = DegToRad(lat);
        double sinphi = Math.Sin(phi);
        double con = ECCENT * sinphi;
        con = Math.Pow(((1.0 - con) / (1.0 + con)), COM);
        double ts = Math.Tan(0.5 * ((Math.PI * 0.5) - phi)) / con;
        return 0 - R_MAJOR * Math.Log(ts);
    }
 
    public static double xToLon(double x)
    {
        return RadToDeg(x) / R_MAJOR;
    }
 
    public static double yToLat(double y)
    {
        double ts = Math.Exp(-y / R_MAJOR);
        double phi = PI_2 - 2 * Math.Atan(ts);
        double dphi = 1.0;
        int i = 0;
        while ((Math.Abs(dphi) > 0.000000001) && (i < 15))
        {
            double con = ECCENT * Math.Sin(phi);
            dphi = PI_2 - 2 * Math.Atan(ts * Math.Pow((1.0 - con) / (1.0 + con), COM)) - phi;
            phi += dphi;
            i++;
        }
        return RadToDeg(phi);
    }
 
    private static double RadToDeg(double rad)
    {
        return rad * RAD2Deg;
    }
 
    private static double DegToRad(double deg)
    {
        return deg * DEG2RAD;
    }
}

PHP

function merc_x($lon)
{
	$r_major = 6378137.000;
	return $r_major * deg2rad($lon);
}
 
function merc_y($lat)
{
	if ($lat > 89.5) $lat = 89.5;
	if ($lat < -89.5) $lat = -89.5;
	$r_major = 6378137.000;
    $r_minor = 6356752.3142;
    $temp = $r_minor / $r_major;
	$es = 1.0 - ($temp * $temp);
    $eccent = sqrt($es);
    $phi = deg2rad($lat);
    $sinphi = sin($phi);
    $con = $eccent * $sinphi;
    $com = 0.5 * $eccent;
	$con = pow((1.0-$con)/(1.0+$con), $com);
	$ts = tan(0.5 * ((M_PI*0.5) - $phi))/$con;
    $y = - $r_major * log($ts);
    return $y;
}
 
function merc($x,$y) {
    return array('x'=>merc_x($x),'y'=>merc_y($y));
}
 
$array = merc(122,11); 
Java Implementation 

Java Implementation by Moshe Sayag, based on the JavaScript code published above, 17:11, 15.1.2008
 

public class Mercator {
    final private static double R_MAJOR = 6378137.0;
    final private static double R_MINOR = 6356752.3142;
 
    public double[] merc(double x, double y) {
        return new double[] {mercX(x), mercY(y)};
    }
 
    private double  mercX(double lon) {
        return R_MAJOR * Math.toRadians(lon);
    }
 
    private double mercY(double lat) {
        if (lat > 89.5) {
            lat = 89.5;
        }
        if (lat < -89.5) {
            lat = -89.5;
        }
        double temp = R_MINOR / R_MAJOR;
        double es = 1.0 - (temp * temp);
        double eccent = Math.sqrt(es);
        double phi = Math.toRadians(lat);
        double sinphi = Math.sin(phi);
        double con = eccent * sinphi;
        double com = 0.5 * eccent;
        con = Math.pow(((1.0-con)/(1.0+con)), com);
        double ts = Math.tan(0.5 * ((Math.PI*0.5) - phi))/con;
        double y = 0 - R_MAJOR * Math.log(ts);
        return y;
    }
}

Python

import math
 
def merc_x(lon):
  r_major=6378137.000
  return r_major*math.radians(lon)
 
def merc_y(lat):
  if lat>89.5:lat=89.5
  if lat<-89.5:lat=-89.5
  r_major=6378137.000
  r_minor=6356752.3142
  temp=r_minor/r_major
  eccent=math.sqrt(1-temp**2)
  phi=math.radians(lat)
  sinphi=math.sin(phi)
  con=eccent*sinphi
  com=eccent/2
  con=((1.0-con)/(1.0+con))**com
  ts=math.tan((math.pi/2-phi)/2)/con
  y=0-r_major*math.log(ts)
  return y

Java

public class Mercator {
    final private static double R_MAJOR = 6378137.0;
    final private static double R_MINOR = 6356752.3142;
 
    public double[] merc(double x, double y) {
        return new double[] {mercX(x), mercY(y)};
    }
 
    private double  mercX(double lon) {
        return R_MAJOR * Math.toRadians(lon);
    }
 
    private double mercY(double lat) {
        if (lat > 89.5) {
            lat = 89.5;
        }
        if (lat < -89.5) {
            lat = -89.5;
        }
        double temp = R_MINOR / R_MAJOR;
        double es = 1.0 - (temp * temp);
        double eccent = Math.sqrt(es);
        double phi = Math.toRadians(lat);
        double sinphi = Math.sin(phi);
        double con = eccent * sinphi;
        double com = 0.5 * eccent;
        con = Math.pow(((1.0-con)/(1.0+con)), com);
        double ts = Math.tan(0.5 * ((Math.PI*0.5) - phi))/con;
        double y = 0 - R_MAJOR * Math.log(ts);
        return y;
    }
}

 

参考资料：http://wiki.openstreetmap.org/wiki/Mercator