1

I have rewritten some Delphi function to ORACLE DB functions for converting from lat-long to utm and mgrs. Anyone care to validate the output with real data they know is correct?

package specification:

create or replace package gedaco as 
    function MGRS(lat in number, Lon in number,  a in number, InverseFlattening in number,  
                  Coding in number,  Digits in number) return varchar2;
    Function MGRSLatZone(lat in number) return varchar2;
    function SquareID(UTMzn in number, Northing in number, Easting in number, 
                      Coding in number) return varchar2;
    function UTM(lat in number, Lon in number, a in number, InverseFlattening in number)
                                return varchar2;
    Function UTMX(UTMs in varchar2) return number;
    Function UTMY(UTMs in varchar2) return number;
    function UTMZone(lat in number, Lon in number) return number;
end gedaco;

package body:

create or replace package body gedaco as 

function MGRS(lat in number, Lon in number,  a in number, InverseFlattening in number,
              Coding in number,  Digits in number) return varchar2 is 
  result varchar2(32);
  UTMs1 varchar2(32);
  E1 number;
  N1 number;
  Zn number;
  Lzn varchar2(32);
  Sq varchar2(32);
begin
  UTMs1 := UTM(lat, Lon, a, InverseFlattening) ;
  E1 := UTMX(UTMs1);
  N1 := UTMY(UTMs1);
  Zn := UTMZone(lat, Lon);
  Lzn := MGRSLatZone(lat);
  Sq := SquareID(Zn, N1, E1, Coding);
  result := replace(
      to_char(Zn,'00') || LZn || Sq || 
     to_char(round(E1 - 100000 * trunc(E1/100000)),'00000') || 
     to_char(round(N1 - 100000 * trunc(N1/100000)),'00000')

      ,' ', '');
  return result;
end MGRS;


Function MGRSLatZone(lat in number) return varchar2 is 
    result varchar2(1);
      GridZones CONSTANT varchar2(20) := 'CDEFGHJKLMNPQRSTUVW';
begin
  If (lat >= 72) Then Result := 'X';
  Else Result := substr(GridZones, Trunc((lat + 88) / 8), 1);
  End If;
  return result;
end MGRSLatZone;

function SquareID(UTMzn in number, Northing in number, Easting in number, 
                  Coding in number) return varchar2 is
    result varchar2(32);
    N number;
    E number;
    ZoneSet number;
    Col varchar2(32);
    Rov varchar2(32);
    Col1 CONSTANT varchar2(20) := 'ABCDEFGH';
    Col2 CONSTANT varchar2(20) := 'JKLMNPQR';
    Col3 CONSTANT varchar2(20) := 'STUVWXYZ' ;
    Row1 CONSTANT varchar2(20) := 'ABCDEFGHJKLMNPQRSTUV';
    Row2 CONSTANT varchar2(20) := 'FGHJKLMNPQRSTUVABCDE';
    Row3 CONSTANT varchar2(20) := 'LMNPQRSTUVABCDEFGHJK';
    Row4 CONSTANT varchar2(20) := 'RSTUVABCDEFGHJKLMNPQ';   
begin
  N := Trunc(Northing / 100000);
  N := N - 20 * Trunc(N / 20);
  E := Trunc(Easting / 100000);
  ZoneSet := UTMzn - 6 * Trunc(UTMzn / 6);
  If ((ZoneSet = 1) Or (ZoneSet = 4)) Then
     Col := SubStr(Col1, E, 1);
  End If;
  If ((ZoneSet = 2) Or (ZoneSet = 5)) Then
     Col := SubStr(Col2, E, 1);
  End If;
  If ((ZoneSet = 3) Or (ZoneSet = 0)) Then
     Col := SubStr(Col3, E, 1);
  End If;

  ZoneSet := ZoneSet - 2 * Trunc(ZoneSet / 2);
  If ((Coding = 1) And (ZoneSet = 1)) Then
    Rov := SubStr(Row1, N + 1, 1);
  End If;
  If ((Coding = 1) And (ZoneSet = 0)) Then
    Rov := SubStr(Row2, N + 1, 1);
  End If;
  If ((Coding = 2) And (ZoneSet = 1)) Then
    Rov := SubStr(Row3, N + 1, 1);
  End If;
  If ((Coding = 2) And (ZoneSet = 0)) Then
    Rov := SubStr(Row4, N + 1, 1);
  End If;

  Result:= Col || Rov;
  return result;
end SquareId;


