Mercator projection routines for OpenStreetMaps.

So here’s a simple class which I cobbled together that converts to and from Mercator projected coordinates. The data stored in our database here at work is projected into Mercator coordinates, with the results stored apparently in meters.

I cobbled it together from the OpenStreetMap Wiki discussion on Mercator projections, and the complete class is below.

/**
 * Mercator projection routines, as used by OSM in their internal database.
 * Data stored in OSM is stored as points in X and Y, in meters projected
 * from a Mercator projected map. This allows conversion to and from those
 * points into lat/lng.
 * 
 * @author woody
 *
 */
public class Mercator
{
    final private static double R_MAJOR = 6378137.0;
    final private static double R_MINOR = 6378137.0;
//    final private static double R_MINOR = 6356752.3142;
 
    public static double[] toMerc(double x, double y) 
    {
        return new double[] {toMercX(x), toMercY(y)};
    }
 
    private static double toMercX(double lon) 
    {
        return R_MAJOR * Math.toRadians(lon);
    }
 
    private static double toMercY(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;
    }
    
    public static double[] fromMerc(double x, double y)
    {
    	return new double[] { fromMercX(x), fromMercY(y) };
    }

	private static double fromMercY(double y)
	{
		double temp = R_MINOR / R_MAJOR;
		double e = Math.sqrt(1.0 - (temp * temp));
		return Math.toDegrees(phi2(Math.exp(-y/R_MAJOR),e));
	}
	
	private static double phi2(double ts, double e)
	{
		int N_ITER=15;
		double HALFPI=Math.PI/2;
 
		double TOL=0.0000000001;
		double eccnth, phi, con, dphi;
		int i;
		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 && (0 != --i));
		return phi;
	}

	private static double fromMercX(double x)
	{
		return Math.toDegrees(x / R_MAJOR);
	}
}

Here’s the funny thing. While on the talk page they discuss using a Mercator projection using values for an elliptical Earth, it turns out that the projected data actually uses a spherical projection. Thus, the two values for R_MINOR, with the elliptical solution commented out.

Feh. Time to start a discussion on the OpenStreetMaps Wiki.

Leave a Reply

Please log in using one of these methods to post your comment:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s