5

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. Step 1. Find out the central meridian of the current projection.

  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). Step 2 and 3. Calculate the angle of rotation of the label.

  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.The result is that the labels are parallel to the latitude lines.

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?

0

1 Answer 1

5

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.

8
  • 1
    Thank you very much, the function grid_convergence( 48, 46, 'EPSG:32637') works perfectly instead of the expression I wrote. As parameters, I used the coordinates of the labels and the CRS variable of the current project: -grid_convergence( "700k_lbl_x", "700k_lbl_y", @map_crs ). The longitude_origin(@layer_crs) function also works great, just note that I used the @map_crs variable instead of @layer_crs. Jun 20, 2022 at 9:21
  • Does this function work for a large area or a custom CRS (Ex: EPSG:54032) ? I'm getting the > Eval Error: Invalid projection: EPSG:54032: (Internal Proj Error: proj_create: crs not found) < when i type the expression -grid_convergence( 20, 40, 'EPSG:54032') after following @Gabriel_De_Luca steps.
    – Iven Pepa
    Sep 3, 2022 at 20:02
  • @IvenPepa, EPSG:54032 doesn't exist. Do you mean ESRI:54032? Sep 3, 2022 at 21:08
  • Yes I forgot that was the ID["ESRI",54032] World_Azimuthal_Equidistant
    – Iven Pepa
    Sep 3, 2022 at 21:20
  • @IvenPepa, Use it as 'ESRI:54032', it must work. In some QGIS installations may be shown as EPSG authority in the database, but its authority was always ESRI for PROJ, therefore for pyproj, therefore for the script. Sep 3, 2022 at 21:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.