function UTM(lat in number, Lon in number, a in number, InverseFlattening in number)
  return varchar2
  is result varchar2(320);
  ZoneWidth CONSTANT number := 6;
  CentralScaleFactor CONSTANT number := 0.9996;
  Zone1CentralMeridian CONSTANT number := -177;
  Zone0WestMeridian number;
  Zone0CentralMeridian number;
  FalseEasting CONSTANT number := 500000;

  Pi number;
  SemiMajorAxis number;
  Flattening number; Eccent2 number; Eccent4 number; Eccent6 number;
  A0 number; A2 number; A4 number; A6 number;
  LatRad  number;
  LonRad  number;
  Sin1Lat  number; Sin2Lat  number; Sin4Lat  number; Sin6Lat  number;
  Rho  number;
  Nu  number;
  Psi  number; Psi2  number; Psi3  number; Psi4   number;
  CosLat  number; CosLat2  number; CosLat3  number; CosLat4  number; CosLat5   number;
  CosLat6  number; CosLat7    number;
  TanLat  number; TanLat2  number; TanLat4  number; TanLat6    number;
  DifLon  number; DifLon2  number; DifLon3  number; DifLon4  number; DifLon5   number;
  DifLon6 number; DifLon7  number; DifLon8   number;
  DistOverMeridian    number;
  Zone  number;
  CentralMeridian  Integer;
  East1  number; East2  number; East3  number; East4  number;
  North1  number; North2  number; North3  number; North4  number;
  X  number;
  Y  number;
  Hemi  varchar2(1);
  FalseNorthing  number;

begin
Zone0WestMeridian := Zone1CentralMeridian - (1.5 * ZoneWidth);
Zone0CentralMeridian := Zone0WestMeridian + ZoneWidth / 2;
Pi := 3.141592653589793238462643383279502884197169399375105820974944592307816406;


SemiMajorAxis := 1000 * a  ;

Flattening := 1.0 / InverseFlattening   ;
Eccent2 := 2.0 * Flattening - (Flattening * Flattening);
Eccent4 := Eccent2 * Eccent2   ;
Eccent6 := Eccent2 * Eccent4 ;
A0 := 1 - (Eccent2 / 4.0) - ((3 * Eccent4) / 64.0) - ((5.0 * Eccent6) / 256.0);
A2 := (3.0 / 8.0) * (Eccent2 + (Eccent4 / 4.0) + ((15.0 * Eccent6) / 128.0)) ;
A4 := (15 / 256) * (Eccent4 + ((3.0 * Eccent6) / 4.0));
A6 := (35.0 * Eccent6) / 3072.0 ;
  --  ' Parameters to radians
    LatRad := lat / 180 * Pi;
    LonRad := Lon / 180 * Pi ;


  --  'Sin of latitude and its multiples
    Sin1Lat := sIn(LatRad) ;
    Sin2Lat := sIn(2 * LatRad) ;
    Sin4Lat := sIn(4 * LatRad);
    Sin6Lat := sIn(6 * LatRad);

  --  'Meridian Distance
    DistOverMeridian := SemiMajorAxis * 
                        (A0 * LatRad - A2 * Sin2Lat + A4 * Sin4Lat - A6 * Sin6Lat);


  --  'Radii of Curvature
    Rho := SemiMajorAxis * (1 - Eccent2) /Power( (1 - 
           (Eccent2 * Sin1Lat * Sin1Lat)) , 1.5);
    Nu := SemiMajorAxis /Power( (1 - (Eccent2 * Sin1Lat * Sin1Lat)) , 0.5);
    Psi := Nu / Rho  ;
    Psi2 := Psi * Psi ;
    Psi3 := Psi * Psi2;
    Psi4 := Psi * Psi3  ;

  --  'Powers of cos latitude
    CosLat := Cos(LatRad);
    CosLat2 := CosLat * CosLat  ;
    CosLat3 := CosLat * CosLat2 ;
    CosLat4 := CosLat * CosLat3 ;
    CosLat5 := CosLat * CosLat4 ;
    CosLat6 := CosLat * CosLat5 ;
    CosLat7 := CosLat * CosLat6 ;


