Review the Mercator projection

[PaulWessel]

See the forum for background on this issue. That issues only affects data projections and not mapping. This feature request simply tries to clarify why there are differences and if we need to make any changes. It is not yet clear if this is a bug or a "feature".

Following complaints about the northing values in a Mercator projection when standard parallel other than Equator was used, we made this change to our codes in 6.1:

vmerc (setting up the projection) added a projected y-value for the standard latitude, previously always zero.

GMT->current.proj.j_yc = (fabs (slat) > 0.0) ? GMT->current.proj.j_x * d_log (GMT, tand (45.0 + 0.5 * slat)) : 0.0;

merc_sph (forward projection) subtracts this value from the result, yielding zero northing at the standard latitude:

*y = (fabs (lat) < 90.0) ? GMT->current.proj.j_x * d_log (GMT, tand (45.0 + 0.5 * lat)) - GMT->current.proj.j_yc : copysign (DBL_MAX, lat);

imerc_sph (inverse projection) then adds the northing back in before recovering the latitude:

*lat = atand (sinh ((y + GMT->current.proj.j_yc) * GMT->current.proj.j_ix));

What this does is to say that the projected y-value ("northing") shall be zero at the standard parallel instead of at the Equator.

This is a big part of the forum poster's differences. But it is not all. Reading Snyder on Mercator, he sets up everything for Equator. At the end of the chapter he says other standard parallels can be used, and then one makes the adjustment by multiplying by cos (stlat) and dividing for the inverse. He also says for ellipsoidal equations the scale is cos (stdlat)/sqrt (1 - e^2sin^2(stdlat)). Note for Equator these scales are 1.

In GMT 3(?) we had separate equations for spherical and ellipsoidal mercator, but @WalterHFSmith figured out that if we convert latitude to conformal latitude and use that in the spherical equations, we get the same result. So since at least GMT 4 we no longer have ellipsoidal versions and do it via the auxiliary latitudes required for the projection at hand. However, there are still some differences with Snyder:

1. We do, and have for a long time, scaled the projection by cos(stdlat), as suggested per Snyder.
2. We also divide by sqrt (1 - e^2sin^2(stdlat)), which should only happen for the ellipsoidal case.

Because maps are unaffected, and very few people use Mercator for data projection (most use UTM), it is possible that this is slightly wrong and nobody really noticed until the forum post. I see a few things that may be problematic or limiting:

• Given we use auxiliary latitudes (here, conformal) in the spherical equations, I am not sure about the 1/sqrt (1 - e^2sin^2(stdlat)), but a vague memory says that this is an overall scaling factor needed when using the conformal latitude. This I am sure @WalterHFSmith can set us straight on.
• Snyder did not discuss a northing adjustment (like ours) so that the projected y-value is zero at the standard latitude; he only added the scale.
• Our introduction of j_yc above means we cannot obtain Snyder's value for a (lon, lat) point when a standard parallel other than Equator is selected.
• If the poster tries to insert Equator as his standard parallel then he still cannot recreate the 6.0.0 and earlier results because now the cos(stdlat) is 1.

Perhaps @joa-quim can help sort this out. I know proj supports a standard parallel, but do they set y-out to zero there? If we want to make that an optional thing (the northing shift) then perhaps that should instead be dealt with via -C in mapproject and grdproject?

[joa-quim]

I know proj supports a standard parallel, but do they set y-out to zero there?

Apparently not because with proj we recover the 6.0 results (see forum post)

[PaulWessel]

If so, I suggest we revert that zero northing from the lib and add that capability to mapproject and grdproject instead. No good reason to be out of step with proj, no?

[joa-quim]

Yes, agree.

[PaulWessel]

#2287 was the reason for the change. Perhaps a special case when giving projected relative coordinates?

[PaulWessel]

I think I can fix this by computing that reference northing but not apply it in the projection, but use it when relative -R is given and add in the missing northing/easting there instead (when parsing -R).

[WalterHFSmith]

