[Map Projection] mocato projection (3)

This article mainly describes how to use different languages to realize the normal calculation and Inverse Calculation of mocato, that is, the conversion of latitude and longitude and plane coordinates.

Due to the rush of writing, there is something I don't understand in this article. I will add comments in turn in a few days.

 

Macato

 

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 mathdef 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));}

Elliptical mocato

 

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, WGS84r_minor:6356752.314245179,//defined as constantf:298.257223563,//1/f=(a-b)/a , a=r_major, b=r_minordeg2rad: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 radvar 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;}});//usagevar mercator=Conv.ll2m(47.6035525, 9.770602);//output mercator.x, mercator.yvar 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#endifinline 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;    }}

 

References: http://wiki.openstreetmap.org/wiki/Mercator