--    'Powers of tan latitude
    TanLat := Tan(LatRad) ;
    TanLat2 := TanLat * TanLat ;
    TanLat4 := TanLat2 * TanLat2  ;
    TanLat6 := TanLat2 * TanLat4  ;

 --   'Zone
 --   'Zone := Int((Lon - Zone0WestMeridian) / ZoneWidth)
    Zone := UTMZone(lat, Lon)   ;



    CentralMeridian := Trunc((Zone * ZoneWidth) + Zone0CentralMeridian ) ;
    DifLon := (Lon - CentralMeridian) / 180 * Pi    ;
    DifLon2 := DifLon * DifLon  ;
    DifLon3 := DifLon * DifLon2 ;
    DifLon4 := DifLon * DifLon3 ;
    DifLon5 := DifLon * DifLon4 ;
    DifLon6 := DifLon * DifLon5 ;
    DifLon7 := DifLon * DifLon6 ;
    DifLon8 := DifLon * DifLon7 ;

    East1 := DifLon * CosLat  ;
    East2 := DifLon3 * CosLat3 * (Psi - TanLat2) / 6.0;
    East3 := DifLon5 * CosLat5 * (4.0 * Psi3 * (1.0 - 6.0 * TanLat2) + Psi2 * 
    (1.0 + 8.0 * TanLat2) -Psi * (2.0 * TanLat2) + TanLat4) / 120.0;
    East4 := DifLon7 * CosLat7 * (61.0 - 479.0 * TanLat2 + 179.0 * TanLat4 - TanLat6) 
    / 5040.0  ;
    X := CentralScaleFactor * Nu * (East1 + East2 + East3 + East4) + FalseEasting  ;

    If (lat >= 0) Then
      Hemi := 'N';
      FalseNorthing := 0;
    Else
      Hemi := 'S';
      FalseNorthing := 10000000;
    end if;

    North1 := Sin1Lat * DifLon2 * CosLat / 2.0 ;
    North2 := Sin1Lat * DifLon4 * CosLat3 * (4.0 * Psi2 + Psi - TanLat2) / 24.0 ;
    North3 := Sin1Lat * DifLon6 * CosLat5 * (8.0 * Psi4 * (11.0 - 24.0 * TanLat2)
              - 28.0 * Psi3 * (1.0 - 6.0 * TanLat2) +
     Psi2 * (1.0 - 32.0 * TanLat2) - Psi * (2.0 * TanLat2) + TanLat4) / 720;
    North4 := Sin1Lat * DifLon8 * CosLat7 * (1385 - 3111 * TanLat2 + 543 * 
              TanLat4 - TanLat6) / 40320.0 ;
    Y := CentralScaleFactor * (DistOverMeridian + Nu *
         (North1 + North2 + North3 + North4)) + FalseNorthing;

  Result := Zone || Hemi || ' ' || 
            to_char(round(X, 3),'0000000.000') || 
            to_char(round(Y, 3),'0000000.000');
  return result;
End UTM;


Function UTMX(UTMs in varchar2) return number
    is result number;
begin
  Result := to_number(substr(UTMs, 6, 11));
  return result;
End UTMX;


Function UTMY(UTMs in varchar2) return number
    is result number;
begin
  Result := to_number(substr(UTMs, 18, 11));
  return result;
End UTMY;


function UTMZone(lat in number, Lon in number) return number is
  result number;
  UTMZone number;
  e number;
  d number;
  ZoneWidth CONSTANT number := 6;
  Zone1CentralMeridian CONSTANT number := -177;
  Zone0WestMeridian number;
begin
  Zone0WestMeridian := Zone1CentralMeridian - (1.5 * ZoneWidth);
  d:=ZoneWidth;

   UTMZone := Trunc((lon - Zone0WestMeridian) / d);
    --Special Cases for Norway & Svalbard
    CASE 
    WHEN (lat > 55) AND (UTMZone = 31) AND (lat < 64) AND (lon >  2) THEN UTMZone := 32;
    WHEN (lat > 71) AND (UTMZone = 32) AND (lon <  9) THEN UTMZone := 31;
    WHEN (lat > 71) AND (UTMZone = 32) AND (lon >  8) THEN UTMZone := 33;
    WHEN (lat > 71) AND (UTMZone = 34) AND (lon < 21) THEN UTMZone := 33;
    WHEN (lat > 71) AND (UTMZone = 34) AND (lon > 20) THEN UTMZone := 35; 
    WHEN (lat > 71) AND (UTMZone = 36) AND (lon < 33) THEN UTMZone := 35;
    WHEN (lat > 71) AND (UTMZone = 36) AND (lon > 32) THEN UTMZone := 37;
    ELSE UTMZone := UTMZone;  
    END CASE;

  Result := UTMZone;
  return result;
end UTMZone;

end gedaco;

Function is used MGRS(:latitude, :longitude, 6378.137, 298.2572236, 1, 5) for WGS84 with 5 digits precision.
list of datums:

Datum           Radius          InverseFlattening
'WGS84',        6378.137,       298.2572236
'NAD27',        6378.2064,      294.9786982
'NAD83',        6378.137,       298.2572221
'WGS66',        6378.145,       298.25
'GRS67',        6378.16,        298.2472
'IAU68',        6378.16,        298.2472
'WGS72',        6378.135,       298.26
'Clarke66',     6378.2064,      294.9786982
'GRS80',        6378.137,       298.2572221
'Krasovsky',    6378.2064,      298.3
'Bessel',       6377.397155,    299.1528128

Just to be completely specific, I have of course done some random checks with Earth Point, but I am asking for someone to control my functions with a sustainable amount of data. A proper reply to my question can simply be test with 10 000 records without errors found, and if errors are found I am of course interested in the location not calculated correctly if that is possible.

1

I'm sorry for the late answer, but you might check out the NGA Gold Data. They have MGRS and UTM test points on various geographic coordinate reference systems (datums).

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.