Skip to content

Instantly share code, notes, and snippets.

Show Gist options
  • Select an option

  • Save Technoboggle/4fc1fdd97482f6214fd333f103793803 to your computer and use it in GitHub Desktop.

Select an option

Save Technoboggle/4fc1fdd97482f6214fd333f103793803 to your computer and use it in GitHub Desktop.
PL/PGSQL function to convert text Grid Reference into OSGB Geometry
-- convert a grid-reference E.g. SH123456 into a northing easting geometry for postgis with SRID=27700 (OSGB36)
-- standing on the shoulders of...
-- http://www.movable-type.co.uk/scripts/latlong-os-gridref.html
-- https://github.com/chrisveness/geodesy/blob/master/osgridref.js [MIT]
-- consider this MIT licensed also.
CREATE OR REPLACE FUNCTION get_geom_from_grid_ref(IN grid_ref character varying)
RETURNS public.geometry
LANGUAGE 'plpgsql'
AS $BODY$
DECLARE
parts text[];
l1 integer; l2 integer;
e100km integer; n100km integer;
easting varchar; northing varchar;
BEGIN
parts := regexp_matches(grid_ref, '([S,N,H,J,O,T])([A,B,C,D,E,F,G,H,J,K,L,M,N,O,P,Q,R,S,T,U,V,W,X,Y,Z]) ?([0-9]{3,5}) ?([0-9]{3,5})');
IF parts is null or array_length(parts, 1) = 0 THEN return null; END IF;
-- // get numeric values of letter references, mapping A->0, B->1, C->2, etc:
-- var l1 = gridref.toUpperCase().charCodeAt(0) - 'A'.charCodeAt(0);
-- var l2 = gridref.toUpperCase().charCodeAt(1) - 'A'.charCodeAt(0);
-- // shuffle down letters after 'I' since 'I' is not used in grid:
-- if (l1 > 7) l1--;
-- if (l2 > 7) l2--;
l1 := ascii(parts[1]) - ascii('A');
l2 := ascii(parts[2]) - ascii('A');
IF l1 > 7 THEN l1 := l1 - 1; END IF;
IF l2 > 7 THEN l2 := l2 - 1; END IF;
-- // convert grid letters into 100km-square indexes from false origin (grid square SV):
-- var e100km = ((l1-2)%5)*5 + (l2%5);
-- var n100km = (19-Math.floor(l1/5)*5) - Math.floor(l2/5);
e100km := ((l1-2)%5)*5 + (l2%5);
n100km := (19-floor(l1/5)*5) - floor(l2/5);
IF (e100km<0 or e100km>6 or n100km<0 or n100km>12) THEN
RAISE EXCEPTION 'Invalid grid reference: %', grid_ref;
END IF;
easting := e100km::varchar || rpad(parts[3], 5, '0');
northing := n100km::varchar || rpad(parts[4], 5, '0');
return public.ST_GeomFROMEWKT('SRID=27700;POINT(' || easting || ' ' || northing || ')');
END
$BODY$;
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment