Thursday, May 14, 2020

Oracle Functions Converting Between Latitude/Longitude and MGRS Coordinates

Let me provide the general context here. I worked on a gig in the recent past where, in addition to my DBA duties, I did some custom Oracle development work in a geospatial environment. For the sake of readability and context, suppose we were building a database of test aircraft bombing runs. We can use latitude (NS) and longitude (EW) to describe the target location. NATO militaries conventionally use an MGRS format 15-character length string (to specify within a meter or so precision):

The Military Grid Reference System (MGRS) is a grid-based system used to represent locations on the universal transverse Mercator (UTM) and universal polar stereographic (UPS) grid systems, expressed as an alphanumeric string. An MGRS coordinate defines an area on the Earth's surface as opposed to a specific point. A fully qualified MGRS string is 15 characters long and consists of the following three components: grid zone designation, 100,000-meter square identifier, and easting/northing.
Feeds might include data in latitude/longitude and/or MGRS format. Also, latitude/longitude data may be expressed in DMS (degrees/minutes/seconds) or decimal format (60 minutes per degree, 60 seconds per minute). I also wrote conversion functions for these, not included in this post. The functional developers I worked with preferred decimal format in working with related software like ArcGIS. For the purposes of the conversion functions below, the inputs or outputs are in decimal format.

One of the issues I had in dealing with a functional project manager was that, despite his decade of geospatial experience, he didn't understand the limitations of the data. For example I might have DMS data in whole numbers, and truncated/rounded fractions of seconds could effect computed coordinates.

I personally did not find Oracle-compatible functions for the conversions and so I adapted Java-like code from sources like Steve Brooker's embedded code and here. I cross-validated against online conversion websites like here.

Another thing I'm exemplifying here is my preferred style of technical documentation:

CONV_FROM_MGRS


Purpose:
To derive latitude/longitude in decimal format from MGRS in standard 15-character format. Among other reasons, in theory this provides an internal consistency check among source latitude, longitude, MGRS.

Syntax:
conv_from_mgrs(’15-character MGRS value’)

Output:
latitude decimal, longitude decimal
*** denotes an input error. There are explicit checks for input length and integer contents for the last 10 characters.

Code:

create or replace function conv_from_mgrs (mgrs in varchar2)
return varchar2
IS outlatlong varchar2(100);

pi number:=asin(1)*2;

type iarray is varray(20) of double precision;
i iarray;
type jarray is varray(20) of double precision;
j jarray;

lat double precision;
long double precision;

b1 varchar2(3);
b2 varchar2(2);
b3 double precision;
b4 double precision;

c double precision;
d varchar2(1);
e double precision;
f varchar2(1);
g integer;
h integer;
k double precision;
l double precision;
m double precision;
n double precision;
o double precision;
p double precision;
q double precision;
r double precision;
s double precision;
t double precision;
u double precision;
v double precision;
w double precision;
x double precision;
y double precision;
z double precision;
aa double precision;
ab double precision;
ac double precision;
ad double precision;
ae double precision;
af double precision;
ag double precision;
ah double precision;
ai double precision;
aj double precision;
ak double precision;
al double precision;

begin
if length(mgrs)=15 and invalid_integer_check(substr(mgrs,6,10))=0 then
b1:=substr(mgrs,1,3);
b2:=substr(mgrs,4,2);
b3:=to_number(substr(mgrs,6,5));
b4:=to_number(substr(mgrs,11,5));
c:=to_number(substr(b1,1,2));
d:=substr(b1,3,1);

e := ((c-1)*6-177)*pi / 180;

case
  when ((c-1) mod 3) = 0 then f:=instr('ABCDEFGH',substr(b2,1,1));
  when ((c-1) mod 3) = 1 then f:=instr('JKLMNPQR',substr(b2,1,1));
  else f:=instr('STUVWXYZ',substr(b2,1,1));
end case;

g:=instr('CDEFGHJKLMNPQRSTUVWXX',d);

case
  when ((c-1) mod 2) = 0 then h:=instr('ABCDEFGHJKLMNPQRSTUV',substr(b2,2,1))-1;
  else h:=instr('FGHJKLMNPQRSTUVABCDE',substr(b2,2,1))-1;
end case;


i:=iarray(1.1,2.0,2.8,3.7,4.6,5.5,6.4,7.3,8.2,9.1,0,0.8,1.7,2.6,3.5,4.4,5.3,6.2,7.0,7.9);
j:=jarray(0,2,2,2,4,4,6,6,8,8,0,0,0,2,2,4,4,6,6,6);
k:=i(g);
l:= j(g) + h / 10;
if (l < k)  then l:=l + 2; end if;
m:=f*100000.0 + b3;
n:=l*1000000 + b4;
m:= m-500000.0;
if (d < 'N') then n:=n-10000000.0; end if;
m:= m/0.9996;
n:=n/0.9996;
o:= n / 6367449.14570093;
p:= o + (0.0025188266133249035*sin(2.0*o)) + (0.0000037009491206268*sin(4.0*o)) + (0.0000000074477705265*sin(6.0*o)) + (0.0000000000170359940*sin(8.0*o));
q:=tan(p);
r:=q*q;
s:=r*r;
t:=cos(p);
u:=0.006739496819936062*power(t,2);
v:=40680631590769/(6356752.314*sqrt(1+u));
w:=v;
x:=1.0/(w*t);
w:=w*v;
y:=q/(2.0*w);
w:=w*v;
z:=1.0/(6.0*w*t);
w:=w*v;
aa:= q / (24.0*w);
w:=w*v;
ab:= 1.0 / (120.0*w*t);
w:=w*v;
ac:= q / (720.0*w);
w:=w*v;
ad:= 1.0 / (5040.0*w*t);
w:=w*v;
ae:= q / (40320.0*w);
af:= -1.0-u;
ag:= -1.0-2*r-u;
ah:= 5.0 + 3.0*r + 6.0*u-6.0*r*u-3.0*(u*u)-9.0*r*(u*u);
ai:= 5.0 + 28.0*r + 24.0*s + 6.0*u + 8.0*r*u;
aj:= -61.0-90.0*r-45.0*s-107.0*u + 162.0*r*u;
ak:= -61.0-662.0*r-1320.0*s-720.0*(s*r);
al:= 1385.0 + 3633.0*r + 4095.0*s + 1575*(s*r);
lat:= p + y*af*(m*m) + aa*ah*power(m,4) + ac*aj*power(m,6) + ae*al*power(m,8);
long:= e + x*m + z*ag*power(m,3) + ab*ai*power(m,5) + ad*ak*power(m,7);
lat:=lat*180/pi;
long:=long*180/pi;
outlatlong:=to_char(lat,'fm999.9999999')||','||to_char(long,'fm999.9999999');
return(outlatlong);
else return('***');
end if;
exception
when others then return('***');
end;
/




CONV_TO_MGRS

Purpose:
To derive latitude/longitude in decimal format from MGRS in standard 15-character format. Among other reasons, in theory this provides an internal consistency check among source latitude, longitude, MGRS.

Syntax:
conv_to_mgrs(’latitude decimal’, ‘longitude decimal’)
Note: single quotes are optional, but unquoted/undefined nonnumeric characters in inputs will result in an explicit SQL*PLUS error. Exceptions for quoted inputs are captured.

Output:
MGRS (standard 15-character string)
*** denotes an invalid input
(SQL*PLUS session serveroutput must be on for any relevant raised displayed messages, e.g., lat/long limits)

Code:
create or replace function conv_to_mgrs (lat in varchar2, lon in varchar2)
return varchar2
IS outmgrs varchar2(15);

latitude double precision;
longitude double precision;

invalid_latitude exception;
invalid_longitude exception;

pi double precision:=asin(1)*2;

outaa varchar2(5);
outab varchar2(5);

c double precision;
e double precision;
k double precision;
l double precision;
m double precision;
n double precision;
o double precision;
p double precision;
q double precision;
r double precision;
s double precision;
t double precision;
u double precision;
v double precision;
w double precision;
x double precision;
y double precision;
z double precision;
aa double precision;
ab double precision;
ad varchar2(1);
ae double precision;
af varchar2(1);
ag double precision;
ah varchar2(1);

begin

latitude:=to_number(lat);
longitude:=to_number(lon);

if abs(latitude) > 90 then raise invalid_latitude; end if;
if abs(longitude) > 180 then raise invalid_longitude; end if;

c:=1+floor((longitude+180)/6);
e:=c*6-183;
k:=latitude*pi/180;
l:=longitude*pi/180;
m:=e*pi/180;
n:=cos(k);
o:=0.006739496819936062*power(n,2);
p:=40680631590769/(6356752.314*sqrt(1+o));
q:=tan(k);
r:=q*q;
s:=(r*r*r)-power(q,6);
t:=l-m;
u:=1-r+o;
v:=5-r+9*o+4*(o*o);
w:=5-18*r + (r*r) + 14*o - 58*r*o;
x:=61-58*r + (r*r) + 270*o - 330*r*o;
y:=61-479*r + 179*(r*r) - (r*r*r);
z:=1385-3111*r + 543*(r*r)- (r*r*r);
aa:=p*n*t+(p/6*power(n,3)*u*power(t,3))+(p/120*power(n,5)*w*power(t,5))+(p/5040*power(n,7)*y*power(t,7));
ab:=6367449.14570093*(k - (0.00251882794504*sin(2*k)) + (0.00000264354112*sin(4*k)) - (0.00000000345262*sin(6*k)) + (0.000000000004892*sin(8*k))) + (q/2.0*p*power(n,2)*power(t,2)) + (q/24.0*p*power(n,4)*v*power(t,4)) + (q/720.0*p*power(n,6)*x*power(t,6)) + (q/40320.0*p*power(n,8)*z*power(t,8));
aa:=aa*0.9996 + 500000.0;
ab := ab*0.9996;
if ab < 0 then ab:=ab + 10000000.0; end if;
ad:=substr('CDEFGHJKLMNPQRSTUVWXZ',floor((latitude/8 + 11)),1);
ae:=floor(aa/100000);
case
  when ((c-1) mod 3) = 0 then af:= substr('ABCDEFGH',ae,1);
  when ((c-1) mod 3) = 1 then af:= substr('JKLMNPQR',ae,1);
  else af:=substr('STUVWXYZ',ae,1);
end case;
ag:=1+(floor(ab/100000)) mod 20;
case
  when ((c-1) mod 2) = 0 then ah:=substr('ABCDEFGHJKLMNPQRSTUV',ag,1);
  else ah:=substr('FGHJKLMNPQRSTUVABCDE',ag,1);
end case;
outaa:=ltrim(to_char((aa mod 100000),'00000'));
outab:=ltrim(to_char((ab mod 100000),'00000'));
outmgrs:=c||ad||af||ah||outaa||outab;
return(outmgrs);
exception
when invalid_latitude then dbms_output.put_line('latitude between -90 and 90'); return('***');
when invalid_longitude then dbms_output.put_line('longitude between -180 and 180'); return('***');
when others then return('***');
end;
/