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.
Best Answer
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).