経緯度とGoogle Mercatorとの変換式を確認してみる（OpenLayers & Proj4js）。

既存の JavaScript プログラムに埋め込まれた、経緯度とGoogle Mercatorとの変換式（特にメルカトルから経緯度への逆変換式）を確認してみました。Safari の Web インスペクタを使って簡単に調べてみます。

基本的方法

Safari や Chrome に付属している Web インスペクタには「コンソール」があります。ここに JavaScript コードを書き込むと、そのページでロードされている JavaScript オブジェクトに簡単にアクセスすることができます。Ruby の irb を使うように、Web インスペクタコンソールを使って、JavaScript プログラムを理解していきます。

OpenLayers の場合

http://trac.openlayers.org/wiki/SphericalMercator を読むと結局、OpenLayers.Layer.SphericalMercator.projectInverse が逆変換関数であるようです。調べてみると、実際には、中で別の関数が呼び出されていました。

> OpenLayers.Layer.SphericalMercator.projectInverse
function (point) {var lonlat=OpenLayers.Layer.SphericalMercator.inverseMercator(point.x,point.y);point.x=lonlat.lon;point.y=lonlat.lat;return point;}
> OpenLayers.Layer.SphericalMercator.inverseMercator
function (x, y) {var lon=(x/20037508.34)*180;var lat=(y/20037508.34)*180;lat=180/Math.PI*(2*Math.atan(Math.exp(lat*Math.PI/180))-Math.PI/2);return new OpenLayers.LonLat(lon,lat);}

これは、http://en.wikipedia.org/wiki/Mercator_projection#Mathematics_of_the_projection　の素直な実装であるように見えます。
一方、順変換は、

> OpenLayers.Layer.SphericalMercator.projectForward
function (point) {var lonlat=OpenLayers.Layer.SphericalMercator.forwardMercator(point.x,point.y);point.x=lonlat.lon;point.y=lonlat.lat;return point;}
> OpenLayers.Layer.SphericalMercator.forwardMercator
function (lon, lat) {var x=lon*20037508.34/180;var y=Math.log(Math.tan((90+lat)*Math.PI/360))/(Math.PI/180);y=y*20037508.34/180;return new OpenLayers.LonLat(x,y);}

で、これも http://en.wikipedia.org/wiki/Mercator_projection#Mathematics_of_the_projection の素直な実装であるように見えます。

数値実験

> var m = OpenLayers.Layer.SphericalMercator
undefined
> var mp = m.forwardMercator(135, 35)
undefined
> mp.toString()
"lon=15028131.255,lat=4163881.1434847"
> m.inverseMercator(mp.lon, mp.lat).toString()
"lon=135,lat=35"

とりあえず、順変換して逆変換すれば元通りであることが確認できました。

Proj4js

Proj4js の場合も同様に、Safari の Web インスペクタを使って調べてみます。
情報源：http://trac.osgeo.org/proj4js/wiki/UserGuide#Basics
Web インスペクタを動作させるページ：http://proj4js.org/

> Proj4js.transform
function (source, dest, point) {
        if (!source.readyToUse || !dest.readyToUse) {
            this.reportError("Proj4js initialization for "+source.srsCode+" not yet complete");
            return point;
        }
        
        // Workaround for Spherical Mercator
        if ((source.srsProjNumber =="900913" && dest.datumCode != "WGS84") ||
            (dest.srsProjNumber == "900913" && source.datumCode != "WGS84")) {
            var wgs84 = Proj4js.WGS84;
            this.transform(source, wgs84, point);
            source = wgs84;
        }

        // Transform source points to long/lat, if they aren't already.
        if ( source.projName=="longlat") {
            point.x *= Proj4js.common.D2R;  // convert degrees to radians
            point.y *= Proj4js.common.D2R;
        } else {
            if (source.to_meter) {
                point.x *= source.to_meter;
                point.y *= source.to_meter;
            }
            source.inverse(point); // Convert Cartesian to longlat
        }

        // Adjust for the prime meridian if necessary
        if (source.from_greenwich) { 
            point.x += source.from_greenwich; 
        }

        // Convert datums if needed, and if possible.
        point = this.datum_transform( source.datum, dest.datum, point );

        // Adjust for the prime meridian if necessary
        if (dest.from_greenwich) {
            point.x -= dest.from_greenwich;
        }

        if( dest.projName=="longlat" ) {             
            // convert radians to decimal degrees
            point.x *= Proj4js.common.R2D;
            point.y *= Proj4js.common.R2D;
        } else  {               // else project
            dest.forward(point);
            if (dest.to_meter) {
                point.x /= dest.to_meter;
                point.y /= dest.to_meter;
            }
        }
        return point;
    }
> new Proj4js.Proj('EPSG:900913').inverse(new Proj4js.Point(15028131.255, 4163881.1434847)).toString()
"x=2.35619448986436,y=0.6108652381235775"
> new Proj4js.Proj('EPSG:900913').inverse(new Proj4js.Point(15028131.255, 4163881.1434847)).x
2.35619448986436
> new Proj4js.Proj('EPSG:900913').inverse(new Proj4js.Point(15028131.255, 4163881.1434847)).x * Proj4js.common.R2D
134.99999998120785

OpenLayers で順変換した座標を Proj4js で逆変換すると、（ほぼ）元通りになることが確認できました。肝心の逆変換式は、次のようになっていました。

> new Proj4js.Proj('EPSG:900913').inverse
function (p) {	

    var x = p.x - this.x0;
    var y = p.y - this.y0;
    var lon,lat;

    if (this.sphere) {
      lat = Proj4js.common.HALF_PI - 2.0 * Math.atan(Math.exp(-y / this.a * this.k0));
    } else {
      var ts = Math.exp(-y / (this.a * this.k0));
      lat = Proj4js.common.phi2z(this.e,ts);
      if(lat == -9999) {
        Proj4js.reportError("merc:inverse: lat = -9999");
        return null;
      }
    }
    lon = Proj4js.common.adjust_lon(this.long0+ x / (this.a * this.k0));

    p.x = lon;
    p.y = lat;
    return p;
  }