Hi. It is nice to hear that all that old work I did on auxiliary latitudes is still useful even in the Proj4 era. I am a little confused about what the issue is. I sort of half understand what Paul and Joaquim are talking about, but I am not sure I can make any logical sense of what the user wants. It is trivial to use a map projection and then subtract a constant from the Y coordinate, of course. If we are first converting geodetic lat to the conformal form of the auxiliary lat and then using any sort of mercator development of the sphere, then we will be conformal. But what exactly the user wants to have happen at a particular latitude parallel L then causes me to have some confusion. I think part of my confusion stems from the fact that the user seems to want (?) to measure y from some given (geodetic) latitude, L ? The user needs to understand that we can’t be both equidistant and conformal, and the Mercator is a conformal development, not an equidistant one. For a north-south true distance development of the sphere we would use the auxiliary meridional latitude, not the conformal one, and then some projection very different from Mercator. If the east-west extent of the mapped area is narrow then a transverse Mercator using the meridional latitude might work. If the east-west extent is larger a conformal conic would be better. If the user wants to make the horizontal (east-west) scale true on a parallel at geodetic latitude L, we can do that. We can also make the scale true on a certain parallel. Picture the sphere cut by a piece of paper rolled into a cylinder and cutting the sphere on the parallel where we want true scale. But note that if the user specifies true scale at latitude L, that is geodetic latitude L and the scale (equivalently, the width of the Mercator mapped area) must be calculated according to the circumference of that parallel of latitude on the ellipsoid. And because that is a calculation where we want true distance, we don’t transform L to the auxiliary conformal latitude to do that. (**Or we can, but see below.) The direct way to do it would be to compute W = 2 * pi * (east-west) * cos(lat_g) * a_e / sqrt(1 - esq * sin(L) * sin(L)) with esq = f*(2-f), f being the ellipsoid flattening and a_e the ellipsoid semi-major axis. Then W is the distance (in a_e units) across the map width from west to east (meridians bounding the mapped area). Then from W and the desired scale or width this defines the projection. (**) Alternatively we use the auxiliary conformal latitude, Lg, and we find an auxiliary radius, rPrime, such that the above W is also equal to W = 2 * pi (east-west) * rPrime. Then the scale can be computed using a spherical earth of a radius rPrime with latitude mapped to auxiliary conformal latitude. But whatever way it is done, the meaning of distance in Y is not going to have anything to do with north-south distance on the Earth, because the Mercator is not the right projection for that. I am sure I have just given a pedantic lecture that everyone knows and is totally useless for the question at hand. The problem is I don’t understand what the user wants, and whether it makes any sense to change GMT so it can do what s/he/they is asking it to do. walter Walter H F Smith Geophysicist, Laboratory for Satellite Altimetry NOAA NCWCP code E/RA-31 5830 University Research Court, room 3752 College Park MD 20740-3818 301-683-3377 (voice, my desk) 301-683-3300 (voice, reception) 301-683-3301 (fax, reception) Walter.HF.Smith@noaa.gov http://www.star.nesdis.noaa.gov/star/Smith_WHF.php

…

