Math Help - Mercator projection and scale.

1. Mercator projection and scale.

I need to draw maps on a computer screen.

From Wikipedia, I found that for a given latitude, the "y" value for a Mercator projection is given by...

y = ln(tan(Latitude) + (1 / Cos(Latitude)))

... so I wrote a program to calculate "y" for latitudes from 0.0 to 89.9.

Now from a "sit back and look at it" position, the "y" value increases the further you move from the equator, so if I have a 1 degree square cell, when I plot it, with say 100 pixels for 1 degree of longitude, the scaling due to the projection should mean that I multiply that same 100 in longitude by the factor I get from the formula above, and that I should end up with a figure larger then 100, i.e. that "y" calculated above should be greater then 1 for all values of latitude.

It is not. I have to set latitude = 50.0 before the result of the above formula passes 1. I am multiplying the latitude in degrees by Pi/180.0 to convert the degrees to Radians. I am using 3.14159265358979323846 as Pi. All of my variables are either declared to be double precision floating points, or explicitly stated with a decoimal point to force float operation by the compiler.

Is the formula wrong? Am I doing something silly?

2. There's nothing wrong with the formula, but it uses radians throughout. So you should convert your "latitude = 50.0" to "latitude = 0.87266..." as the point where the formula passes 1.

3. I am multiplying the latitude in degrees by Pi/180.0 to convert the degrees to Radians. I am using 3.14159265358979323846 as Pi.

When my loop variable is 50.0 it is converted to Radians by multiplying by Pi/180.0 giving 0.87266462599716477 which goes into the formula and comes out as 1.0106831886830212.

4. What I am saying is that the formula for y gives the answer in radians. If you want it in degrees then you should use the formula $y = \tfrac{180}\pi\ln(\tan(\text{Latitude}) + (1 / \cos(\text{Latitude})))$.

5. That is even worse. When I do that, the value returned by the formula is less then 1 for all values.

The Wikipedia article says...

The following equations determine the x and y coordinates of a point on a Mercator map from its latitude φ and longitude λ (with λ0 being the longitude in the center of map):
This is the inverse of the Gudermannian function:

This is the Gudermannian function:

The scale is proportional to the secant of the latitude φ, getting arbitrarily large near the poles, where φ = ±90°. Moreover, as seen from the formulas, the pole's y is plus or minus infinity.
I read that as this being a formula which scales Latitude and Longitude to x and y coordinates on a screen, (or whatever).

What I see is that, assuming the centre of the map is zero, x = Longitude - zero, so I draw the point 00:00:00N 00:00:00E at the origin. The point 00:00:00N 01:00:00E at (1 - 0) * (how many pixels per degree). and so on. Latitude being Zero.

The point 01:00:00N 00:00:00E is still drawn at x = 0, (from above), but the y point is now...

rLat = 1.0 * (180.0 / Pi)
y = ln(tan(rLat) + (1 / cos(rLat))) * (how many pixels per degree)

... to generate the necessary exageration as I move away from the equator.

To acheive the exageration y must surely be >1 for all non-zero values, but it is not.

6. Okay, let's use that example of Latitude = 50°. If $\phi = 50^\circ$ then $\tan\phi\approx1.19$ and $\sec\phi\approx1.56$ (just working to 2 decimal places). Then my formula for y gives $y = \tfrac{180}\pi\ln(1.19 + 1.56)\approx57.91$. That is somewhat greater than 50, as one would wish.

In fact, the "exaggeration" in y at latitiude $\phi$ is equal to the derivative of $\ln(\tan\phi + \sec\phi)$, which is $\sec\phi$, and this is always greater than 1.

7. Arrgh!

When I do that with a calculator, it works as you say, however, the program I have here does not give the same answer. The 1.19 and 1.56 agree, I can see them in the debugger, but then the next step gives a different result.

I will look at it further.

8. Fixed!

Works now, you won't believe how stupid I was being, when you said convert the degrees to Radians, I knew I had done that, but also knew I had "the conversion factor", already calculated as I didn't want the overhead of calculating it every time.

Trouble is, the factor I had was Pi / 180.0 which is fine for where I was using it, but not where I added it later on your reco!

Suffice to say, I now have two factors r2d and d2r.

I hope I don't need another r2d factor otherwise I will have r2d2 and risk Imperial intervention.