สูตรที่ใช้ในการคำนวณ ผมอ้างอิงมาจาก http://www.uwgb.edu/dutchs/UsefulData/UTMFormulas.htm เขียนโดย Steve Dutch กล่าวถึงสูตรที่นำมาใช้คำนวณแปลงค่าพิกัด โดยอ้างถึง U.S. Army Corps of Engineer และ USGS (U.S. Geological Survey Professional) ถ้าเป็นการแปลงจาก Lat/Long ไป UTM เอกสารของ U.S. Army Corps of Engineer ทำได้ดีกว่าชัดเจนกว่า แต่ถ้าเป็นการแปลงพิกัดจาก UTM ไป Lat/Long USGS ทำได้ดีกว่า เราจะใช้สูตรที่ Steve Dutch ได้อ้างถึงและแสดงไว้ในเอกสาร นำมาโค๊ดเป็น Lazarus (Delphi) และที่น่าสนใจคือมีไฟล์ของ Excel สำหรับการคำนวณไว้ให้ดาวน์โหลดได้ด้วย ทำไว้ได้ดีมากครับ อยู่ตรงนี้ http://www.uwgb.edu/dutchs/UsefulData/UTMConversions1.xls UTM Zone UTM Zone No GeoEllipsoids Unit เราจะเริ่มด้วย Unit “GeoEllipsoids” (เก็บอยู่ในไฟล์ geoellipsoids.pas ที่ Lazarus พยายามให้เราตั้งชื่อเป็น lower case) เป็นยูนิตที่ declare class รูปแบบของ Ellipsoid ชื่อ TEllipsoid และมีสร้าง collection สำหรับเก็บ Ellipsoid ไว้ด้วยในยูนิตเดียวกัน unit GeoEllipsoids;
{$mode objfpc}{$H+}
interface
uses
Classes, Contnrs, SysUtils;
type
TEllipsoid = class
private
FEllipName : string;
FMajorAxis : double;
FInvFlat : double;
public
property EllipsoidName : string read FEllipName write FEllipName;
property MajorAxis : double read FMajorAxis write FMajorAxis;
property InvFlattening : double read FInvFlat write FInvFlat;
constructor Create ;
destructor Destroy ; override ;
end;
TEllipsoidList = class (TStringList)
private
protected
function GetEllipsoid(idx : Integer ) : TEllipsoid ;
public
constructor Create ;
destructor Destroy ; override ;
procedure AddEx (const Name : string; const MajorAxisA : double; const InvFlat : double);
function FindEx(const Name : string): TEllipsoid ;
end;
implementation
//TEllipsoid
constructor TEllipsoid.Create;
begin
inherited create;
end;
destructor TEllipsoid.Destroy;
begin
inherited Destroy;
end;
//TEllipsoidList - Collection of TEllipsoid.
constructor TEllipsoidList.Create ;
begin
inherited Create;
Sorted := True ;
AddEx('WGS84', 6378137, 298.257223563);
AddEx('GRS80', 6378137, 298.257222101);
AddEx('WGS72', 6378135, 298.260);
AddEx('Australian 1965', 6378160, 298.250);
AddEx('Krasovsky 1940', 6378245, 298.3);
AddEx('International (1924)', 6378388, 297);
AddEx('Clake 1880', 6378249.1, 293.460);
AddEx('Clarke 1866', 6378206.4, 294.980);
AddEx('Airy 1830', 6377563.4, 299.320);
AddEx('Bessel 1841', 6377397.2, 299.150);
AddEx('Everest 1830', 6377276.345, 300.8017);
AddEx('Hayford 1909', 6378388.0, 296.999362081575);
AddEx('North American 1927', 6378206.4, 294.978698213898);
AddEx('NAD 83', 6378137, 298.257223563);
end ;
destructor TEllipsoidList.Destroy ;
var
i : Integer ;
begin
for i:= Count -1 downto 0 do
begin
if Assigned( Objects[i] ) then
TEllipsoidList( Objects[i] ).Free ;
end ;
inherited Destroy;
end ;
procedure TEllipsoidList.AddEx(const Name : string ;
const MajorAxisA : double ;
const InvFlat : double);
var
e : TEllipsoid ;
pos : Integer ;
begin
if Find(Name, pos ) then begin
if Assigned(Objects[pos]) then
TEllipsoid(Objects[pos]).Free ;
end ;
pos := Add(Name) ;
e := TEllipsoid.Create ;
with e do
begin
EllipsoidName := Name;
MajorAxis := MajorAxisA;
InvFlattening := InvFlat;
end ;
Objects[pos] := e;
end;
function TEllipsoidList.FindEx(const Name : String): TEllipsoid ;
var
pos : Integer ;
begin
if Find( Name, pos ) then
Result := TEllipsoid(Objects[pos])
else
Result := nil ;
end ;
function TEllipsoidList.GetEllipsoid(idx : Integer) : TEllipsoid ;
begin
Result := TEllipsoid( Objects[ idx ] ) ; ;
end ;
end.จากโค๊ตด้านบนจะเห็นผม Declare class ชื่อ TEllipsoid TEllipsoidList ผมเพิ่มเมธอด FindEx สำหรับค้นหาและดึง object ของทรงรีออกมา และยังเพิ่ม AddEx สำหรับเผื่อไว้ในกรณีต้องการเพิ่มทรงรีอื่นในภายหลังได้ GeoCompute Unit ต่อไปเป็น Unit ที่ 2 ชื่อ “GeoCompute ” (เก็บอยู่ในไฟล์ geocompute.pas) unit GeoCompute;
{$mode objfpc}{$H+}
interface
uses
Classes, SysUtils, GeoEllipsoids, Math;
type
TZoneHemisphere = (hemEast, hemWest, hemNorth, hemSouth);
TUTMProjection = record
K0 : double;
FE : double;
FN : double;
ZoneNo : integer;
CM : double;
LatHem : TZoneHemisphere;
Longhem : TZoneHemisphere;
end;
TCoordinate = record
X : double;
Y : double;
end;
TGeo2UTM = class
private
FEllipsoid : TEllipsoid;
FGeoCoor : TCoordinate;
FUTMCoor : TCoordinate;
FUTMProj : TUTMProjection;
procedure SetGeoCoordinate (const geocoor : TCoordinate);
procedure SetEllipsoid(const ellipse : TEllipsoid);
public
property Ellipsoid : TEllipsoid write SetEllipsoid;
property GeoCoordinate : TCoordinate write SetGeoCoordinate; //input
property UTMCoordinate : TCoordinate read FUTMCoor; //output
property UTMProjection : TUTMProjection read FUTMProj; //output
procedure Compute;
constructor Create;
destructor Destroy; override;
end;
//On UTM Grid coordianates we need to know the zone no, latitude hemisphere and
//longitude hemisphere also. Then we can convert utm coordinates to geographic
//coordinate (latitude/longitude).
TUTM2Geo = class
private
FEllipsoid : TEllipsoid;
FGeoCoor : TCoordinate;
FUTMCoor : TCoordinate;
FUTMProj : TUTMProjection;
procedure SetUTMProjection (const utmproj : TUTMProjection);
procedure SetUTMCoordinate (const utmcoor : TCoordinate);
procedure SetEllipsoid(const Ellipsoid : TEllipsoid);
public
property UTMProjection : TUTMProjection read FUTMProj write SetUTMProjection;
property Ellipsoid : TEllipsoid write SetEllipsoid;
property UTMCoordinate : TCoordinate write SetUTMCoordinate; //input
property GeoCoordinate : TCoordinate read FGeoCoor; //output
procedure Compute;
constructor Create;
destructor Destroy; override;
end;
implementation
//TGeo2UTM
procedure TGeo2UTM.SetEllipsoid(const ellipse : TEllipsoid);
begin
FEllipsoid.EllipsoidName := ellipse.EllipsoidName;
FEllipsoid.MajorAxis := ellipse.MajorAxis;
FEllipsoid.InvFlattening := ellipse.InvFlattening;
end;
procedure TGeo2UTM.SetGeoCoordinate(const geocoor : TCoordinate);
var
dLat, dLong : double;
begin
FGeoCoor := geocoor;
dLat := FGeoCoor.Y;
dLong := FGeoCoor.X;
//set hemisphere
if (dLat >=0) then
FUTMProj.LatHem := hemNorth
else
FUTMProj.LatHem := hemSouth;
if (dLong >=0) then
FUTMProj.LongHem := hemEast
else
FUTMProj.LongHem := hemWest;
//compute zone
if (FUTMProj.LongHem = hemWest) then
FUTMProj.ZoneNo := trunc((180 + dLong) / 6) + 1
else if (FUTMProj.LongHem = hemEast) then
FUTMProj.ZoneNo := trunc(dLong/6) + 31;
//compute zone
if (FUTMProj.LatHem = hemNorth) then
FUTMProj.FN := 0
else if (FUTMProj.LatHem = hemSouth) then
FUTMProj.FN := 10000000;
FUTMProj.CM := 6 * FUTMProj.ZoneNo - 183;
end;
procedure TGeo2UTM.Compute;
var
K0, FE, FN : extended;
lat, long : extended; //in radian.
a, b, f, e, et2, rm, n, rho, nu, p, cm : extended;
M, A0, B0, C0, D0, E0 : extended;
Ki, Kii, Kiii, Kiv, Kv, A6 : extended;
begin
//Assign to new variable for easy looking at mathematic formula.
K0 := FUTMProj.K0;
FE := FUTMProj.FE;
FN := FUTMProj.FN;
a := FEllipsoid.MajorAxis;
f := 1 / FEllipsoid.InvFlattening;
long := PI/180 * FGeoCoor.X; //convert latitude to radian.
lat := PI/180 * FGeoCoor.Y;
cm := PI/180 * FUTMProj.CM;
p := long - cm; //already in radian.
b := a * (1 - f); //Minor axis;
// =.08 approximately. This is the eccentricity of the earth's elliptical crosssection.
e := sqrt(1 - b*b / (a*a));
//= .007 approximately. The quantity e' only occurs in even powers so it need only be calculated as e'2
et2 := power((e * a/b), 2);
n := (a - b) / (a + b);
//This is the radius of curvature of the earth in the meridian plane.
rho := a * (1 - e*e) / power (1 - power(e*sin(lat),2), 3/2);
//This is the radius of curvature of the earth perpendicular to the meridian plane.
nu := a / sqrt(1 - sqr(e*sin(lat)));
//=================== compute Meridian Arc===================================
//A' = a(1 - n + (5/4)(n2- n3) + (81/64)(n4- n5))
A0 := a * (1 - n + (5.0/4.0)*(n*n - n*n*n) + (81.0/64.0)*(power(n,4) - power(n,5)));
//B' = (3an/2)[1 - n + (7/8)(n2- n3) + (55/64)(n4- n5) ...]
B0 := (3*a*n/2) * (1 - n + (7/8)*(n*n - n*n*n) + (55/64)* power(n,4) - power(n,5));
//C' = (15 tan2/16)[1 - n + (3/4)(n2- n3) ...]
C0 := (15*a*n*n/16) * (1 - n + (3/4)*(n*n - n*n*n));
//D' = (35 tan3/48)[1 - n + (11/16)(n2- n3) ...]
D0 := (35*a*n*n*n/48) * (1 - n + (11/16)*(n*n - n*n*n));
//E' = (315 tan4/512)[1 - n ...]
E0 := (315*a*power(n,4)/512) * (1-n);
//M = A'lat - B'sin(2lat) + C'sin(4lat) - D'sin(6lat) + E'sin(8lat).
M := A0*lat - B0*sin(2*lat) + C0*sin(4*lat) - D0*sin(6*lat) + E0*sin(8*lat);
//===========================================================================
//===================Converting Lat/Long to UTM Grid Coordinate==============
Ki := M * K0;
Kii := K0 * nu * sin(2*lat) / 4;
Kiii := (K0 * nu * sin(lat) * power(cos(lat),3) / 24) * ((5 - tan(lat)*tan(lat)
+ 9*et2 * cos(lat)*cos(lat) + 4 * et2*et2 * power(cos(lat),4)));
Kiv := K0 * nu * cos(lat);
Kv := (K0 * nu * power(cos(lat),3)/6) * (1 - tan(lat)*tan(lat) + et2*cos(lat)*cos(lat));
//Now we 'v got Grid UTM Coordinates.
FUTMCoor.Y := FN + Ki + Kii*p*p + Kiii*power(p,4);
FUTMCoor.X := FE + Kiv*p + Kv*p*p*p;
end;
constructor TGeo2UTM.Create;
begin
inherited create;
FUTMProj.K0 := 0.9996;
FUTMProj.FE := 500000;
FEllipsoid := TEllipsoid.Create;
end;
destructor TGeo2UTM.Destroy;
begin
FEllipsoid.Free;
inherited Destroy;
end;
//TUTM2Geo
procedure TUTM2Geo.SetEllipsoid(const Ellipsoid : TEllipsoid);
begin
FEllipsoid.EllipsoidName := Ellipsoid.EllipsoidName;
FEllipsoid.MajorAxis := Ellipsoid.MajorAxis;
FEllipsoid.InvFlattening := Ellipsoid.InvFlattening;
end;
procedure TUTM2Geo.SetUTMProjection(const utmproj : TUTMProjection);
begin
FUTMProj.LongHem := utmproj.LongHem;
FUTMProj.LatHem := utmproj.LatHem;
if (utmproj.LatHem = hemNorth) then
FUTMProj.FN := 0
else if (utmproj.LatHem = hemSouth) then
FUTMProj.FN := 10000000;
FUTMProj.CM := 6 * utmproj.ZoneNo - 183;
if (FUTMProj.CM < 0) then
FUTMProj.Longhem := hemWest
else
FUTMProj.Longhem := hemEast;
end;
procedure TUTM2Geo.SetUTMCoordinate(const utmcoor : TCoordinate);
var
dN, dE : double;
begin
FUTMCoor := utmcoor;
dN := FGeoCoor.Y;
dE := FGeoCoor.X;
end;
procedure TUTM2Geo.Compute;
var
K0, FE, FN, a, b, f, e, et2, cm : double;
mu, M, e1 : double;
J1, J2, J3, J4 : double;
C1, T1, R1, N1, D : double;
fp, Q1, Q2, Q3, Q4, Q5, Q6, Q7 : double;
lat, long : double;
x, y : double;
begin
//Assign to new variable for easy looking at mathematic formula.
K0 := FUTMProj.K0;
FE := FUTMProj.FE;
FN := FUTMProj.FN;
a := FEllipsoid.MajorAxis;
f := 1 / FEllipsoid.InvFlattening;
x := (FE - FUTMCoor.X);
if (FUTMProj.LatHem = hemNorth) then
y := FUTMCoor.Y
else if (FUTMProj.LatHem = hemSouth) then
y := (10000000 - FUTMCoor.Y); //temporary
cm := PI/180 * FUTMProj.CM;
b := a * (1 - f); //Minor axis;
// =.08 approximately. This is the eccentricity of the earth's elliptical crosssection.
e := sqrt(1 - b*b / (a*a));
//= .007 approximately. The quantity e' only occurs in even powers so it need only be calculated as e'2
et2 := power((e * a/b), 2);
M := y / K0; // M is Meridian Arc;
//Compute footprint Latitude.
//mu = M/[a(1 - e2/4 - 3e4/64 - 5e6/256...)
mu := M/(a*(1 - e*e/4 - 3/64*power(e,4) - 5/256*power(e,6)));
//e1 = [1 - (1 - e2)1/2]/[1 + (1 - e2)1/2]
e1 := (1 - sqrt(1 - e*e))/(1 + sqrt(1 - e*e));
//J1 = (3e1/2 - 27e13/32 ..)
J1 := (3*e1/2 - 27*e1*e1*e1/32);
//J2 = (21e12/16 - 55e14/32 ..)
J2 := 21*e1*e1/16 - 55*power(e1,4)/32;
//J3=(151e13/96 ..)
J3 := 151*e1*e1*e1/96;
//J4 = (1097e14/512 ..)
J4 := 1097*power(e1,4)/512;
//Footprint Latitude fp = mu + J1sin(2mu) + J2sin(4mu) + J3sin(6mu) + J4sin(8mu)
fp := mu + J1*sin(2*mu) + J2*sin(4*mu) + J3*sin(6*mu) + J4 * sin(8*mu);
//C1=e'2cos2(fp)
C1 := et2*cos(fp)*cos(fp);
//T1 = tan2(fp)
T1 := tan(fp)*tan(fp);
//R1 = a(1-e2)/(1-e2sin2(fp))3/2
R1 := a*(1-e*e)/power(1-e*e*sin(fp)*sin(fp), 1.5);
//N1 = a/(1-e2sin2(fp))1/2
N1 := a / sqrt(1 -e*e*sin(fp)*sin(fp));
//D = x/(N1k0)
D := x / (N1 * K0);
//Compute Latitude
//Q1 = N1 tan(fp)/R1
Q1 := N1 * tan(fp)/R1;
//Q2 = (D2/2)
Q2 := D*D/2;
//Q3 = (5 + 3T1 + 10C1 - 4C12-9e'2)D4/24
Q3 := (5 + 3*T1 + 10*C1 - 4*C1*C1-9*et2)*power(D,4)/24;
//Q4 = (61 + 90T1 + 298C1 +45T12- 3C12-252e'2)D6/720
Q4 := (61 + 90*T1 + 298*C1 + 45*T1*T1 - 3*C1*C1 - 252*et2)* power(D,6)/720;
Q5 := D;
//Q6 = (1 + 2T1 + C1)D3/6
Q6 := (1 + 2*T1 + C1) *D*D*D/6;
//Q7 = (5 - 2C1 + 28T1 - 3C12+ 8e'2+ 24T12)D5/120
Q7 := (5 - 2*C1 + 28*T1 - 3*C1*C1 + 8*et2 + 24*T1*T1)*power(D,5)/120;
//lat
lat := 180/PI *(fp - Q1*(Q2 + Q3 + Q4));
//long = long0 + (Q5 - Q6 + Q7)/cos(fp)
long := 180/PI *(CM - (Q5 - Q6 + Q7)/cos(fp));
FGeoCoor.Y := lat;
FGeoCoor.X := long;
end;
constructor TUTM2Geo.Create;
begin
inherited create;
FEllipsoid := TEllipsoid.Create;
FUTMProj.K0 := 0.9996;
FUTMProj.FE := 500000;
end;
destructor TUTM2Geo.Destroy;
begin
FEllipsoid.Free;
inherited Destroy;
end;
end.
.จากโค๊ดด้านบนผม declare type ที่เป็น record ไว้ 2 แบบ คือ
TUTMProjection TCoordinate
ถัดลงมาอีกผม declare class ไว้ 2 class คือ TUTM2Geo และ TGeo2UTM
TGeo2UTM คำนวณค่าพิกัดจาก Geographic ไปเป็น UTM
property สำคัญคือ Ellipsoid เราต้อง set โดยการส่ง object ของ class TEllipsoid ที่ผู้ใช้เลือกมาให้
property อีกตัวคือ GeoCoordinate เป็นค่าพิกัด Lat, Long ที่เราส่งไปให้เพื่อจะแปลงค่าพิกัดเป็นพิกัดยูทีเอ็ม
Compute เป็น method สำคัญที่ทำหน้าที่คำนวณแปลงค่าพิกัด ในเมธอดนี้จะมีสูตรที่ผมอ้างอิงถึงดังที่กล่าวไปแล้วส่วนผลลัพธ์ส่งออกมาทาง property UTMCoordinate
และ property UTMProjection ก็เป็นผลลัพธ์อีกตัวที่ได้จากการคำนวณ ข้อดีของค่าพิกัดในรูป Geographic (Latitude/Longitude) สามารถคำนวณหาได้ว่าอยู่ Zone ไหนของยูทีเอ็ม, หา Central Meridian ได้
TUTM2Geo คำนวณค่าพิกัดจาก UTM ไปเป็น Geographic
Ellipsoid property เช่นเดียวกันกับ class TGeo2UTMUTMProjection property ค่าพารามิเตอร์ของยูทีเอ็มที่เราส่งไปให้ เพื่อทำการคำนวณแปลงพิกัดจาก UTM ไปยัง Geographic จำเป็นต้องรู้พารามิเตอร์ของยูทีเอ็ม เช่น Zone No., Hemisphere อยู่เหนือหรือใต้เส้นศูนย์สูตร (ค่าพิกัดยูทีเอ็มค่าพิกัดเดียวกัน อาจจะมีซ้ำกันหลายๆโซน ดังนั้นถ้าบอกพิกัดยูทีเอ็มต้องบอก Zone No. และ Hemisphere ด้วย)UTMCoordinate property เป็นค่าพิกัดยูทีเอ็มที่ต้องการแปลงค่าพิกัดไปยัง GeographicCompute method เช่นเดียวกันเมื่อ input คือ Ellipsoid, UTMprojection และ UTMCoordinate พร้อมก็เรียกใช้ method นี้ได้เลยเพื่อทำการคำนวณGeoCoordinate property เป็น output คือค่าพิกัด Latitude และ Longitude ที่เราต้องการนั่นเอง
ทิ้งท้ายไว้ตรงนี้ครับ ตอนหน้าเป็นตอนสุดท้าย 🙂
Related