On Aug 28, 2020, at 5:44 PM, Paul Wessel ***@***.***> wrote: See the forum for background on this issue. That issues only affects data projections and not mapping. This feature request simply tries to clarify why there are differences and if we need to make any changes. It is not yet clear if this is a bug or a "feature". Following complaints about the northing values in a Mercator projection when standard parallel other than Equator was used, we made this change to our codes in 6.1: vmerc (setting up the projection) added a projected y-value for the standard latitude, previously always zero. GMT->current.proj.j_yc = (fabs (slat) > 0.0) ? GMT->current.proj.j_x * d_log (GMT, tand (45.0 + 0.5 * slat)) : 0.0; merc_sph (forward projection) subtracts this value from the result, yielding zero northing at the standard latitude: *y = (fabs (lat) < 90.0) ? GMT->current.proj.j_x * d_log (GMT, tand (45.0 + 0.5 * lat)) - GMT->current.proj.j_yc : copysign (DBL_MAX, lat); imerc_sph (inverse projection) then adds the northing back in before recovering the latitude: *lat = atand (sinh ((y + GMT->current.proj.j_yc) * GMT->current.proj.j_ix)); What this does is to say that the projected y-value ("northing") shall be zero at the standard parallel instead of at the Equator. This is a big part of the forum poster's differences. But it is not all. Reading Snyder on Mercator, he sets up everything for Equator. At the end of the chapter he says other standard parallels can be used, and then one makes the adjustment by multiplying by cos (stlat) and dividing for the inverse. He also says for ellipsoidal equations the scale is cos (stdlat)/sqrt (1 - e^2sin^2(stdlat)). Note for Equator these scales are 1. In GMT 3(?) we had separate equations for spherical and ellipsoidal mercator, but @WalterHFSmith figured out that if we convert latitude to conformal latitude and use that in the spherical equations, we get the same result. So since at least GMT 4 we no longer have ellipsoidal versions and do it via the auxiliary latitudes required for the projection at hand. However, there are still some differences with Snyder: • We do, and have for a long time, scaled the projection by cos(stdlat), as suggested per Snyder. • We also divide by sqrt (1 - e^2sin^2(stdlat)), which should only happen for the ellipsoidal case. Because maps are unaffected, and very few people use Mercator for data projection (most use UTM), it is possible that this is slightly wrong and nobody really noticed until the forum post. I see a few things that may be problematic or limiting: • Given we use auxiliary latitudes (here, conformal) in the spherical equations, I am not sure about the 1/sqrt (1 - e^2sin^2(stdlat)), but a vague memory says that this is an overall scaling factor needed when using the conformal latitude. This I am sure @WalterHFSmith can set us straight on. • Snyder did not discuss a northing adjustment (like ours) so that the projected y-value is zero at the standard latitude; he only added the scale. • Our introduction of j_yc above means we cannot obtain Snyder's value for a (lon, lat) point when a standard parallel other than Equator is selected. • If the poster tries to insert Equator as his standard parallel then he still cannot recreate the 6.0.0 and earlier results because now the cos(stdlat) is 1. Perhaps @joa-quim can help sort this out. I know proj supports a standard parallel, but do they set y-out to zero there? If we want to make that an optional thing (the northing shift) then perhaps that should instead be dealt with via -C in mapproject and grdproject? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[WalterHFSmith]

I guess I am thrown off by using the term “northing” to refer to a projected Y coordinate. “Northing” is a term of art in UTM but AFAIK it isn’t meaningful for a generic Mercator. And if it is to be used to indicate northward distance from some parallel, then the issue is how we compute the scale (equivalently, rPrime in my previous message), not how we use the auxiliary conformal latitudes. w Walter H F Smith Geophysicist, Laboratory for Satellite Altimetry NOAA NCWCP code E/RA-31 5830 University Research Court, room 3752 College Park MD 20740-3818 301-683-3377 (voice, my desk) 301-683-3300 (voice, reception) 301-683-3301 (fax, reception) Walter.HF.Smith@noaa.gov http://www.star.nesdis.noaa.gov/star/Smith_WHF.php

…

