Skip to content

Instantly share code, notes, and snippets.

@robshep
Last active March 12, 2019 13:37
Show Gist options
  • Save robshep/bd61b868c9d68ce6baa81d64d8e4149c to your computer and use it in GitHub Desktop.
Save robshep/bd61b868c9d68ce6baa81d64d8e4149c 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]{5}) ?([0-9]{5})$');
if parts is null or array_length(parts, 1) = 0 then
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]{4}) ?([0-9]{4})$');
end if;
if parts is null or array_length(parts, 1) = 0 then
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}) ?([0-9]{3})$');
end if;
-- abandon.. not a grid ref
if parts is null or array_length(parts, 1) = 0 then
RAISE NOTICE 'Invalid grid reference: %', grid_ref;
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$;
-- convert SRID:27700-based OSGB36 coordinate point geometry to a OS Grid Reference E.g. SH123456
-- 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_grid_ref_from_geom(geom public.geometry)
RETURNS character varying
LANGUAGE plpgsql
AS $function$
DECLARE
e integer; n integer;
e100k integer; n100k integer;
l1 integer; l2 integer;
letterPair varchar;
begin
e := ST_X(geom)::integer;
n := ST_Y(geom)::integer;
-- // get the 100km-grid indices
e100k := floor(e/100000);
n100k := floor(n/100000);
-- // translate those into numeric equivalents of the grid letters
-- var l1 = (19-n100k) - (19-n100k)%5 + Math.floor((e100k+10)/5);
-- var l2 = (19-n100k)*5%25 + e100k%5;
l1 := (19-n100k) - (19-n100k)%5 + floor((e100k+10)/5);
l2 := (19-n100k)*5%25 + e100k%5;
-- // compensate for skipped 'I' and build grid letter-pairs
-- if (l1 > 7) l1++;
-- if (l2 > 7) l2++;
IF l1 > 7 THEN l1 := l1 + 1; END IF;
IF l2 > 7 THEN l2 := l2 + 1; END IF;
-- var letterPair = String.fromCharCode(l1+'A'.charCodeAt(0), l2+'A'.charCodeAt(0));
letterPair := chr( ascii('A') + l1 ) || chr( ascii('A') + l2 );
-- // strip 100km-grid indices from easting & northing, and reduce precision
-- e = Math.floor((e%100000)/Math.pow(10, 5-digits/2));
-- n = Math.floor((n%100000)/Math.pow(10, 5-digits/2));
e := floor(e%100000);
n := floor(n%100000);
-- // pad eastings & northings with leading zeros (just in case, allow up to 16-digit (mm) refs)
-- e = ('00000000'+e).slice(-digits/2);
-- n = ('00000000'+n).slice(-digits/2);
-- return letterPair + ' ' + e + ' ' + n;
return letterPair || e || n;
END
$function$
@atph
Copy link

atph commented Aug 22, 2018

Note to future me!
Please consider using the updated version of function_get_gridref_from_geom.sql here: https://gist.github.com/atph/0829d8d645720e679d8d04cbd7cfd5de#file-function_get_gridref_from_geom-sql

The code listed here does not provide the correct grid references under certain conditions:
346728, 800834 becomes NJ46728834 which should be NJ4672800834
346997, 804999 becomes NJ469974999 which should be NJ4699704999

@Radders500
Copy link

Hi, I was after a lat long to easting northing function in MS SQL. I've searched the web but couldn't find one can you help at all?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment