[OSM-dev] Algorithm for converting lat/lon into Mercator as used in OSM PostGIS?

[Artem Pavlenko]

Nick,

see below

( ellipsoid params:
	p.a = 6378137.0
         p.b = 6356752.3142
    	p.k0=1.0
	t = 1.0 - b/a
         es = 2*t - t*t
	p.e = sqrt(es)
)

HTH
Artem


#ifndef MERCATOR_HPP
#define MERCATOR_HPP

#include <cmath>
#include "ellipsoid.hpp"
#include "projection.hpp"

#define EPS10 1e-10
#define HALFPI 1.5707963267948966
#define N_ITER 15

namespace proj {

     // determine small t
     inline double small_t (double phi, double sinphi, double e)
     {
         sinphi *=e;
         return ( tan(0.5 * (HALFPI - phi))/
                  pow ((1.0 - sinphi) / (1.0 + sinphi), 0.5 * e));
     }

     // determine latitude angle
     inline double phi2 ( double ts, double e)
     {
         double eccnth = 0.5 * e;
         double Phi = HALFPI - 2.0 * atan(ts);
         double dphi;
         int i = N_ITER;

         do {
             double con = e * sin(Phi);
             dphi = HALFPI - 2.0 *
                 atan(ts * pow((1.0 - con)/(1.0 + con), eccnth)) - Phi;
             Phi += dphi;

         } while ( fabs(dphi) > EPS10 && --i) ;
         return Phi;
     }

     bool merc_forward (double& x, double& y, params const& p)
     {
         if (fabs(fabs(y) - HALFPI) <= EPS10) return false;
         x = p.k0 * x;
         y = -p.k0 * log(small_t(y, sin(y), p.e));
         x *= p.a;
         y *= p.a;
         return true;
     }

     bool merc_inverse (double& x, double& y, params const& p)
     {
         x /= p.a;
         y /= p.a;
         x = x / p.k0;
         y =  phi2(exp(- y/p.k0),p.e);
         if ( y == std::numeric_limits<double>::infinity())
             return false;
         return true;
     }
}

#endif // MERCATOR_HPP

On 10 Mar 2007, at 22:25, Nick Whitelegg wrote:

> Hello everyone,
>
> Can anyone point me to any code which will do the conversion from  
> WGS84
> lat/lon into the Mercator system used in the OSM PostGIS database?
>
> I've tried using the Mercator conversion code present in various  
> places on
> SVN, mostly for utilities written pre  Mapnik, e.g:
>
> double mercator(double normalLat)
> {
>     double a = log (tan ( (M_PI/4.0) + (normalLat/180.0*M_PI/2.0))) *
>             180/M_PI;
>     return a  * 20037508.34 / 180;
> }
>
> but this gives me the wrong answer (out by about 30km at 51 N).
>
> I've also tried examining the PROJ.4 code used by Mapnik but it's  
> hideously
> convoluted. If there's a quicker route I'd be grateful of it :-)
>
> I'm aware Mapnik can do this for me, but (in Freemap), to avoid  
> server load
> and unnecessary projection conversions on the server I want to have  
> the bbox
> of the Freemap renderer in Mercator, and do all my conversions  
> client side.
>
> Thanks,
> Nick
>
> _______________________________________________
> dev mailing list
> dev at openstreetmap.org
> http://lists.openstreetmap.org/cgi-bin/mailman/listinfo/dev
>