3

I have various point and polygon layers (shapefiles and geopackages) which I use regularly in QGIS (3.6.3), and often need to work their location as an 8 figure grid reference (British co-ordinates, EPSG-code 27700, same code used by each project). For polygons this means working out the centroid first and then converting this to a grid ref.

At the moment I am using the FSC plugin tools, 'OSGR Tool' and 'Add GRs to Layers Tool' to do this, but this means I have to go through a lengthy and somewhat risky process every time I have a new/edited feature which I need the info for.

I would prefer the grid ref to be automatically generated every time a point/polygon is added/edited, & then added to the 'Grid Ref' field for that feature, within the existing layer (see below). Is this possible?

Attributes Table

I presume it could be done using the 'default value' option for the existing field, but I cannot work out what expression I would need to do this.

10
  • Is there a separate grid layer, or is the centroid a part of the polygons attributes?
    – Erik
    Jan 8, 2020 at 13:08
  • I'm probably getting my terminology muddled here - by centroid I mean the central point inherent in any polygon - there is currently no separate grid layer
    – JimS-W
    Jan 8, 2020 at 14:09
  • And what do you mean by "grid reference"? The coordinates of said centroid?
    – Erik
    Jan 8, 2020 at 14:19
  • 2
    Well, problem: geom_to_wkt(centroid($geometry)) gives you the coordinates of the centroid - but only as numeric values. I have no idea how to convert these to the characteristic grid letters, sorry.
    – Erik
    Jan 8, 2020 at 16:08
  • 3
    This post has some code for converting grid reference to eastings and northings, as well as links to some other references. One of the links might help figure out how to do it. Or the other solution would be to download the British grid as a layer, and use an expression that references that layer.
    – csk
    Jan 8, 2020 at 16:43

1 Answer 1

4
+50

I wrote a custom function get_osgb36 which retrieves the geometry's centroid and computes the Ordnance Survey Grid code (or returns "Out of bounds").

It requires one argument, the precision; this defines the length of X/Y coordinates used for the number part of the code, e.g.:

get_osgb36(2) -> TQ3877
get_osgb36(0) -> TQ

The precision has a minimum of 1 and a maximum of 4, i.e. an 8-figure grid reference.

You can setup this function (follow the link for instructions on how to do this), using the code below.

Once set-up, the function will be available as an option within a new 'GIS SE' category in the Field Calculator’s Expressions lists. Choose the function and then enter the precision within the expression brackets (see above). Then set the expression as the default value for the relevant field, and select ‘Apply default value on update’ if required.

from qgis.core import *
from qgis.gui import *
from string import ascii_uppercase

GRID_LETTERS = [ascii_uppercase.replace("I", "")[i : i + 5] for i in range(0, 25, 5)]

@qgsfunction(args='auto', group='GIS SE', referenced_columns=[])
def get_osgb36(precision, feature, parent):
    """Get Ordnance Survey National Grid first letters (OSGB36)
    <h2>Argument:</h2>
    <ul>
      <li>Precision: if 1, you'll have the first 2 letters + 1 number for the X + 1 number for the y</li>
    </ul>
    <h2>Example usage:</h2>
    <ul>
      <li>get_osgb36(0) -> TQ</li>
      <li>get_osgb36(2) -> TQ3877</li>
    </ul>
    """
    oob_msg = "Out of bounds"

    if precision < 0 or precision > 4 or not isinstance(precision, int):
        precision = 3

    geom = feature.geometry()
    x = geom.centroid().asPoint().x()
    y = geom.centroid().asPoint().y()

    if y < -500000 or y >= 2000000:
        return oob_msg

    if x < -1000000 or x >= 1500000:
        return oob_msg

    # Ordnance Survey National Grid don't have I
    grid_letters = ascii_uppercase.replace("I", "")

    new_x = x + 1000000
    new_y = 2000000 - y
    x_letter_pos = int(new_x // 500000)
    y_letter_pos = int(new_y // 500000)
    y_letter_pos = y_letter_pos if y_letter_pos >= 0 else 0

    x_mod = int(new_x % 500000) // 100000
    y_mod = int(new_y % 500000) // 100000

    first_letter = GRID_LETTERS[y_letter_pos][x_letter_pos]
    second_letter = GRID_LETTERS[y_mod][x_mod]
    return f"{first_letter}{second_letter}{str(int(x))[1:1 + precision]}{str(int(y))[1:1 + precision]}"
1
  • That does the job perfectly
    – JimS-W
    Jan 31, 2023 at 9:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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