On Aug 28, 2020, at 5:44 PM, Paul Wessel ***@***.***> wrote: See the forum for background on this issue. That issues only affects data projections and not mapping. This feature request simply tries to clarify why there are differences and if we need to make any changes. It is not yet clear if this is a bug or a "feature". Following complaints about the northing values in a Mercator projection when standard parallel other than Equator was used, we made this change to our codes in 6.1: vmerc (setting up the projection) added a projected y-value for the standard latitude, previously always zero. GMT->current.proj.j_yc = (fabs (slat) > 0.0) ? GMT->current.proj.j_x * d_log (GMT, tand (45.0 + 0.5 * slat)) : 0.0; merc_sph (forward projection) subtracts this value from the result, yielding zero northing at the standard latitude: *y = (fabs (lat) < 90.0) ? GMT->current.proj.j_x * d_log (GMT, tand (45.0 + 0.5 * lat)) - GMT->current.proj.j_yc : copysign (DBL_MAX, lat); imerc_sph (inverse projection) then adds the northing back in before recovering the latitude: *lat = atand (sinh ((y + GMT->current.proj.j_yc) * GMT->current.proj.j_ix)); What this does is to say that the projected y-value ("northing") shall be zero at the standard parallel instead of at the Equator. This is a big part of the forum poster's differences. But it is not all. Reading Snyder on Mercator, he sets up everything for Equator. At the end of the chapter he says other standard parallels can be used, and then one makes the adjustment by multiplying by cos (stlat) and dividing for the inverse. He also says for ellipsoidal equations the scale is cos (stdlat)/sqrt (1 - e^2sin^2(stdlat)). Note for Equator these scales are 1. In GMT 3(?) we had separate equations for spherical and ellipsoidal mercator, but @WalterHFSmith figured out that if we convert latitude to conformal latitude and use that in the spherical equations, we get the same result. So since at least GMT 4 we no longer have ellipsoidal versions and do it via the auxiliary latitudes required for the projection at hand. However, there are still some differences with Snyder: • We do, and have for a long time, scaled the projection by cos(stdlat), as suggested per Snyder. • We also divide by sqrt (1 - e^2sin^2(stdlat)), which should only happen for the ellipsoidal case. Because maps are unaffected, and very few people use Mercator for data projection (most use UTM), it is possible that this is slightly wrong and nobody really noticed until the forum post. I see a few things that may be problematic or limiting: • Given we use auxiliary latitudes (here, conformal) in the spherical equations, I am not sure about the 1/sqrt (1 - e^2sin^2(stdlat)), but a vague memory says that this is an overall scaling factor needed when using the conformal latitude. This I am sure @WalterHFSmith can set us straight on. • Snyder did not discuss a northing adjustment (like ours) so that the projected y-value is zero at the standard latitude; he only added the scale. • Our introduction of j_yc above means we cannot obtain Snyder's value for a (lon, lat) point when a standard parallel other than Equator is selected. • If the poster tries to insert Equator as his standard parallel then he still cannot recreate the 6.0.0 and earlier results because now the cos(stdlat) is 1. Perhaps @joa-quim can help sort this out. I know proj supports a standard parallel, but do they set y-out to zero there? If we want to make that an optional thing (the northing shift) then perhaps that should instead be dealt with via -C in mapproject and grdproject? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[PaulWessel]

Yes, I think the user should use UTM for what he is doing. Since per Snyder the standard parallel does not mean resetting the y-values, we will revert to what we had and deal with the case that caused the initital change (#2287) in a different way.

[PaulWessel]

@WalterHFSmith commented:

We can also make the scale true on a certain parallel. Picture the sphere cut by a piece of paper rolled into a cylinder and cutting the sphere on the parallel where we want true scale. But note that if the user specifies true scale at latitude L, that is geodetic latitude L and the scale (equivalently, the width of the Mercator mapped area) must be calculated according to the circumference of that parallel of latitude on the ellipsoid. And because that is a calculation where we want true distance, we don’t transform L to the auxiliary conformal latitude to do that. (*Or we can, but see below.) The direct way to do it would be to compute
W = 2 * pi * (east-west) * cos(lat_g) * a_e / sqrt(1 - esq * sin(L) * sin(L))
with esq = f(2-f), f being the ellipsoid flattening and a_e the ellipsoid semi-major axis.
Then W is the distance (in a_e units) across the map width from west to east (meridians bounding the mapped area). Then from W and the desired scale or width this defines the projection.

This is exactly what GMT does. Good, so there is agreement between Walter's brain and GMT projection core. All is well.
I agree with the user being confused to think Mercator is a good projection to deal with distances. The only place in GMT where we consider relative Mercator distances is in cases like this:

gmt psbasemap -JM173/-42/16c -R-200/200/-200/200+uk -Baf -BWSne -P > t.ps

where the user wants a 200x200km map centered on the specified central longitude and standard latitude. This we now do under the hood without affecting the output of mapproject on actual coordinates.

[PaulWessel]

Actually, I must have gotten confused in 6.1. The other reason the numbers are different is that I used the conformal standard latitude in the equation Walter highlighted, but as he says, that should be geodetic. We then use the conformal in the spherical Mercator equation to get the y-value at that latitude. Fixing this in #4120.