Parameters of the coordinate system for calculating the rotation of labels along latitude lines

Question

I'm making a map of a rather large region and I need to set the angle of rotation for the labels so that they are parallel to the lines of latitude. I do it in three steps:

1. I take the central meridian of the current CRS of my project. For example, the WGS 84 / UTM zone 37N projection covers 36°E to 42°E, which means that its central meridian is in the middle of 36 and 42, that is, the central meridian is located at 39 degrees.

2. In the label rotation settings, I subtract from the central meridian (39 degrees) the X coordinate of the label (for example 47 degrees) and get a rotation angle of -8 degrees. (The label layer stores data in the EPSG:4326 WGS 84 coordinate system).

3. I multiply the result of the second step (-8 degrees) by some coefficient. For UTM zone 37N, the coefficient is approximately 0.72 (selected manually). For other projections, the coefficient may differ slightly.

As a result, everything is fine - the labels are parallel to the latitude lines in the middle and along the edges of the map.

But I would like to know two things:

1. Is it possible to automatically get the central meridian of the current projection in QGIS?

2. On what projection parameter does the multiplication coefficient in the third step depend? Is it possible to calculate it automatically based on some projection parameters?

Answer 1

The PROJ library provides those values, and you can access it in QGIS through Python functions importing the pyproj module.

I would think that versions of pyproj higher than 3.0.0 should support it. To find out the version installed with QGIS, open a Python console (Ctrl+Alt+P) and type:

import pyproj; pyproj.show_versions()

Following is a script that creates two new functions for the expression engine: longitude_origin( crs) and grid_convergence( lon, lat, crs):

---

from qgis.core import *
from qgis.gui import *
import pyproj

@qgsfunction(args='auto', group='Custom')
def longitude_origin(crs, feature, parent):
    """
    Returns the <b>Longitude of natural origin</b> (lon_0) parameter of the CRS.
    <h2>Example usage:</h2>
    <ul>
      <li>longitude_origin('EPSG:32637') -> 39</li>
      <li>longitude_origin(@layer_crs) -> Longitude of natural origin of 
      the current layer's CRS.</li>
    </ul>
    """
    crs_object = pyproj.CRS(crs)
    crs_coordinate_operation = crs_object.coordinate_operation
    if not crs_coordinate_operation:
        return f'No coordinate operation for {crs}'
    operation_parameters = crs_coordinate_operation.to_json_dict().get('parameters')
    lon_0_parameter = next(item for item in operation_parameters if
                            item['name'] == 'Longitude of natural origin')
    if not lon_0_parameter:
        return f'No lon_0 parameter for {crs}'
    lon_0_value = lon_0_parameter.get('value')
    return lon_0_value

@qgsfunction(args='auto', group='Custom')
def grid_convergence(lon, lat, crs, feature, parent):
    """
    Returns the <b>Grid convergence angle</b> (in degrees) for <i>Longitude</i> and
    <i>Latitude</i> coordinates in a <i>CRS</i>.<br>
    Grid convergence is the angle from True North to Grid North, positive clockwise.
    <h2>Example usage:</h2>
    <ul>
      <li>grid_convergence( 48, 46, 'EPSG:32637') -> 6.500059984899736</li>
    </ul>
    """
    crs_object = pyproj.CRS(crs)
    projection_object = pyproj.Proj(crs_object)
    convergence_angle = projection_object.get_factors(lon, lat, False, True).meridian_convergence
    return convergence_angle

---

Save it as a new file in the Function Editor tab of the expressions builder.

Then, you will find the new functions under the Custom group of functions.

1. If you type: longitude_origin( @layer_crs) the function outputs the longitude of natural origin for the current layer's CRS.

2. Not your multiplication parameter, but the angle from True North to Grid North is called Meridian convergence (or Grid convergence) angle and it depends of the coordinates of the position and the projection used. It is fully documented at The Mercator Projections. Peter Osborne, 2013..
   Now you can use the function: grid_convergence( 48, 46, 'EPSG:32637') to get the angle from True North to Grid North (positive clockwise), for the Latitude 46, Longitude 48 position in the EPSG:32637 projection. As allways, you can compute the coordinates using other functions. Since you want the angle for the label text, and it is positive counterclokwise, use the inverse of computed angle (prepending a negative sign to the function).

---

The generalized expression to rotate points labels perpendicular to the meridian convergence (in some cases this means aligning them to the parallel lines), regardless of the CRS of the layer to be labeled and the CRS of the map, would be:

- with_variable(
'convergence',
grid_convergence(
  x(
    transform(
      $geometry,
      @layer_crs,
      'EPSG:4326')),
  y(
    transform(
      $geometry,
      @layer_crs,
      'EPSG:4326')),
  @map_crs),
if(
  @convergence > 0,
  @convergence,
  @convergence + 360))

---

Links to documentation of pyproj's crs.CoordinateOperation class and proj.Factors class.