How to do the inverse Mercator projection

Question

Here is some code I started writing, the Mercator projection is based on this answer.

using System;
using System.Drawing;

namespace CoordinatesTool
{
    public class GeoPoint
    {
        public double Longitude { get; set; }
        public double Latitude { get; set; }

        public string ToString()
        {
            return Latitude + "," + Longitude;
        }

        public PointF ToMercator(int width, int height)
        {
            var x = (float)((Longitude + 180) * width / 360);

            var latRadians = Latitude * Math.PI / 180;
            var yTransformed = Math.Log(Math.Tan((latRadians / 2) + (Math.PI / 4)));
            var yScaled = (float)((height / 2.0) - (width * yTransformed / (2 * Math.PI)));

            return new PointF(x, yScaled);
        }

        public static GeoPoint FromMercator(PointF point, int width, int height)
        {
            return FromMercator(point.X, point.Y, width, height);
        }

        public static GeoPoint FromMercator(double x, double y, int width, int height)
        {
            // No clue what to do here
        }
    }
}

My goal is to use this utility class in a WinForms application. I'm using it with this map: http://en.wikipedia.org/wiki/File:Mercator-projection.jpg (width: 2048, height: 1588).

The Mercator inversion is working quite well (however, I suspect it's not very accurate in the arctic/antarctic reagions).

But the inverse Mercator projection really leaves me puzzled. I played around with the solution proposed in another question, but couldn't get anywhere. Especially I don't understand the Gudermannian (and inverse) function, the DEGREES_PER_RADIAN and RADIANS_PER_DEGREE constants, and how I should convert the y value into a latitude to call the GudermannianInv() function with.

EDIT: Here is how I tried how to do the inverse projection:

Starting with yScaled (parameter y in FromMercator function):

var yTransformed = 2 * Math.PI * (height / 2.0 - yScaled) / width;
var latRadians = Math.Arctan(Math.Pow(Math.E, yTransformed) - Math.PI / 4) / 2;
// ...

Answer 1

Here are some bits of what you seek:

1. radians * degrees/radians == degrees : degrees_per_radian is just a way of expressing degrees/radians in 'English' rather than in 'maths'. radians_per_degree is left as an exercise for the reader. So these two constants represent the numbers you use when converting between angles in degrees and angles in radians.

2. Looking at the code you've posted, you have the lines to convert latRadians to y. It looks straightforward to implement code to invert those operations. You need to 'un'-scale and 'un'-transform 'yScaled'. For example, treating the line (part):

   yScaled = ((height / 2.0) - (width * yTransformed / (2 * Math.PI)))

   as a mathematical equation, you need to solve for yTransformed in terms of yScaled.

3. As for the Gudermannian, the question you refer to implements that in one line of code and the inverse Gudermannian in 3 lines. The Gudermannian is just a way of transforming a circular measurement (such as a measurement in degrees or radians) into a linear measurement (such as one in centimetres on a chart to be published on paper). In particular, and pertinent here, the Gudermannian will transform a latitude into a linear distance from 0.

EDIT

OK, so looking a bit closer, in your original code you transform the latitude from an angular measurement into a linear one with the line:

yTransformed = Math.Log(Math.Tan((latRadians / 2) + (Math.PI / 4)))

I think that you should probably replace this with a call to the inverse Gudermannian, as indicated in the answer to the question you link to. I suspect that your home-brewed transformation is stumbling over extreme points in the tan/arctan functions. You would then use the direct Gudermannian in the un-transformation of course.