การเขียนโปรแกรมคำนวณการแปลงค่าพิกัดระหว่าง UTM Grid และ Geographic (Lat/Long) ด้วย Lazarus และ Delphi (ตอนที่ 2)

สูตรที่ใช้ในการคำนวณ

  • ผมอ้างอิงมาจาก 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
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 มี property เก็บค่าของรูปทรงรี เนื่องจากค่าของรูปทรงรีที่คนใช้กันเป็นค่ามาตรฐานผมจึงนำ object ของ class “TEllipsoid” มาสร้างเก็บไว้ใน class ที่สืบทอดจาก TStringList ให้ชื่อว่า TEllipsoidList ทำไมต้อง derive จาก TStringList เนื่องจาก TStringList ออกแบบมาให้เก็บชุดของตัวอักษรแต่สามารถเก็บ Object ได้  มี Index ที่สามารถค้นหาได้ด้วยเมธอด Find ความหมายง่ายๆคือเราเอา Object ของ class “TEllipsoid” มายัดใส่ใน list ของ TEllipsoidList(derive มาจาก TStringList) โดยอาศัยชื่อของทรงรี (EllipsoidName) มาเป็นตัว index เอาไว้ค้นหา Object
  • ผมเพิ่มเมธอด 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 แบบ คือ
    1. TUTMProjection เพื่อช่วยกำหนดตำแหน่งของ UTM Zone เช่น Zone No., CM (Central Meridian), FE (False Easting เท่ากับ 500000 เสมอ), FN(False Northing = 0 ถ้าอยู่ด้านเหนือเส้นศูนย์สูตร และเท่ากับ 10000000 ถ้าอยู่ด้านใต้เส้นศูนย์สูตร), K0 (Scale Factor = 0.9996)
    2. TCoordinate เพื่อเก็บค่าพิกัดซึ่งจะมี 2 ค่าคือ Latitude, Longitude หรือ X,Y (UTM)
  • ถัดลงมาอีกผม declare class ไว้ 2 class คือ TUTM2Geo และ TGeo2UTM มาดูรายละเอียด
    1. 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 ได้
    2. TUTM2Geo ทำหน้าที่คำนวณค่าพิกัดจาก UTM ไปเป็น Geographic
      • Ellipsoid property เช่นเดียวกันกับ class TGeo2UTM
      • UTMProjection property ค่าพารามิเตอร์ของยูทีเอ็มที่เราส่งไปให้ เพื่อทำการคำนวณแปลงพิกัดจาก UTM ไปยัง Geographic จำเป็นต้องรู้พารามิเตอร์ของยูทีเอ็ม เช่น Zone No., Hemisphere อยู่เหนือหรือใต้เส้นศูนย์สูตร (ค่าพิกัดยูทีเอ็มค่าพิกัดเดียวกัน อาจจะมีซ้ำกันหลายๆโซน ดังนั้นถ้าบอกพิกัดยูทีเอ็มต้องบอก Zone No. และ Hemisphere ด้วย)
      • UTMCoordinate property เป็นค่าพิกัดยูทีเอ็มที่ต้องการแปลงค่าพิกัดไปยัง Geographic
      • Compute method เช่นเดียวกันเมื่อ input คือ Ellipsoid, UTMprojection และ UTMCoordinate พร้อมก็เรียกใช้ method นี้ได้เลยเพื่อทำการคำนวณ
      • GeoCoordinate property เป็น output คือค่าพิกัด Latitude และ Longitude ที่เราต้องการนั่นเอง
  • ทิ้งท้ายไว้ตรงนี้ครับ ตอนหน้าเป็นตอนสุดท้าย 🙂

Leave a Reply

Your email address will not be published.