Calculate UTM grid convergence using UTM coordinates

Question

What would be the simplest way to arrive at a value for UTM grid convergence for a point given in UTM coordinates?

All advice I can find requires that the point be given in latitude/longitude. For example, there's a good post on this forum here:

Calculating grid convergence (True North to Grid North)

But that would first require conversion of point information from UTM to lat/lon. That's easy to do manually using any of the many online tools available. Coding an automated solution is also doable, though tedious. I managed to put together an Excel spreadsheet, benefiting from guidance offered here:

https://www.quora.com/How-do-I-convert-UTM-into-longitude-and-latitude-without-using-software

My particular application is in QGIS, where I'm trying to use the Expression String Builder to automatically set Print Layout > Main Properties > Map Rotation to orient the map so that up is true north. Ideally I'd like to be able to use x(@map_extent_center) and y(@map_extent_center) to get grid convergence at map center, then rotate the map by that amount. The ultimate goal there is to create maps similar to USGS 7.5-minute quadrangles, to be used for orienteering in the field.

Note that I did check out the QGIS plugin Lat Lon Tools. It includes functions for use in the Expression String Builder, but those functions convert lat/lon to UTM, not visa versa.

Answer 1

For an exhaustive mathematical discussion see Peter Osborne's "The Mercator Projections". Osborne's material was used to create the more succinct Wikipedia article here:

https://en.wikipedia.org/wiki/Transverse_Mercator_projection#Convergence

That discussion gives convergence angle (in radians!) as

arctan( tanh(x/k₀a) * tan(y/k₀a) )

where for UTM

k₀ = 0.9996

and

a = 6378137

But that doesn't go into details about the position of x relative to the central meridian. If we just plug in our UTM easting value, the result will be incorrect. For UTM, x starts counting 500000 meters to the west of the central meridian. So our equation would look more like

arctan( tanh((x-500000)/k₀a) * tan(y/k₀a) )

We have to know whether the UTM zone we're working with is in the northern or southern hemisphere because

y = northing in the north hemisphere

but

y = northing - 10,000,000 in the southern hemisphere

Note that if working in QGIS we don't have tanh(), but can substitute tanh()=sin()*cos().

arctan( sin((x-500000)/k₀a) * cos((x-500000)/k₀a) * tan(y/k₀a) )

Also, if wanting an answer in degrees instead of radians (i.e., so we can rotate our map), then we have to multiply by 180/pi()

Using the above, we can create an expression in QGIS to automatically rotate the map to true north. The expression might look something like this:

atan(sin((x(@map_extent_center)-500000)/6375585.745)*cos((x(@map_extent_center)-500000)/6375585.745)*tan(y(@map_extent_center)/6375585.745))*180/pi()

I have used that expression with success, but in the end it has seemed more convenient to use the idea proposed by Houska in the answer here:

Calculating convergence angle for point data in QGIS

where all the math is done under the hood by QGIS, so that we don't have to worry about what hemisphere we're in or even what CRS we're using. However, I could not get Houska's code to work as is, but had to adapt it as follows:

azimuth( make_point(x(@map_extent_center),y(@map_extent_center)), transform(translate(transform( make_point(x(@map_extent_center),y(@map_extent_center)),@map_crs,'EPSG:4326'),0,0.1),'EPSG:4326',@map_crs))/-pi()*180