Archive

Archive for ตุลาคม, 2009

แนะนำการใช้ฐานข้อมูล SQLite กับ Lazarus

ข้อดีของ SQLite

  • พัฒนาโดย D. Richard Hipp ด้วยภาษา C จำนวนโค๊ดรวมๆแล้วประมาณสามหมื่นกว่าบรรทัด ซึ่งผู้พัฒนาได้รับคำชมว่าเป็นผู้ที่เข้าใจในวิทยาการด้านคอมพิวเตอร์อย่างลึกซึ้ง
  • สำหรับ SQLite น่าจะเป็นฐานข้อมูลที่นิยมใช้กันมากที่สุดในโลก เนื่องจาก เล็ก เร็ว แรง และที่สำคัญมากคือ เสถียร และข้อดีอีกที่ไม่พูดไม่ได้คือ ฟรีและ cross-platform เป็นฐานข้อมูลที่จัดเก็บในไฟล์เดียว ไม่ต้องมีฝั่ง Server ไม่ต้องตั้งค่าใดๆให้ยุ่งยาก
  • สามารถใช้ภาษาโปรแกรมได้เกือบทุกภาษาในโลกนี้ โดยเฉพาะ Lazarus จะมี component ที่ติดตั้งมาให้พร้อมใช้งานได้เลย

สถานการณ์ที่เหมาะและไม่เหมาะสำหรับ SQLite

  • ถ้าใช้เป็นฐานข้อมูลสำหรับ website ที่คนเข้าไม่มากนักก็ใช้ได้ดี หรือเป็นฐานข้อมูลสำหรับโปรแกรมสำหรับเครื่อง PDA ก็สะดวกเพราะมีขนาดเล็ก สำหรับความคิดของผมแล้วฐานข้อมูลไหนก็ตามที่ใช้ Microsoft Access ใช้ SQLite แทนได้ทั้งนั้น แถมดีกว่าทุกๆอย่าง
  • สถานการณ์ที่ไม่เหมาะได้แก่โปรแกรม Client/Server ที่ระบบข้อมูลมีการส่งผ่านข้อมูลสูงมาก เช่น website ที่มีคนเข้ามาก ควรหันไปใช้ RDBMS ตัวใหญ่ๆจะดีกว่าเช่น MySQL, MS SQL

Download และ ติดตั้ง

  • สำหรับ windows ไปดาวน์โหลดไบนารีได้ที่ http://www.sqlitedll.org/sqlite-3_6_19.zip (ขณะที่เขียน blog นี้เป็น version 3.6.19) เมื่อแตกไฟล์ zip ออกมาจะได้ sqlite3.dll และไฟล์  sqlite3.def  ให้ copy ไปไว้ที่ c:\windows\systems32 หรือ copy ไปไว้ที่โฟลเดอร์ของโปรแกรมที่กำลังพัฒนาด้วย lazarus
  • ถ้าต้องการ tools สำหรับ admin ใช้บน command line เพื่อจัดการฐานข้อมูลก็ดาวน์โหลดได้ที่ http://www.sqlite.org/sqlite-3_6_19.zip
  • สำหรับ Linux ที่ผมใช้อยู่คือ Ubuntu และ PCLinuxOS ใช้คำสั่งเดียวกันคือ
$ sudo apt-get install sqlite

สุดยอด Admin Tools สำหรับ SQLite บนวินโดส์

  • ลองมาหลายอย่างแล้วไม่มี tools ตัวไหนดีเท่ากับตัวนี้ ชื่อโปรแกรมก็คือ SQLite Administrator สมกับที่ผู้พัฒนาเรียกโปรแกรมเขาว่า power tool สามารถออกแบบ สร้าง บริหารฐานข้อมูล SQLite ได้
  • ดาวน์โหลดได้ที่ http://download.orbmu2k.de/files/sqliteadmin.zip
  • โปรแกรมมีขนาดเล็กไม่ต้องติดตั้ง unzip ไปวางไว้ตรงโฟลเดอร์ไหนก็ได้ มารันโปรแกรมกันดู

sqliteadmin1

เริ่มต้นสร้างฐานข้อมูลด้วย SQLite Administrator

  • การสร้างฐานข้อมูล SQLite ใหม่สามารถทำได้หลายวิธี วิธีที่ยากกว่า(แต่ก็ไม่ยากเท่าไหร่) ก็คือสร้างด้วยโค๊ด ที่เราสนใจคือสร้างด้วย Lazarus โดยเรียกใช้ผ่านภาษา sql ส่วนวิธีที่ง่ายกว่า ก้สร้างด้วยโปรแกรม SQLite Administrator ที่แนะนำกันนี้แหละ
  • ตัวอย่างที่จะแนะนำกันได้แก่สร้างฐานข้อมูลเพื่อเก็บรูปทรงรีเพื่อมาใช้แทนยูนิตของ Lazarus ที่ผมเสนอไปหลายตอนแล้วคือ GeoEllipsoids (geoellipsoids.pas) ผมจะเริ่มจาก text file (.csv) ชื่อไฟล์ Ellipsoids.csv ที่เก็บรูปทรงรีดังลิสต์ด้านล่าง (ถ้าผู้อ่านสนใจก็ลอง copy ไปวางที่ text editor แล้วเซฟชื่อเป็น ellipsoids.csv)

    Airy 1830,AA,6377563.396,6356256.909,299.3249646
    Modified Airy,AM,6377340.189,6356034.448,299.3249646
    Australian National,AN,6378160,6356774.719,298.25
    Bessel 1841(Namibia),BN,6377483.865,6356165.383,299.1528128
    Bessel 1841,BR,6377397.155,6356078.963,299.1528128
    Clarke 1866,CC,6378206.4,6356583.8,294.9786982
    Clarke 1880,CD,6378249.145,6356514.87,293.465
    Everest  1830,EA,6377276.345,6356075.413,300.8017
    Everest (E. Malasia & Brunei),EB,6377298.556,6356097.55,300.8017
    Everest 1956,EC,6377301.243,6356100.228,300.8017
    Everest 1969,ED,6377295.664,6356094.668,300.8017
    Everest 1948,EE,6377304.063,6356103.039,300.8017
    Everest (Pakistan),EF,6377309.613,6356109.571,300.8017
    Mod. Fischer 1960(South Asia),FA,6378155,6356773.32,298.3
    Helmert 1906,HE,6378200,6356818.17,298.3
    Hough 1960,HO,6378270,6356794.343,297
    Indonesian 1974,ID,6378160,6356774.504,298.247
    International 1924,IN,6378388,6356911.946,297
    Krassovsky 1940,KA,6378245,6356863.019,298.3
    GRS 80,RF,6378137,6356752.314,298.2572221
    South American 1969,SA,6378160,6356774.719,298.25
    WGS 72,WD,6378135,6356750.52,298.26
    WGS 84,WE,6378137,6356752.314,298.2572236

  • ที่เมนเมนูของโปรแกรมคลิกที่ Database > New… หรือคลิกที่ icon “Create Database” ก็ได้ ป้อนชื่อฐานข้อมูลเป็น Datums.s3db แล้วเซฟไว้

สร้าง Table ใหม่

  • ตอนนี้เราได้ฐานข้อมูลมาแล้วหนึ่งไฟล์ ต่อไปเราจะสร้าง Table เพื่อเก็บรูปทรงรีจะตั้งชื่อ Table นี้ว่า Ellipsoids คลิกขวาที่ Tables > Create Table ดังรูปด้านล่าง

sqliteadmin3

  • จะเห็น dialog ขึ้นมา ป้อนชื่อ table แล้วคลิกที่ Add field ป้อนชื่อ field , field type ดังรูปด้านล่าง

sqliteadmin4

  • ทำการ Add field เข้าไปอีกให้ครบ

sqliteadmin5

  • จากรูปด้านบนผมตั้งให้ field ชื่อ Code ของทรงรี เป็น primary key และไม่ซ้ำกัน ต่อไปเราจะ import หรือปั๊มข้อมูลเข้า table ซึ่งยังว่างๆในตอนนี้

การปั๊มข้อมูลเข้า Table

  • ที่เมนูหลักคลิกที่ Data > Import data… จะเห็น dialog เป็นดังรูปด้านล่าง

sqliteadmin6

  • จากรูปด้านบน จุดที่ 1 คลิกที่ Open File เลือกไฟล์ ellipsoids.csv ที่เราเตรียมไว้ โปรแกรมจะถามว่าใช้อะไรเป็นเครื่องหมายคั่น ให้พิมพ์เครื่องหมายจุลภาค(comma) และถามต่อว่าไฟล์ ellipsoids.csv มีหัวไฟล์ (header column) ตอบ No ดูที่จุดที่ 2 เลือก Target Table เลือกตาราง Ellipsoids จุดที่ 3 ทำการ map column ของไฟล์ให้เข้ากับ field ของฐานข้อมูลดังรูปด้านบน จุดที่ 4 คลิกที่ Import Data เพื่อนำข้อมูลเข้าได้เลย
  • ตรวจสอบข้อมูลได้ ดูรูปด้านล่างประกอบ

sqliteadmin7

  • จากรูปด้านบนตรงแท็ปด้านขวามือจะเห็น SQL Query คลิกเพื่อลองทดสอบการ query พิมพ์ sql  ดังรูปด้านล่าง

sqliteadmin8

  • แล้วคลิกที่เมนูหลัก Query > Execute with Result คลิกไปที่แท็ป Result เพื่อดูผลลัพธ์ดังรูปด้านล่าง

sqliteadmin9

การใช้เครื่องมือ Adimin ที่เป็น command line

  • สำหรับความรู้สึกผมแล้วการใช้ command line ในบางครั้งกลับสะดวกกว่าโปรแกรมที่ใช้ GUI เสียอีก ตัวอย่างก็หนีไม่พ้นใน linux การใช้ command line เป็นเรื่องสะดวกมากๆ เช่นคำสั่ง ls, dmesg, df, env, lspci, lsmod อะไรประมาณนี้ การใช้เครื่องมือ admin ที่แนะนำให้ดาวน์โหลดไปตอนต้นคือ sqlite3.exe ผมจะแนะนำการใช้คร่าวๆ ผม copy ไปไว้ที่โฟลเดอร์ c:\sqlite3 ถ้าใช้ windows กดคีย์บอร์ด ALT + R พิมพ์คำสั่ง cmd กด Enter แล้วเรียกใช้คำสั่งดังรูปด้านล่างเพื่อเรียกโปรแกรม
  • cmd1ผมเปิดฐานข้อมูลด้วยคำสั่ง sqlite3 e:\sourcecodes\lazarus\sqliteproj1\datums.s3db ถ้าสงสัยหรือจำคำสั่งไม่ได้ก็ให้พิมพ์ .help เข้าไป
  • ผมใช้คำสั่ง .databases เพื่อดูว่า database อะไรที่เราใช้งานอยู่ ตามด้วย .header ON เพื่อให้แสดงคอลัมภ์ชื่อฟิลด์ และ .show เพื่อแสดงพารามิเตอร์ แล้วลอง query  ดูด้วย syntax “select” ดังรูปด้านล่าง
  • cmd2จากรูปด้านบนจะเห็นผลลัพธ์ของ sql  ด้วย syntax “select” ที่จริงยังมีคำสั่งอีกมากที่ทำได้ในสภาวะ command line แต่ตอนนี้ขอพอแค่นี้ก่อน ตอนหน้าจะลงลึกไปอีกนิดว่าจะนำ SQLite มาใช้กับ Lazarus อย่างไร
Categories: Linux, Programming, Windows ป้ายกำกับ:, , ,

การเขียนโปรแกรมเพื่อคำนวณระยะทางและอะซิมัท (Distance/Azimuth) บน Ellipsoid ด้วย Lazarus (ตอนที่ 2)

  • ตอนที่แล้วได้แนะนำสูตรที่จะใช้ในการคำนวณและแสดงยูนิต GeodesicCompute พร้อมทั้งยูนิต GeoEllipsoids ที่เคยแสดงไปแล้วเรื่องการแปลงค่าพิกัดระหว่าง UTM และ Geographic
  • เปิด Lazarus คลิกที่เมนเมนู Project > New Project… คลิกเลือก Application คลิก OK

ตั้งค่า Compiler Options ก่อน

  • ผมขอย้ำอีกครั้งว่าเมื่อสร้าง Project ใหม่แล้วสิ่งแรกที่ต้องทำคือตั้ง Compiler Options ให้ตรงกับ Platform และ Widget ที่เราใช้อยู่ ที่เมนู Project > Compiler options… ที่ path เลือก LCL Widget เป็น win32/win64 ดังรูปด้านล่าง
      geodesicsb1
  • ที่ code คลิกเลือก Target OS เป็น win32/win64 และ Target CPU Family เป็น i386 ส่วน Target processor เลือกเป็น Default ดังรูปด้านล่าง

geodesicsb2

เพิ่มยูนิตเข้าไปใน Project

  • จากที่ได้กล่าวไปข้างต้นเรามี 2 ยูนิตที่เตรียมไว้แล้วคือ GeodesicCompute (geodesiccompute.pas) และ GeoEllipsoids (geoellipsoids.pas) ที่ Project Inspector คลิกที่เครื่องหมายบวกเพื่อเพิ่มยูนิต คลิกที่แท็ป Add files จากนั้นคลิกที่ Browse เพื่อเลือกไฟล์ที่เราเตรียมไว้

projectInspector

  • เมื่อเลือกไฟล์ได้แล้วจะเห็นหน้าตาของ Project Inspector ดังรูปด้านล่าง

projectInspector2

  • คลิกที่เมนู Project > Save Project As… ตั้งชื่อเป็น Geodesics จะเห็นหน้าตาของ Lazarus และ project Geodesics ดังรูปด้านล่าง geodesicsb4
  • ต่อไปจะแปะพวก object ทั้งหลายเฃ่น Label, Edit, ComboBox, Button ลงไปและ ตั้งฃื่อ (Name) ดังรูปด้านล่าง

geodesicsb5

  • ส่วน Label และ GroupBox ใ้ส่ชื่อดังรูปด้านล่าง

geodesicsb7

โค๊ดของฟอร์ม Unit1

unit Unit1;

{$mode objfpc}{$H+}

interface

uses
 Classes, SysUtils, FileUtil, LResources, Forms, Controls, Graphics, Dialogs,
 StdCtrls, ExtCtrls, GeoEllipsoids, GeodesicCompute;

type

 { TForm1 }

 TForm1 = class(TForm)
   cmdEx1: TButton;
   cmdCompute: TButton;
   cmdEx2: TButton;
   cmdClose: TButton;
   cmbEllipsoids: TComboBox;
   grbIn: TGroupBox;
   grbOut: TGroupBox;
   labin2: TLabel;
   labin1: TLabel;
   labout3: TLabel;
   labout1: TLabel;
   Label8: TLabel;
   labout2: TLabel;
   radCT: TRadioGroup;
   txtIn1: TEdit;
   txtOut1: TEdit;
   txtLong1: TEdit;
   txtLat1: TEdit;
   GroupBox1: TGroupBox;
   Label1: TLabel;
   Label2: TLabel;
   txtIn2: TEdit;
   txtOut3: TEdit;
   txtOut2: TEdit;
   procedure cmdCloseClick(Sender: TObject);
   procedure cmdComputeClick(Sender: TObject);
   procedure cmdEx1Click(Sender: TObject);
   procedure cmdEx2Click(Sender: TObject);
   procedure FormCreate(Sender: TObject);
   procedure radCTClick(Sender: TObject);
 private
 { private declarations }
   FEllipses : TEllipsoidList;
 public
{ public declarations }
 end;

var
  Form1: TForm1;

implementation

{ TForm1 }

function Coordinate(const X, Y : double) : TCoordinate;
begin
  result.X := X;
  result.Y := Y;
end;

procedure TForm1.cmdCloseClick(Sender: TObject);
begin
  Close;
end;

procedure TForm1.cmdComputeClick(Sender: TObject);
var
  inverse : TGeodesicInverse;
  direct : TGeodesicDirect;
  ell : TEllipsoid;
  x1, y1, x2, y2, az, ds : double;
  geocoor : TCoordinate;
begin
  ell := FEllipses.FindEx(cmbEllipsoids.Text);
  if (radCT.ItemIndex = 0) then
  begin
    x1 := strtofloat(txtlong1.Text);
    y1 := strtofloat(txtlat1.Text);
    x2 := strtofloat(txtIn2.Text);
    y2 := strtofloat(txtIn1.Text);
    try
      inverse := TGeodesicInverse.Create;
      inverse.Ellipsoid := ell;
      inverse.GeoCoordinate1 := Coordinate(x1, y1);
      inverse.GeoCoordinate2 := Coordinate(x2, y2);
      inverse.Compute;
      txtOut1.Text := format('%11.8f', [inverse.InitialAzimuth]);
      txtOut2.Text := format('%11.8f', [inverse.FinalAzimuth]);
      txtOut3.Text := format('%10.3f m.', [inverse.Distance]);
    finally
      inverse.Free;
    end;
  end
  else if (radCT.ItemIndex = 1) then //Direct
  begin
    x1 := strtofloat(txtlong1.Text);
    y1 := strtofloat(txtlat1.Text);
    ds := strtofloat(txtIn2.Text);
    az := strtofloat(txtIn1.Text);
    try
      direct := TGeodesicDirect.Create;
      direct.Ellipsoid := ell;
      direct.GeoCoordinate1 := Coordinate(x1, y1);
      direct.Azimuth := az;
      direct.Distance := ds;
      direct.Compute;
      x2 := direct.GeoCoordinate2.X;
      y2 := direct.GeoCoordinate2.Y;
      txtOut1.Text := format('%11.8f', [y2]);
      txtOut2.Text := format('%11.8f', [x2]);
    finally
      direct.Free;
    end;
  end;
end;

procedure TForm1.cmdEx1Click(Sender: TObject);
begin
  radCT.ItemIndex := 0; //Geodesic Inverse
  txtLat1.Text := '53.150555556';
  txtLong1.Text :='-1.844444444';
  txtIn1.Text := '52.205277778';
  txtIn2.Text :='-0.142500000';
  txtOut1.Text := '';
  txtOut2.Text := '';
  txtOut3.Text := '';
end;

procedure TForm1.cmdEx2Click(Sender: TObject);
begin
  radCT.ItemIndex := 1; //Geodesic Direct
  txtLat1.Text := '-37.9510334';
  txtLong1.Text := '144.424868';
  txtIn1.Text :='306.8681583';
  txtIn2.Text :='54972.271';
  txtOut1.Text := '';
  txtOut2.Text := '';
  txtOut3.Text := '';
end;

procedure TForm1.FormCreate(Sender: TObject);
var
 i : integer;
begin
  FEllipses := TEllipsoidList.Create;
  for i := 0 to FEllipses.Count - 1 do
  begin
    cmbEllipsoids.Items.Add(TEllipsoid(FEllipses.Objects[i]).EllipsoidName);
  end;
  cmbEllipsoids.ItemIndex := FEllipses.Count - 1; //Default with WGS84
end;

procedure TForm1.radCTClick(Sender: TObject);
begin
  if (radCT.ItemIndex = 0) then //Inverse.
  begin
    grbin.Caption := 'Point 2';
    labin1.Caption := 'Latitude : ';
    labin2.Caption := 'Longitude : ';
    grbout.Caption := 'Azimuth && Distance';
    labout1.Caption := 'Initial Azimuth : ';
    labout2.Caption := 'Final Azimuth : ';
    labout3.Show;
    txtout3.Show;
    labout3.Caption := 'Distance : ';
  end
  else if (radCT.ItemIndex = 1) then
  begin
    grbin.Caption := 'Azimuth && Distance';
    labin1.Caption := 'Azimuth : ';
    labin2.Caption := 'Distance : ';
    grbout.Caption := 'Point 2';
    labout1.Caption := 'Latitude : ';
    labout2.Caption := 'Longitude : ';
    labout3.Hide;
    txtout3.Hide;
  end;
end;

initialization
 {$I unit1.lrs}

end.

อธิบายโค๊ด

  • ที่นี้มาลองไล่ดูโค๊ดเฉพาะบางส่วนที่สำคัญ เริ่มแรกตรงเมธอด FormCreate
<pre>procedure TForm1.FormCreate(Sender: TObject);
var
  i : integer;
begin
  FEllipses := TEllipsoidList.Create;
  for i := 0 to FEllipses.Count - 1 do
  begin
    cmbEllipsoids.Items.Add(TEllipsoid(FEllipses.Objects[i]).EllipsoidName);
  end;
  cmbEllipsoids.ItemIndex := FEllipses.Count - 1; //Default with WGS84
end;</pre>
  • เ่ริ่มจากการสร้าง object สำหรับเก็บทรงรี TEllipsoidList นั้นออกแบบมาคล้ายๆกับ Collection เมื่อสร้าง Object เรียบร้อยแล้วก็นำมาใ่ส่ให้ ComboBox ชื่อ cmbEllipsoids เ้สร็จแล้วให้สถานะของทรงรีเป็น WGS84 (ทรงรี WGS84 จะอยู่ท้ายสุดใน list)
  • ต่อไปผมสร้าง Event เมื่อคลิกด้วยเมาส์เพื่อโยงให้ radio group box ชื่อ radCT (property Items ผมใส่ string ให้ 2 string คือ Geodesic Inverse และ Geodesic Direct ตามลำดับ) ตัว radCT จะทำหน้าที่ให้ผู้ใช้โปรแกรมเลือกว่าจะคำนวณแบบไหน เมื่อเลือกแบบใดแบบหนึ่งโปรแกรมจะเปลี่ยน label ให้สอดคล้อง input และ output
<pre>procedure TForm1.radCTClick(Sender: TObject);
begin
  if (radCT.ItemIndex = 0) then //Inverse.
  begin
    grbin.Caption := 'Point 2';
    labin1.Caption := 'Latitude : ';
    labin2.Caption := 'Longitude : ';
    grbout.Caption := 'Azimuth && Distance';
    labout1.Caption := 'Initial Azimuth : ';
    labout2.Caption := 'Final Azimuth : ';
    labout3.Show;
    txtout3.Show;
    labout3.Caption := 'Distance : ';
  end
  else if (radCT.ItemIndex = 1) then
  begin
    grbin.Caption := 'Azimuth && Distance';
    labin1.Caption := 'Azimuth : ';
    labin2.Caption := 'Distance : ';
    grbout.Caption := 'Point 2';
    labout1.Caption := 'Latitude : ';
    labout2.Caption := 'Longitude : ';
    labout3.Hide;
    txtout3.Hide;
  end;
end;</pre>
  • เพื่อให้การใช้โปรแกรมนั้นง่าย ผมจะใส่ตัวอย่างข้อมูสำหรับป้อนเข้า ไว้สองตัวอย่างคือคำนวณแบบ Inverse (Example 1) และคำนวณแบบ Direct (Example 2) เมื่อผู่ใช้ทำการคลิกที่ปุ่ม Example 1 จะเรียก Event ที่สร้างไว้คือ cmdEx1Click ส่วนปุ่ม Example 2 จะเรียก event cmdEx2Click
<pre>
<pre>procedure TForm1.cmdEx1Click(Sender: TObject);
begin
  radCT.ItemIndex := 0; //Geodesic Inverse
  txtLat1.Text := '53.150555556';
  txtLong1.Text :='-1.844444444';
  txtIn1.Text := '52.205277778';
  txtIn2.Text :='-0.142500000';
  txtOut1.Text := '';
  txtOut2.Text := '';
  txtOut3.Text := '';
end;

procedure TForm1.cmdEx2Click(Sender: TObject);
begin
  radCT.ItemIndex := 1; //Geodesic Direct
  txtLat1.Text := '-37.9510334';
  txtLong1.Text := '144.424868';
  txtIn1.Text :='306.8681583';
  txtIn2.Text :='54972.271';
  txtOut1.Text := '';
  txtOut2.Text := '';
  txtOut3.Text := '';
end;</pre>
</pre>
  • สุดท้ายก็มาดูที่เมธอดสำคัญคือ cmdComputeClick สร้างไว้เพื่อดัก event เมื่อผู้ใช้คลิกที่ปุ่ม Compute
<pre>procedure TForm1.cmdComputeClick(Sender: TObject);
var
  inverse : TGeodesicInverse;
  direct : TGeodesicDirect;
  ell : TEllipsoid;
  x1, y1, x2, y2, az, ds : double;
  geocoor : TCoordinate;
begin
  ell := FEllipses.FindEx(cmbEllipsoids.Text);
  if (radCT.ItemIndex = 0) then
  begin
    x1 := strtofloat(txtlong1.Text);
    y1 := strtofloat(txtlat1.Text);
    x2 := strtofloat(txtIn2.Text);
    y2 := strtofloat(txtIn1.Text);
    try
      inverse := TGeodesicInverse.Create;
      inverse.Ellipsoid := ell;
      inverse.GeoCoordinate1 := Coordinate(x1, y1);
      inverse.GeoCoordinate2 := Coordinate(x2, y2);
      inverse.Compute;
      txtOut1.Text := format('%11.8f', [inverse.InitialAzimuth]);
      txtOut2.Text := format('%11.8f', [inverse.FinalAzimuth]);
      txtOut3.Text := format('%10.3f m.', [inverse.Distance]);
    finally
      inverse.Free;
    end;
  end
  else if (radCT.ItemIndex = 1) then //Direct
  begin
    x1 := strtofloat(txtlong1.Text);
    y1 := strtofloat(txtlat1.Text);
    ds := strtofloat(txtIn2.Text);
    az := strtofloat(txtIn1.Text);
    try
      direct := TGeodesicDirect.Create;
      direct.Ellipsoid := ell;
      direct.GeoCoordinate1 := Coordinate(x1, y1);
      direct.Azimuth := az;
      direct.Distance := ds;
      direct.Compute;
      x2 := direct.GeoCoordinate2.X;
      y2 := direct.GeoCoordinate2.Y;
      txtOut1.Text := format('%11.8f', [y2]);
      txtOut2.Text := format('%11.8f', [x2]);
    finally
      direct.Free;
    end;
  end;
end;
</pre>
  • จากโค๊ดด้านบน เมื่อผู้ใช้เลือกแบบการคำนวณโดยการคลิกที่ radio group โปรแกรมเราจะเลือกการคำนวณให้ที่เหมาะสม โดยจะมี 2 ตัวแปรสำคํญคือ direct เรา assign จาก class ชื่อ TGeodesicDirect และตัวแปรอีกตัวคือ inverse โดย derive มาจาก class ชื่อ TGeodesicInverse ลองดู source code ก็ไม่มีอะไรยากเลย เริ่มต้นจากการสร้าง object => inverse หรือ direct จากนั้นก็เอาทรงรีมา assign ให้จากผู้ใช้คลิกเลือกทรงรีที่ ComboBox ต่อไปก็ input ค่าพิกัด Geographic ให้ 2 จุด (ถ้าคำนวณ inverse) หรือป้อนค่าพิกัดให้ 1 จุด แล้วป้อน ระยะทาง และ อะซิมัท จากนั้นจะเรียกเมธอด Compute และส่งค่าผลลัพธ์มาทาง property

การคำนวณแบบ Inverse

  • เมื่อทุกสิ่งทุกอย่างพร้อมกดปุ่มคีย์บอร์ด F9 เพื่อให้ Lazarus ทำการ compile และ build เราตั้งชื่อ project เราว่า Geodesics เมื่อ build แล้วจะได้ไฟล์ Geodesics.exe ส่วนใน Linux จะไม่มี extension ให้ เวลารันถ้าใช้ command ใน linux จะเป็นดังนี้ $./geodesics (พิมพ์จุดและเครื่องหมาย fore slash แล้วตามด้วยชื่อโปรแกรม) ตัวอย่างเมื่อรันโปรแกรม คลิกที่ปุ่ม Example 1 แล้วคลิก Compute เพื่อทำการคำนวณ จะได้ผลลัพธ์ดังรูปด้านล่าง

geodesicsb8

  • สังเกตว่าที่เราทำการคำนวณคือระยะทางและอะซิมัทไปตามเส้น Geodesic ซึ่งเป็นเส้นโค้งอะซิมัทจะค่อยๆเปลี่ยนไปตามเส้นโค้ง ดังน้น Initial Azimuth จึงเป็นค่าอะซิมัทเริ่มต้น และ Final Azimuth คือค่าอะซิมัทสุดท้ายที่เส้น geodesic ไปแตะจุดที่สองพอดี

การคำนวณแบบ Direct

  • รันโปรแกรม คลิกที่ปุ่ม Example 2 แล้วคลิก Compute เพื่อทำการคำนวณ จะได้ผลลัพธ์ดังรูปด้านล่าง

geodesicsb9

ข้อจำกัดของโปรแกรม

  • เนื่องจากเป็นโปรแกรมตัวอย่างที่ให้ดูกันง่ายๆ จึงไม่มีการดักข้อมูล error จากการพิมพ์ใส่ข้อมูล ที่ EditBox เช่นถ้าผู้ใช้ป้อนแทนที่จะเป็นตัวเลขแต่ไปใส่อักขระอื่นแทนก็จะ error และไม่มีการดักจับค่าพิกัด 2 จุดจะต้องไม่เท่ากัน ถ้าเท่ากันก็จะ error เช่นเดียวกันขณะรันโปรแกรม
  • ก็ขอจบกันไว้เพียงเท่านี้ ตอนหน้าจะมาว่าเรื่องฐานข้อมูลฉบับกระเป๋าเล้กพริกขี้หนู SQLite ส่วนในตอนต่อๆไปอาจจะเป็นฐานข้อมูลเจ้านกไฟ FireBird

Download sourcecode

  • sourcecode ของโปรแกรมตอนนี้สามารถ download ได้ที่นี่ GeodesicProj.zip

เทคนิคการใช้โปรแกรมแปลงค่าพิกัด GeoCalc

  • ใน Blog ของ WordPress ที่ผมเขียนอยู่หา Theme ที่ถูกใจยาก ไม่ใช่เรื่องความสวยงามแต่เป็นขนาดของคอลัมภ์ ส่วนใหญ่จะแคบมากพอวางรูปแล้วจะต้องย่อมากๆ และอีกเรื่องคือวาง source code แล้วจะโดนปัดบรรทัดลงทำให้ดู code ยาก งง ตอนนี้ก็เปลี่ยนไปเปลี่ยนมายังไม่ถูกใจสักที
  • ผมเห็นคน search เข้ามาใน blog มากเรื่องการแปลงพิกัด บางคนมีแต่พิกัดยูทีเอ็ม ต้องการค่าพิกัดเป็น Lat/Long เพื่อไปจุดลง Google Earth แต่ไม่รู้จะใช้โปรแกรมอะไร คือโปรแกรมในอินเทอร์เน็ตมันมากเกินจนสับสน ไม่รู้จะเลือกอะไร สำหรับผมแล้ว Coordinate Calculator ของ Trimble คือสุดยอดของโปรแกรมแปลงค่าพิกัดพวกนี้ รองลงมาก็ได้แก่ GeoCalc ของ GeoComp ซึ่งเป็นของฟรี ส่วน Coordinate Calculator ของ Trimble จะมาพร้อมกับ Terramodel (ที่ทีมงานผมใช้แทน Autodesk Land Desktop มานานแล้ว) หรือ HydroPro Navigation (ทีมงานผมใช้อยู่แต่ค่อนข้างหน่อมแน้มเมื่อเทียบกับคู่แข่งคือ HyPack) โปรแกรมสำหรับงาน Hydrographic Survey หรือ Trimble Geomatic Office สำหรับงาน GPS ขั้นสูง
  • GeoCalc ผมเคยแนะนำโปรแกรมนี้ไปก่อนหน้านี้ เป็นของฟรี ดาวน์โหลดได้ที่ http://www.geocomp.com.au/geocalc/gcalc420.exe ถึงจะเป็นโปรแกรมเก่าไปหน่อยแต่ก็ใช้ดี โปรแกรมมีขนาดเล้ก เรียบง่าย แต่ก่อนจะใช้งานต้องปรับจูนกันนิดหน่อย ให้เหมาะกับ Datum และทรงรีที่เราใช้งานกันอยู่ เมื่อดาวน์โหลดมาก็ทำการติดตั้งให้เป็นที่เรียบร้อย

แก้ค่าสัณฐานทรงรี Everest 1830 (Indian 1975 datum) ให้ถูกต้อง

  • ประเทศไทยก่อนหน้านี้ใช้ Indian 1975 เป็น Datum มีทรงรีเรียกว่า Everest 1830 มีค่า semi-major axis (a) =  6,377,276.345 m. ค่า flattening (1/ f) = 300.8017 หลังๆมาเราหันมาใช้ทรงรี WGS84 กันมากขึ้น ดังนั้นระหว่างระบบพิกัดเดิม (Indian 1975) กับระบบพิกัดใหม่ (WGS84) ก็อาจจะต้องมีโปรแกรมที่ทำหน้าที่คำนวณแปลงพิกัดระหว่าง datum
  • โปรแกรม GeoCalc มี mistake อยู่เล็กน้อยคือค่าพารามิเตอร์ของทรงรี Everest ไม่ถูกต้อง เมื่อติดตั้งเรียบร้อยลองรันโปรแกรมดู เราจะทำการแก้ไขค่าพารามิเตอร์ของทรงรี Everest  1830 ให้ถูกต้องเสียก่อน
เมนูแก้ไขค่าพารามิเตอร์ทรงรี

เมนูแก้ไขค่าพารามิเตอร์ทรงรี

  • จากรูปด้านบนคลิกเมนู Edit > Spherod Definitions จะเห็นไดอะล็อก ขึ้นมาจะทำการแก้ไขค่าพารามิเตอร์ของทรงรี Everest 1830
ตั้งค่าพารามิเตอร์ Everest 1830 ให้ถูกต้อง

ตั้งค่าพารามิเตอร์ Everest 1830 ให้ถูกต้อง

  • จากรูปด้านบนที่ combo box “Edit/Select Spheroid” คลิกไปหา “Thailand/Vietnam (Indian Everest)” ที่นี้ถ้าติดตั้งโปรแกรมครั้งแรกจะเห็นค่า Semi-Major Axis = 6,377,267.345 ซึ่งผิด ดังนั้นจุดที่ 1 แก้ไขค่าเป็น 6,377,276.345 (เลข 7 กับเลข 6 สลับตำแหน่งกัน) จุดที่ 2 ตั้งค่า Translation X, Y และ Z เป็น 206, 837 และ 295 ตามลำดับ
  • ต่อไปจะเพิ่ม Datum ที่ประเทศไทยใช้งานอยู่ เหมือนเดิมที่เมนู Edit > Map System Definitions
เมนูแก้ไข Coordinate Systemsเมนูแก้ไข Coordinate Systems
  • โปรแกรม Geocalc เรียกว่า Map System ท่านผู้อ่านอาจจะพบคำว่า Coordinate Systems มากกว่า มีความหมายคือระบบพิกัดนั่นเอง

เพิ่มระบบพิกัด Thailand UTM 47N Indian 1975

เพิ่ม Coordinate Systems ให้กับประเทศไทย

เพิ่ม Coordinate Systems ให้กับประเทศไทย

  • โปรแกรมนี้มีทรงรี Everest 1830 อยู่แต่ระบบพิกัดของประเทศไทยไม่ได้เตรียมไว้ แต่ข้อเสียของโปรแกรมคือไม่สามารถเพิ่ม (Add) เข้าไปได้ เราจะใช้วิธีหักเอาด้วยกำลัง คือเขียนทับ Coordinate Systems ของประเทศอื่นๆ ที่ผู้ใช้ไม่คิดว่าจะได้ใช้ระบบพิกัดของเขา ตอนนี้ผมเล็งที่ช่อง Selected Mapping System ดังรูปด้านบนผมเลือกเอาบรรทัดที่ห้า พิมพ์ THAILAND UTM47N INDIAN 1975 ทับของเดิมเข้าไป

เพิ่มระบบพิกัด Indian 1975 ให้กับประเทศไทย

เพิ่มระบบพิกัด Indian 1975 ให้กับประเทศไทย

  • จากรูปด้านบน จุดที่ 1 พิมพ์ THAILAND UTM47N INDIAN 1975 จุดที่ 2 เลือกทรงรีเป็น Thailand/Vietnam(Indian Everest) จุดที่ 3 เลือก Projection Type เป็น Transverse Mercator จุดที่ 4 ตั้งค่า Origin Parameters ตามที่เห็น จุดที่ 5 ตั้ง Central Meridian เป็น Fixed และข้างล่างติดกันเลือก Zone 0 Longitude เท่ากับ 47 (ก็คือ UTM Zone 47)

เพิ่มระบบพิกัด Thailand UTM 47N WGS84

  • ตอนนี้เรามีระบบพิกัดของไทยที่ใช้ทรงรี Everest 1830 เป็นที่เรียบร้อยที่ผมให้ชื่อไว้ THAILAND UTM 47 N INDIAN 1975 ที่นี้เรามาสร้างระบบพิกัดของไทยที่ใช้ทรงรี WGS84 กันบ้างวิธีการก็เหมือนเดิมผมเลือกเอาบรรทัดที่หก ตรง ComboBox “Selected Mapping System” พิมพ์ THAILAND UTM 47N WGS84 ตั้งค่าต่างๆดังรูปด้านล่าง
เพิ่มระบบพิกัด UTM 47N บนทรงรี WGS84

เพิ่มระบบพิกัด UTM 47N บนทรงรี WGS84

คำนวณการแปลงพิกัดจาก Indian 1975 ไป WGS84

  • ที่เมนูหลักคลิก Process > Manual Input
เลือกเมนูคำนวณ

เลือกเมนูคำนวณ

  • โจทย์ของผมวันนี้คือเรามีค่าพิกัด UTM อยู่ในระบบ Indian 1975 มีค่า Easting(E) = 391024.838 ค่า Northing(N) = 1576384.958 ต้องการแปลงเป็นค่า Lat/Long ในระบบพิกัด WGS84 ดูรูปด้านล่าง
การคำนวณการแปลงค่าพิกัดจาก Indian 1975 ไปยัง WGS84

การคำนวณการแปลงค่าพิกัดจาก Indian 1975 ไปยัง WGS84

  • จากรูปด้านบน จุดที่ 1 ที่ช่อง From และ To ตั้งค่าดังรูป จะสอดคล้องกับ Group Box ด้านซ้ายจะเป็นช่องนำเข้าข้อมูล Point in From Map System ด้านขวามือจะมี Point in to Map System จุดที่ 2 ป้อนค่าพิกัดที่ตั้งโจทย์ไว้ จุดที่ 3 คือผลลัพธ์ค่าพิกัดในรูป Lat/Long ที่เราต้องการ โปรดสังเกตด้านล่างที่ขีดเส้นใต้เป็นฟอร์แมตของค่าพิกัด Lat/Long เป็น DDD.MMSSssss นอกจากนั้น Group Box ด้านซ้ายจะเห็นค่าพิกัด Lat/Long จะเป็นค่าพิกัดแลตติจูดและลองจิจูดบน Indian 1975 หรือจะเรียกว่าบนทรงรี Everest 1830 ก็ได้ และเช่นเดียวกัน Group Box ด้านขวาจะเห็น Easting และ Northing ที่ผมไม่ได้พูดถึงเป็นค่าพิกัด UTM บนทรงรี WGS84 คลิกที่ปุ่ม “Calculate” เราคำนวณครั้งเดียวจะได้ค่ามาหมดเลย ขี้นอยู่กับว่าต้องการค่าไหนไปใช้
  • จากรูปด้านบนเราได้ค่าพิกัด Lat = 14.15349881 Long = 97.59120712 ถ้าดูรูปฟอร์แมต DDD.MMSSssss กระจายค่า latitude 14.15349881 ได้ 14 องศา 15 ลิปดา 34.9881 ฟิลิปดา แต่ถ้าคิดเป็น degree (ฟอร์แมต DDD.dddddddd) จะได้ 14+15/60+34.9881/3600 = 14.25971891
  • ส่วนค่า longitude 97.59120712 เขียนกระจายได้ 97 องศา 59 ลิปดา 7.12 ฟิลิปดา เขียนเป็น degree ได้ = 97+59/60+7.12/3600 = 97.985311

ทดสอบการจุดค่าพิกัดบน Google Earth

  • ตอนนี้เรามีพิกัดอยู่สองแบบของฟอร์แม็ต แบบแรกคือ DDD.MMSSssss ดังที่โปรแกรม GeoCalc แสดงผลลัพธ์ให้เรา ส่วนอีกแบบคือเป็นค่า degree เรามาดูแบบแรกกัน Lat = 14.15349881 Long = 97.59120712 (ค่า Lat = 14 องศา 15 ลิปดา 34.9881 ฟิลิปดา ค่า Long = 97 องศา 59 ลิปดา 7.12 ฟิลิปดา) เขียนกระจายใหม่เพื่อนำไปใช้กับ Google Earth ได้ 14 15 349881N,97 59 07.12E ข้อสำคัญคือเราต้องเติมค่า N ถ้าอยู่ด้านเหนือของเส้นศูนย์สูตร และลองจิจูดใส่ค่า E ไปด้านท้ายค่าลองจิจูดให้เพราะเราอยู่ซีกโลกด้านตะวันออก และไม่ลืมที่จะคั่นแลตติจูดกับลองจิจูดด้วยเครื่องหมายคอมม่า (,)
  • ที่โปรแกรม Google Earth เปิดที่เมนูคลิก Tools > Options… ตั้งค่าตรงประเภทของมุมให้เป็น Degrees,Minutes,Seconds ตามรูปด้านล่าง
ตั้งค่าชนิดของมุมให้ Google Earth

ตั้งค่าชนิดของมุมให้ Google Earth

  • จากรูปด้านบนคลิก OK แล้ว Copy ค่า 14 15 34.9881N,97 59 7.12E ไปวางไว้ตรง Search (จุดที่ 1 ของรูป) ของ Google Earth กด Enter ดูผลลัพธ์รูปด้านล่างตรงจุดที่ 2

ผลลัพธ์ของการจุดค่าพิกัดบน Google Earth

ผลลัพธ์ของการจุดค่าพิกัดบน Google Earth

  • ถ้าจะป้อนค่า Lat/Long บน Google Earth เป็นค่า degree ต้องไปเปลี่ยน ที่เมนู Tools > Options… โดยคลิกรูปแบบมุม (Show Lat/Long) ที่ Decimal Degree แล้วค่าพิกัดให้เอาค่า degree ที่ผมกล่าวไว้ข้างต้นคือ 14.25971891N,97.985311E ไปวางที่ Search ของ Google Earth จะให้ผลเช่นเดียวกัน
  • เสริมอีกนิดหน่อยครับ ตรง Search ของ Google Earth ที่เราป้อนค่าที่ต้องการค้นหา นอกจากค่า Lat/Long แล้ว โปรแกรมจะเอาประโยคที่เราป้อนส่งไป process ถ้าเป็นที่อยู่ต่างประเทศเช่น อเมริกาหรือยุโรป เช่นถ้าเราป้อนที่อยู่จากบ้านเลขที่ ถนน เขต เมือง โปรแกรมจะทำการ GeoCode เพื่อคำนวณหา Lat/Long ให้ เพราะที่ๆอยู่ประเทศเหล่านี้เขาวางผังมาดี แต่บ้านเราทำไม่ได้ครับ ขนาดบุรุษไปรษณีย์ไปบ่อยยังหลงเลยครับ

ข้อดีและข้อจำกัดของ GeoCalc

  • ข้อจำกัดของโปรแกรมคือผู้อ่านจะเห็นว่ามุมทั้ง input และ output ป้อนได้ฟอร์แม็ตเดียวคือ DDD.MMSSssss น่าจะสนับสนุนแบบ degree ข้อจำกัดข้อที่สองคือไม่สนับสนุนการคำนวณ Geoid (ทอนความสูงจากบนทรงรีมาเป็น Orthometric Height หรือ MSL) ข้อดีคือฟรี และเล็กๆ ขนโปรแกรม (Copy โฟลเดอร์ของ Geocalc ก็พอ) ใส่ Thumb drive ก็คำนวณได้

การเขียนโปรแกรมเพื่อคำนวณระยะทางและอะซิมัท (Distance/Azimuth) บน Ellipsoid ด้วย Lazarus (ตอนที่ 1)

  • ตอนก่อนหน้านี้ ผมเขียนโปรแกรมแปลงพิกัดระหว่าง UTM และ Geographic (Lat/Long) และและถ้าไม่เขียนการหาระยะทางและอะซิมัท (เมื่อกำหนดจุด Latitude, Longitude ให้สองจุด) ก็ดูจะขาดอะไรไปอย่าง

Model ที่ใช้ในการคำนวณ

  • สัณฐานหรือรูปทรงที่ใช้แทนโลก ใช้กันอยู่ 2 แบบ คือ ทรงกลม(Spherical)และทรงรี(Ellipsoid) แต่ทรงรีจะใกล้เคียงกับกับโลกเรามากกว่า ดังนั้นการหาระยะทางและทิศทาง(อะซิมัท) บนทรงรีจะให้ความละเอียดถูกต้องมากกว่า
รูปทรงที่ใช้แทนสัณฐานของโลก

รูปทรงที่ใช้แทนสัณฐานของโลก (ภาพอ้างอิงจาก http://www.codeguru.com/cpp/cpp/algorithms/article.php/c5115/)

  • การคำนวณระยะทางและอะซิมัทบนทรงกลมนั้นง่ายกว่าการคำนวณบนทรงรี การคำนวณระยะทางบนทรงกลมเช่นต้องการทราบระยะทางจากจุด A ไปจุด B เราจะใช้ plane คือแผ่นเรียบตัดผ่านจุด A และจุด B และที่สำคัญคือต้องผ่านจุดศูนย์กลางทรงกลมด้วย เส้นที่เกิดจากแผ่น plane มาตัดกับทรงกลมจะเรียกว่า  Great Circle ระยะที่สั้นที่สุดบนผิวทรงกลมก็คือระยะที่ไปตาม Great circle นี่เอง
Great circle

Great circle อ้างอิงจาก http://www.codeguru.com/cpp/cpp/algorithms/article.php/c5115/

  • ส่วนทรงรีที่สัณฐานใกล้เคียงกับโลกมากกว่า จะยุบลงที่ขั้วโลก โป่งออกที่เส้นศูนย์สูตร เราจะหาระยะทางและอะซิมัทโดยอิงสัณฐานของทรงรีโดยใช้สูตรอันลือชื่อของ Thaddeus Vincenty

สูตรที่ใช้อ้างอิงในการคำนวณ (Vincenty’s formulae)

  • Thaddeus Vincenty เกิดที่โปแลนด์ อพยพเข้าสหรัฐอเมริกาเมื่อสงครามโลกครั้งที่สอง ทำงานให้กับ U.S. Air Force และ  National Geodetic Survey ตามลำดับ ที่เป็นที่รู้จักกันมากที่สุดก็คือ Vincenty’s Formulae นี้เอง เขาเผยแพร่สูตรนี้เมื่อปี 1975 ซึ่งให้ความละเอียดสูงมากจนถึงครึ่งมิลลิเมตร
  • สำหรับสูตรของคุณ Thaddeus Vincenty ถ้าเห็นสูตรยาวๆแล้วสารอะดรีนาลีนหลั่งออกมาก็ ไปดาวน์โหลดได้ที่ http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
  • สำหรับสูตรที่จะนำเสนอในที่นี้ผมแปลมาจาก Javascript ดูได้ที่ http://www.movable-type.co.uk/scripts/latlong-vincenty.html สำหรับสูตรของ Vincenty จะมีการอินทิเกรต แต่ในทางโปรแกรมมิ่งจะใช้เทคนิค iteration แทนทำให้สูตรที่จะโปรแกรมดูกระชับมาก
  • สำหรับสูตรที่นำสรุปจาก Vincenty’s formulae คำนวณหาระยะทางและอะซิมัทจากจุด 2 จุด แสดงให้เห็นดังข้างล่าง
a, b = major & minor semiaxes of the ellipsoid
f = flattening (ab)/a
φ1, φ2 = geodetic latitude
L = difference in longitude
U1 = atan((1−f).tanφ1) (U is ‘reduced latitude’)
U2 = atan((1−f).tanφ2)
λ = L (first approximation)
iterate until change in λ is negligible (e.g. 10-12 ≈ 0.06mm) {
sinσ = √[ (cosU2.sinλ)² + (cosU1.sinU2 − sinU1.cosU2.cosλ)² ] (14)
cosσ = sinU1.sinU2 + cosU1.cosU2.cosλ (15)
σ = atan2(sinσ, cosσ) (16)
sinα = cosU1.cosU2.sinλ / sinσ (17)
cos²α = 1 − sin²α (trig identity; §6)
cos2σm = cosσ − 2.sinU1.sinU2/cos²α (18)
C = f/16.cos²α.[4+f.(4−3.cos²α)] (10)
λ′ = L + (1−C).f.sinα.{σ+C.sinσ.[cos2σm+C.cosσ.(−1+2.cos²2σm)]} (11)
}
u² = cos²α.(a²−b²)/b²
A = 1+u²/16384.{4096+u².[−768+u².(320−175.u²)]} (3)
B = u²/1024.{256+u².[−128+u².(74−47.u²)]} (4)
Δσ = B.sinσ.{cos2σm+B/4.[cosσ.(−1+2.cos²2σm) − B/6.cos2σm.(−3+4.sin²σ).(−3+4.cos²2σm)]} (6)
s = b.A.(σ−Δσ) (19)
α1 = atan2(cosU2.sinλ, cosU1.sinU2 − sinU1.cosU2.cosλ) (20)
α2 = atan2(cosU1.sinλ, −sinU1.cosU2 + cosU1.sinU2.cosλ) (21)

where :

    • s is the distance (in the same units as a & b)
    • α1 is the initial bearing, or forward azimuth
    • α2 is the final bearing (in direction p1→p2)

สิ่งที่ต้องเตรียมก่อนจะโปรแกรม

เริ่มต้นโปรแกรม

  • ผมจะสร้างยูนิตที่ทำหน้าที่คำนวณชื่อ GeodesicCompute (อยู่ในไฟล์ geodesiccompute.pas) โค๊ดดังที่กล่าวไปแล้วจาก Javascript แต่ดัดแปลงให้เป็นโค๊ด OOP ในยูนิตนี้มี 2 class คือ TGeodesicInverse ทำหน้าที่คำนวณระยะทางและอะซิมัทเมื่อป้อนค่าพิกัดเป็น Latitude,Longitude 2 จุด เรียกการคำนวณนี้ว่า Inverse
  • และอีกยูนิตหนึ่งคือ TGeodesicDirect(อยู่ในไฟล์ geodesiccompute.pas) ทำหน้าที่คำนวณหาค่าพิกัดจุดที่สอง (Lat/Long) เมื่อกำหนดค่าพิกัดจุดที่หนึ่งให้และ กำหนดอะซิมัทและระยะทาง เรียกการคำนวณนี้สั้นๆว่า Direct

GeodesicCompute Unit

  • โค๊ดที่ผมเขียนดัดแปลงมาจาก Javascript เป็นไปดัง source code ด้านล่าง
unit GeodesicCompute;

{$mode objfpc}{$H+}

interface

uses
 Classes, SysUtils, GeoEllipsoids, Math;

 type
 TZoneHemisphere = (hemEast, hemWest, hemNorth, hemSouth);

 TCoordinate = record
   X : double;
   Y : double;
 end;

//
//  Vincenty Inverse Solution of Geodesics on the Ellipsoid.
//  Calculate geodesic distance (in m) between two points
//  specified by latitude/longitude  (in numeric degrees) using
//   Vincenty inverse formula for ellipsoids.
//

TGeodesicInverse = class
 private
   FEllipsoid : TEllipsoid;
   FGeoCoor1  : TCoordinate;
   FGeoCoor2  : TCoordinate;
   FEllDist, FInitAzi, FFinalAzi : double;
   procedure SetGeoCoordinate1 (const geocoor1 : TCoordinate);
   procedure SetGeoCoordinate2 (const geocoor2 : TCoordinate);//input
   procedure SetEllipsoid(const Ellipsoid : TEllipsoid);//input
   function  GeodesicInverse(lat1, lon1, lat2, lon2 : double) : double;
 public
//input
   property GeoCoordinate1 : TCoordinate write SetGeoCoordinate1;
   property GeoCoordinate2 : TCoordinate write SetGeoCoordinate2;
//output
   property Distance : double read FEllDist;
   property InitialAzimuth  : double read FInitAzi;
   property FinalAzimuth : double read FFinalAzi;
//input
   property Ellipsoid : TEllipsoid write SetEllipsoid;
//main computation
   procedure Compute;

   constructor Create;
   destructor Destroy; override;
 end;

//  Calculate destination point given start point lat/long
//  (numeric degrees), azimuth or bearing (numeric degrees)
//  & distance (in m).
//
//  from: Vincenty direct formula - T Vincenty,
//       "Direct and Inverse Solutions of Geodesics on the
//        Ellipsoid with application of nested equations",
//        Survey Review, vol XXII no 176, 1975
//        http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf

TGeodesicDirect = class
 private
   FEllipsoid : TEllipsoid;
   FGeoCoor1  : TCoordinate;
   FGeoCoor2  : TCoordinate;
   FAzimuth   : double;
   FDistance  : double;
   procedure SetGeoCoordinate1 (const geocoor1 : TCoordinate);
   procedure SetEllipsoid(const Ellipsoid : TEllipsoid);
   procedure SetAzimuth(const azi : double);
   procedure SetDistance(const dist : double);
   function  GeodesicDirect(lat, lon, azimuth, distance : double) : TCoordinate;
 public
   //input
   property Ellipsoid : TEllipsoid read FEllipsoid write SetEllipsoid;
   property GeoCoordinate1 : TCoordinate write SetGeoCoordinate1;
   property Azimuth : double write SetAzimuth;
   property Distance : double write SetDistance;
  //Output
   property GeoCoordinate2 : TCoordinate read FGeoCoor2;
  //main computation
   procedure Compute;

   constructor Create;
   destructor Destroy; override;
 end;
implementation

// TGeodesicInverse
procedure TGeodesicInverse.SetGeoCoordinate1 (const geocoor1 : TCoordinate);
begin
  FGeoCoor1 := geocoor1;
end;

procedure TGeodesicInverse.SetGeoCoordinate2 (const geocoor2 : TCoordinate);
begin
  FGeoCoor2 := geocoor2;
end;

procedure TGeodesicInverse.SetEllipsoid(const Ellipsoid : TEllipsoid);
begin
  FEllipsoid.EllipsoidName := Ellipsoid.EllipsoidName;
  FEllipsoid.MajorAxis := Ellipsoid.MajorAxis;
  FEllipsoid.InvFlattening := Ellipsoid.InvFlattening;
end;

procedure TGeodesicInverse.Compute;
begin
  FEllDist := GeodesicInverse(FGeoCoor1.Y, FGeoCoor1.X, FGeoCoor2.Y, FGeoCoor2.X);
end;

function TGeodesicInverse.GeodesicInverse(lat1, lon1, lat2, lon2 : double) : double;
var
  a, b, f , e : double;
  dlat1, dlon1, dlat2, dlon2 : double;
  L, U1, U2 : double;
  sinU1, cosU1, sinU2, cosU2 : double;
  lambda, lambdaP, sinlambda, coslambda, sinsigma : double;
  cossigma, sinalpha, sigma, cosSqalpha, cos2sigmaM : double;
  cossigmaM, C, uSq, AA, BB, deltaSigma : double;
  az1, az2 : double;
  iterLimit : integer;
  diff : double;
begin
  dlat1 := PI/180 * lat1;
  dlon1 := PI/180 * lon1;
  dlat2 := PI/180 * lat2;
  dlon2 := PI/180 * lon2;

  a := FEllipsoid.MajorAxis;
  f := 1/FEllipsoid.InvFlattening;
  b := a*(1 - f);
  L := dlon2-dlon1;
  U1 := arctan((1-f) * tan(dlat1));
  U2 := arctan((1-f) * tan(dlat2));
  sinU1 := sin(U1);
  cosU1 := cos(U1);
  sinU2 := sin(U2);
  cosU2 := cos(U2);
  lambda := L;
  iterLimit := 100;

  repeat
    sinlambda := sin(lambda);
    coslambda := cos(lambda);
    sinsigma  := sqrt((cosU2*sinLambda) * (cosU2*sinlambda) +
                      (cosU1*sinU2-sinU1*cosU2*coslambda) *
                      (cosU1*sinU2-sinU1*cosU2*coslambda));
    if (sinsigma = 0) then
    begin
      result := 0; // co-incident points
      exit;
    end;
    cossigma := sinU1*sinU2 + cosU1*cosU2*coslambda;
    sigma := arctan2(sinsigma, cossigma);
    sinalpha := cosU1 * cosU2 * sinlambda / sinsigma;
    cosSqAlpha := 1 - sinAlpha*sinAlpha;
    cos2SigmaM := cosSigma - 2*sinU1*sinU2/cosSqAlpha;
    if (IsNaN(cos2sigmaM)) then cosSigmaM := 0; // equatorial line: cosSqAlpha=0 (6)
    C := f/16*cosSqalpha*(4+f*(4-3*cosSqalpha));
    lambdaP := lambda;
    lambda := L + (1-C) * f * sinAlpha *
              (sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*
              (-1+2*cos2SigmaM*cos2SigmaM)));
    dec(iterLimit);
    diff := abs(lambda-lambdaP);
  until ((diff < 1e-12) and (iterLimit > 0));

  if (iterLimit = 0) then
  begin
    result := NaN;   // formula failed to converge
    exit;
  end;

  uSq := cosSqAlpha * (a*a - b*b) / (b*b);
  AA := 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
  BB := uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));
  deltaSigma := BB*sinsigma*(cos2sigmaM+BB/4*
               (cossigma*(-1+2*cos2sigmaM*cos2sigmaM)-
                BB/6*cos2sigmaM*(-3+4*sinsigma*sinsigma)*
               (-3+4*cos2sigmaM*cos2sigmaM)));
  result := b*AA*(sigma-deltaSigma);  //return Ellipse Distance.
  //initial azimuth.
  az1 := arctan((cosU2*sinlambda)/(cosU1*sinU2 - sinU1*cosU2*coslambda));
  //final azimuth.
  az2 := arctan((cosU1*sinlambda)/(-sinU1*cosU2 + cosU1*sinU2*coslambda));

  FInitAzi := az1 * 180/PI;
  if (FInitAzi < 0) then FInitAzi := 360 + FInitAzi;

  FFinalAzi := az2 * 180/PI;
  if (FFinalAzi < 0) then FFinalAzi := 360 + FFinalAzi;
end;

constructor TGeodesicInverse.Create;
begin
  inherited create;
  FEllipsoid := TEllipsoid.Create;
end;

destructor TGeodesicInverse.Destroy;
begin
  FEllipsoid.Free;
  inherited Destroy;
end;

// TGeodesicDirect
procedure TGeodesicDirect.SetGeoCoordinate1 (const geocoor1 : TCoordinate);
begin
  FGeoCoor1 := geocoor1;
end;

procedure TGeodesicDirect.SetAzimuth (const azi : double);
begin
  FAzimuth := azi;
end;

procedure TGeodesicDirect.SetDistance (const dist : double);
begin
  FDistance := dist;
end;

procedure TGeodesicDirect.SetEllipsoid(const Ellipsoid : TEllipsoid);
begin
  FEllipsoid.EllipsoidName := Ellipsoid.EllipsoidName;
  FEllipsoid.MajorAxis := Ellipsoid.MajorAxis;
  FEllipsoid.InvFlattening := Ellipsoid.InvFlattening;
end;

procedure TGeodesicDirect.Compute;
begin
  FGeoCoor2 := GeodesicDirect(FGeoCoor1.Y, FGeoCoor1.X, FAzimuth, FDistance);
end;

function TGeodesicDirect.GeodesicDirect(lat, lon, azimuth, distance : double) : TCoordinate;
var
  a, b, f , e, s : double;
  dlat1, dlon1, dlat2, dlon2 : double;
  sinU1, alpha1, sinAlpha, sinAlpha1, cosAlpha1, tanU1, cosU1, sigma1 : double;
  cos2SigmaM, sinSigma, cosSigma, cosSqAlpha, uSq, deltaSigma : double;
  AA, BB, sigma, sigmaP, tmp, lambda, C, L, lat2, revAz : double;

begin
  dlat1 := PI/180 * lat;
  dlon1 := PI/180 * lon;

  a := FEllipsoid.MajorAxis;
  f := 1/FEllipsoid.InvFlattening;
  b := a*(1 - f);

  s := distance;
  alpha1 := azimuth * PI/180.0;
  sinAlpha1 := sin(alpha1);
  cosAlpha1 := cos(alpha1);

  tanU1 := (1-f) * tan(dlat1);
  cosU1 := 1 / sqrt((1 + tanU1*tanU1));
  sinU1 := tanU1*cosU1;
  sigma1 := arctan2(tanU1, cosAlpha1);
  sinAlpha := cosU1 * sinAlpha1;
  cosSqAlpha := 1 - sinAlpha*sinAlpha;
  uSq := cosSqAlpha * (a*a - b*b) / (b*b);
  AA := 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
  BB := uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));

  sigma := s / (b*AA);
  sigmaP := 2*PI;
  while (abs(sigma-sigmaP) > 1e-12) do
  begin
    cos2SigmaM := cos(2*sigma1 + sigma);
    sinSigma := sin(sigma);
    cosSigma := cos(sigma);
    deltaSigma := BB*sinSigma*(cos2SigmaM+BB/4*
                 (cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-
                  BB/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*
                 (-3+4*cos2SigmaM*cos2SigmaM)));
    sigmaP := sigma;
    sigma := s / (b*AA) + deltaSigma;
  end;

  tmp := sinU1*sinSigma - cosU1*cosSigma*cosAlpha1;
  dlat2 := arctan2(sinU1*cosSigma + cosU1*sinSigma*cosAlpha1,
          (1-f)*sqrt(sinAlpha*sinAlpha + tmp*tmp));
  lambda := arctan2(sinSigma*sinAlpha1,
  cosU1*cosSigma - sinU1*sinSigma*cosAlpha1);
  C := f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha));
  L := lambda - (1-C) * f * sinAlpha *
       (sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*
       (-1+2*cos2SigmaM*cos2SigmaM)));

  revAz := arctan2(sinAlpha, -tmp);  // final bearing
  result.X := lon + L * 180/PI;
  result.Y := dlat2 * 180/PI;
end;

constructor TGeodesicDirect.Create;
begin
  inherited create;
  FEllipsoid := TEllipsoid.Create;
end;

destructor TGeodesicDirect.Destroy;
begin
  FEllipsoid.Free;
  inherited Destroy;
end;

end.

  • สำหรับ source code ที่แสดงดังข้างต้น เป็น class ง่ายๆ เช่น TGeodesicInverse จะมี property สำหรับ input อยู่ 3 ตัวคือ ค่ารูปทรงสัณฐาน Ellipsoid, GeoCoordinate1 และ GeoCoordinate2 เมื่อ input ค่า lat,long 2 ค่าพิกัดให้เสร็จ แล้วก็เรียกเมธอด Compute เพื่อทำการคำนวณ จากนั้นเรียกค่าจาก 2 property คือ Azimuth และ Distance ซึ่งจะเป็นผลลัพธ์การคำนวณ
  • ส่วนคลาส TGeodesicDirect ก็เหมือนๆกัน มี property 4 ตัว คือ Ellipsoid, GeoCoordinate1, Azimuth และ Distance เมื่อรับค่ามาก็จะเรียกเมธอด Compute และให้ผลลัพธ์เป็นค่าพิกัด (lat,long) ของจุดที่ต้องการหา

ปัญหาความเข้ากันของ Widget ในแต่ละ Platform

  • รูปที่จะแสดงให้ดูเป็นปัญหาของ Widget ผมเพิ่งจะประสบความสำเร็จติดตั้ง Widget ของ Qt บน PCLinuxOS เมื่อรันโปรแกรมด้วยชุด source code เดียวกันเปรียบเทียบกับ Win32 บนวินโดส์ และ Gtk2 บน Ubuntu  พบว่า ขนาดของ Lable, EditBox, ComboBox, Groupbox ทั้งหลายบน Gtk2 มีขนาดใหญ่กว่าเพื่อน แต่พบว่าขนาดของ object พวกนี้บน Win32 และ Qt มีขนาดใกล้เคียงกัน เพื่อให้เข้ากับ Gtk2 จึงต้องมาขยับ เช่นตัว Label กับ EditBox จะต้องอยู่ห่างกันพอสมควร ไม่งั้นตัว Label จะล้ำเข้าไปหา EditBox แต่ความสวยงามดูๆ QT จะสวยเนียนกว่าเพื่อน รองลงมาเป็น Gtk2 กับ Win32 ดูรูปผลลัพธ์ของโปรแกรมที่ผมเขียนในตอนนี้
Geodesic computation on Qt widget.

โปรแกรมคำนวณหาระยะทางและอะซิมัทบน Ellipsoid รันด้วย Qt (PCLinuxOS)

โปรแกรมคำนวณหาระยะทางและอะซิมัทบน Ellipsoid รันด้วย Gtk2 (Ubuntu)

โปรแกรมคำนวณหาระยะทางและอะซิมัทบน Ellipsoid รันด้วย Gtk2 (Ubuntu)

โปรแกรมคำนวณหาระยะทางและอะซิมัทบน Ellipsoid รันด้วย Windows

โปรแกรมคำนวณหาระยะทางและอะซิมัทบน Ellipsoid รันด้วย Windows

  • ทิ้งท้ายกันก่อนในตอนนี้ ติดตามกันตอนหน้า
Categories: Linux, Programming, Windows ป้ายกำกับ:, , , , ,

การสร้าง (Create/Generate) DEM จาก DTM (3D Points) ด้วย Global Mapper

  • ตัว Global Mapper สำหรับความรู้สึกผมนั้นดูเรียบง่ายและทรงพลัง การสร้าง DEM (Digital Elevation Model) จากโปรแกรมดังๆอย่าง Erdas Imagine หรือ ArcGIS ยังดูยากไปนิด มาดูกันว่า Global Mapper นั้นเรียบง่ายเพียงใด

Global Mapper V.10

  • ผมมีเหตูต้องศึกษาเรื่องน้ำที่แม่น้ำตะนาวศรี(Tanintharyi River) ชื่อน่าจะคุ้นคนไทยเรา หลักสูตรวิชาภูมิศาสตร์จะมีเทือกเขาตะนาวศรีคั่นไทยกับพม่าอยู่ และก็แม่น้ำที่ไหลเคียงข้างกับเทือกเขาตะนาวศรีด้านพม่าก็คือแม่น้ำตะนาวศรีนี้เอง คนพม่าเรียกแม่น้ำนี้ว่า ตาเน็นตายี่
  • ผมใช้ข้อมูลความสูง ซึ่งมีอยู่ 2 sources ต้องนำมาผสมกัน ข้อมูลแรกเป็น DTM (Digital Terrain Model) ของพม่า โดยที่ DTM ได้ข้อมูลมาเป็นจุดซึ่งวัดมาจากภาพถ่าย Stereo (พม่าจัดทำแผนที่ สเกล 1:50000 ใหม่ทั้งประเทศด้วยกระบวนการภาพถ่ายทางอากาศเมื่อปี 2005 เป็นแผนที่ภาพถ่ายทางอากาศ (Orthophoto map) และแผนที่ลายเส้น) ข้อมูลที่สองเป็น DEM ของ SRTM DEM ที่ผมเคยเขียนวิธีการ download ไปแล้ว
3D points as DTM (Digital Terrain Model)

3D points as DTM (Digital Terrain Model)

  • จากรูปด้านบนจะเห็นจุดที่เรียงเป็นกริดนั้น import มาจาก SRTM DEM ส่วนจุดที่กระจายไม่เป็นระเบียบเป็น 3D points ที่ได้จากก  DTM  ของพม่า
  • ที่เมนูหลักคลิกที่ Tools > Control Center… หรือใช้ shortcut ก็ได้ด้วยการกด Alt + C ผมมี layer ที่เก็บ points รวมไว้แล้วเป็น layer เดียว ที่เลเยอร์นี้คลิกขวาเลือก Create Elevation Grid form 3D Vector Data…
Create elevation grid.

Create elevation grid.

  • จะมีไดอะล็อก ให้ตั้งเงื่อนไขของ Elevation Grid ก่อน
ตั้งค่าเพื่อสร้าง Elevation grid

ตั้งค่าเพื่อสร้าง Elevation grid

  • จากรูปด้านบนตั้งชื่อที่ Description ผมตั้งเป็น Tanintaryi_DEM ระยะ spacing ของกริดผมเลือกเป็น 10 เมตรทั้ง x และ y จากนั้นตัวอื่นทิ้งไว้ตามค่าปริยายคลิกที่ OK เพื่อเริ่มสร้าง elevation grid โปรแกรม Global Mapper จะเริ่มสร้าง TIN (Triangulation Irregular Network) ใช้เวลาประมาณ 3-4 นาทีต่อจำนวนจุดประมาณ 1.4 ล้านจุดก็ถือว่าเร็ว ต่อจากนั้นก็ฟอร์มเป็น Elevation Grid ซึ่งพร้อมจะสร้างเป็น DEM ต่อไป
Generated DEM

Generated DEM

  • จากรูปด้านบนเมื่อกด ALT+C ตัว Control Center จะพ๊อพอัพออกมา ติ๊กที่ layer 3d points ออก ให้เห็น layer ใหม่คือ Tanintaryi_Dem อย่างเดียวจะเห็นรูป DEM สวยงามดังรูปข้างบน

การ Export เป็น DEM

  • ที่เมนูหลัก คลิกที่ File > Export Raster and Elevation Data > Export GeoTiff… เราจะ Export เป็น GeoTiff DEM
Export to DEM

Export to DEM

  • จะเห็นไดอะล็อกถาม options ของฟอร์แม็ต GeoTiff ตั้งค่าดังรูปด้านล่าง
ตั้งค่า GeoTiff

ตั้งค่า GeoTiff

  • จากรูปด้านบนเมื่อคลิก OK โปรแกรมจะให้ป้อนชื่อไฟล์ GeoTiff เมื่อป้อนชื่อไฟล์ รอประมาณ 3 นาทีก็จะเสร็จ ทดสอบดูว่ามีอะไรผิดพลาดหรือไม่ลองไปเปิดด้วยโปรแกรมอื่นผมใช้ Viewer ของ Erdas Imagine ก็ OK
Viewer by Erdas Imagine

Viewer by Erdas Imagine

Categories: GIS, Windows ป้ายกำกับ:, , , , , ,

การติดตั้ง Lazarus แบบ Subversion บน Windows

  • Blog ตอนก่อนผมเขียนเรื่องติดตั้ง lazarus แบบ subversion บน linux ละเอียดพอสมควร ติดค้างเรื่องติดตั้ง lazarus บนวินโดส์ และก็เหมือนเดิมผมยังแนะนำให้ใช้ subversion เหมือนบนลินุกซ์ จะได้ version ที่ update ตลอด ไม่ต้องห่วงเรื่อง bug แก้ไขได้เร็วพอสมควร
  • เป้าหมายการติดตั้งของ lazarus แบบ subversion ของเราจะอยู่ที่ fpc 2.5.1 และ lazarus อยู่ที่เวอร์ชั่น 0.9.29 beta

Lazarus บน Windows ติดตั้งด้วย subversion

Lazarus บน Windows ติดตั้งด้วย subversion

สิ่งที่ต้องเตรียมก่อนจะติดตั้ง Lazarus

  1. Download โปรแกรม Subversion สำหรับ Window ชื่อ TortoiseSVN ที่ http://tortoisesvn.net/downloads มีทั้งเวอร์ชั้น 32 bit และ 64  bit
  2. Download Free pascal(fpc) ที่ build เป็นไบนารีพร้อมจะติดตั้ง มาเรียบร้อยแล้วที่ ftp://ftp.freepascal.org/pub/fpc/dist/i386-win32-2.2.2/fpc-2.2.2.i386-win32.exe ก็เหมือนใน linux เวลาเราดาวน์โหลด freepascal และ lazarus แบบ subversion จะได้ source code มาทั้งดุ้น จะคอมไพล์ก็ทำไม่ได้เนื่องจากไม่มีไฟล์ไบนารีของ fpc.exe, ppc386.exe, … มาคอมไพล์ จึงต้องดาวน์โหลดตัวเก่าคือ fpc v.2.2.2 มาติดตั้งก่อน เมื่อคอมไพล์ source code ของ fpc แบบ subversion ผ่าน เราจะได้ไฟล์ไบนารี  fpc เวอร์ชั่น 2.5.1 (ในขณะที่เขียนอยู่) ที่พร้อมจะนำไปคอมไพล์ source code แทนที่ fpc เวอร์ชั่นเดิม
  3. ทำการติดตั้ง TortoiseSVN และ fpc
  4. ร้างโฟลเดอร์สำหรับ source code ของ fpc และ lazarus ผมเองสร้างโฟลเดอร์ fpcsvn และ lazarussvn อยู่ใต้โฟลเดอร์ D:/Sourcecodes/Lazarus/ ก็คือ D:/Sourcecodes/Lazarus/fpcsvn และ D:/Sourcecodes/Lazarus/lazarussvn

การใช้ TortoiseSVN เพื่อดาวน์โหลด fpc และ lazarus

  • ใช้ Explorer เข้าไปที่โฟลเดอร์ที่เราสร้างไว้ เช่นที่ผมกล่าวไว้ข้างต้น คลิกขวาที่โฟลเดอร์ fpcsvn ก่อน แล้วคลิกที่ SVN Checkout… ดังรูปด้านล่าง

  • ตรง URL of repository  ป้อน http://svn.freepascal.org/svn/fpc/trunk fpc ตรง Checkout directory ป้อนโฟลเดอร์ที่จะใช้เก็บ fpc svn จากนั้นคลิก OK

เรียกโปรแกรม TortoiseSVN

เรียกโปรแกรม TortoiseSVN

  • โปรแกรม TortoiseSVN จะเริ่มดาวน์โหลด fpc svn ดังรูปด้านล่าง
้ป้อน URL ของ fpc svn

ป้อน URL ของ fpc svn

  • รอนานพอสมควรประมาณ 37 นาทีเมื่อได้ source code ของ fpc svn ครบ จากรูปด้านล่าง จะเห็นหมายเลข revision คือ 13842 (r13842) เราก็คลิกที่ OK เพื่อออกมา
ดาวน์โหลด fpc svn เสร็จสิ้น

ดาวน์โหลด fpc svn เสร็จสิ้น

  • เราจะทิ้งโฟลเดอร์ fpcsvn ไว้ก่อนเราจะมาต่อที่ lazarus เช่นเดียวกันที่นี้คลิกขวาที่โฟลเดอร์ของ lazarussvn แล้วเลือกจากป๊อปอัพเมนู “SVN Checkout…”  ป้อน URL of repository เป็น
ใส่รูปภาพ้ป้อน URL และ โฟลเดอร์ปลายทางของ lazarus svn

้ป้อน URL และ โฟลเดอร์ปลายทางของ lazarus svn

  • คลิกที่ OK โปรแกรม TortoiseSVN จะเริ่มตรวจสอบและดาวน์โหลด lazarus svn ให้เราดังรูปด้านล่าง
Checkout and download lazarus svn.

Checkout and download lazarus svn.

  • เมื่อ TortoiseSVN ดาวน์โหลด source code ไฟล์ subversion มาให้ครบทั้ง fpc และ lazarus แล้ว เราจะเริ่มคอมไพล์ จะเริ่มจาก fpc ก่อน

คอมไพล์ fpc ก่อน

  • เราจะมาใช้ไฟล์ไบนารีของ fpc version 2.2.2 ที่เราดาวน์โหลดมาก่อนหน้าแล้วติดตั้งไปแล้ว เพื่อความแน่ใจเราใช้ Windows Explorer เปิดดูโฟลเดอร์ C:\fpc\2.2.2\bin\i386-win32 และตรวจดูไฟล์มีไฟล์ fpc.cfg หรือไม่ เป็นไฟล์ที่สำคัญมากสำหรับการคอมไพล์ lazarus ให้สำเร็จในเวลาต่อไป ตรวจดูจะเห็นไบนารีเช่น fp.exe, fpc.exe, ppc386.exe,…
โฟลเดอร์ fpc เวอร์ชั่น 2.2.2

โฟลเดอร์ fpc เวอร์ชั่น 2.2.2

  • การคอมไพล์จะต้องทำใน command line ซึ่งไม่ยากอยู่แล้วสำหรับโปรแกรมเมอร์ทั้งสมัครเล่นและอาชีพ ผมจะลดวิธีการอ้าง path ในวินโดส์ ซึ่งไม่มีอะไรวิธีการก็ง่ายๆ เริ่มด้วยการคลิกที่เมนู Start > Run… > พิมพ์ cmd กด Enter พิมพ์คำสั่งเพื่อเปลี่ยนโฟลเดอร์ดังรูปด้านล่าง

เปลี่ยนโฟลเดอร์ไปที่ fpc v.2.2.2

เปลี่ยนโฟลเดอร์ไปที่ fpc v.2.2.2

  • ที่ command line เราจะสลับ drive ไปที่โฟลเดอร์ที่เก็บ fpc svn โดยที่อย่าเปลี่ยน path ของ drive c: ที่เก็บ fpc อย่างเด็ดขาด พิมพ์คำสั่งเปลี่ยนไป drive d: และเริ่มคอมไพล์ ดังรูปด้านล่าง

เปลี่ยน drive ทำการคอมไพล์ fpc svn

เปลี่ยน drive ทำการคอมไพล์ fpc svn

  • ไม่มีปัญหาอะไร การคอมไพล์จะเริ่มขึ้น ดังรูปข้างล่าง การคอมไพล์ fpc บนวินโดส์จะช้ากว่าบน linux เล็กน้อย
คอมไพล์ fpc svn

ภาพตัวอย่างตอนคอมไพล์ fpc svn


  • เมื่อคอมไพล์ fpc svn เสร็จแล้วด้วยคำสั่ง make clean all ที่กล่าวไปข้างต้น เราจะ install ด้วยคำสั่ง c:make install (ไม่ลืมว่าตอนนี้อิงคำสั่งจาก drive c: ที่เราเปิดพาธของ c:\fpc\2.2.2\bin\i386-win32 อยู่)
make install

make install

  • เมื่อ make install เราจะได้โฟลเดอร์ที่ถูกสร้างขณะ install มาให้อัตโนมัติคือ โฟลเดอร์ pp จะอยู่ที่ drive d: => d:\pp มีโฟลเดอร์ย่อยอยู่หลายโฟลเดอร์ รวมทั้งไฟล์ไบนารีคอมไพเลอร์ (โฟลเดอร์ pp จะขึ้นอยู่กับ drive ของ fpc svn ถ้าอยู่ drive e: เมื่อ compile & install จะสร้างโฟลเดอร์ที่ drive e: ให้เป็น e:\pp)
โฟลเดอร์ที่ได้จากการ make install

โฟลเดอร์ที่ได้จากการ make install

ไม่ต้อง Uninstall fpc เวอร์ชั่นเก่า

  • ในวินโดส์จะต่างจาก linux อยู่คือไฟล์ไบนารี make จะมาพร้อมกับ OS เลย แต่ในวินโดส์ถ้าผู้ใช้ไม่ได้ติดตั้งคอมไพเลอร์ตัวใดตัวหนึ่งเช่น Visual Studio, Delphi จะไม่มีไฟล์ make.exe ให้ใช้งาน ถ้าเรา uninstall fpc ตัวเก่าออก ไฟล์ make.exe จะหายไปด้วย เราจะไม่ uninstall แต่จะใช้วิธี copy วิธีปฏิบัติเปิด windows explorer ที่โฟลเดอร์ d:\pp copy โฟลเดอร์ภายใต้โฟลเดอร์ pp (รูปข้างบนมีโฟลเดอร์ bin, doc, examples, msg, units) ไป paste ทับโฟลเดอร์ c:\fpc\2.2.2\ โดยใช้ windows explorer
  • เมื่อ paste ทับโฟลเดอร์ c:\fpc\2.2.2\  ตอบ Overwrite all ที่นี้ภายใต้โฟลเดอร์ fpc จะมี binary ไฟล์ที่เป็น fpc.exe, fp.exe, ppc386.exe ซึ่งเป็นเวอร์ชั่นใหม่พร้อมใช้งาน) ลอง double click ppc386.exe จะแสดงเวอร์ชั่นบน dos box ดังรูปด้านล่าง
ตรวจสอบเวอร์ชั่นของ ppc386.exe

ตรวจสอบเวอร์ชั่นของ ppc386.exe

  • จากรูปด้านบนจะเห็นเป็นเวอร์ชั่น fpc version 2.5.1 พร้อมไปต่อได้ ขั้นต่อไปเราจะหักเอาด้วยกำลังโดยการเปลี่ยนชื่อโฟลเดอร์จาก c:\fpc\2.2.2 เปลี่ยนเป็น c:\fpc\2.5.1 เมื่อเปลี่ยนชื่อโฟลเดอร์แล้ว ขั้นต่อไปจะเข้าไปแก้ไขไฟล์ fpc.cfg แล้วค่อยทำการ compile lazarus svn (อย่าลืมว่าอนาคตถ้า fpc เปลีียนเวอร์ชั่น เช่นเปลีียนไปเป็น 2.7.1 ก็ให้มาเปลี่ยนชื่อโฟลเดอร์จาก 2.5.1 เป็น 2.7.1)

แก้ไขไฟล์ fpc.cfg เพื่อคอมไพล์ lazarus

  • ในการคอมไพล์ lazarus svn ต้องการไฟล์ไบนารีและ units ของ fpc กุญแจสำคัญคือไฟล์ fpc.cfg ที่เก็บอยู่ในโฟลเดอร์ c:\fpc\2.5.1\bin\i386-win32
ตรวจสและแก้ไขไฟล์ fpc.cfg
  • จากรูปด้านบนผมเปิด fpc.cfg ด้วย Notepad++ ซึ่งเป็นฟรีแวร์ที่ใช้ดีมาก ผมแนะนำให้แก้ไขไฟล์ ด้วยการเปลี่ยนคำว่า 2.2.2 เป็น 2.5.1 ให้หมด (ใช้เมนู replace… )
แก้ไขไฟล์ fpc.cfg ให้อ้างอิง fpc เวอร์ชั่นใหม่

แก้ไขไฟล์ fpc.cfg ให้อ้างอิง fpc เวอร์ชั่นใหม่

  • ในการคอมไพล์ lazarus จะมองหา unit ของ fpc โดยที่ lazarus จะอ่านจาก fpc.cfg ซึ่งจะเห็นโฟลเดอร์ที่อ้างถึง units ที่อยู่ในโฟลเดอร์ c:\FPC\2.5.1\units\i386-win32  ($FPCTARGET จะแทน i386-win32 ในวินโดส์)

คอมไพล์ lazarus svn

  • ใช้วิธีการเหมือนคอมไพล์ fpc เปิด dos box ด้วย Start > Run… > cmd พิมพ์คำสั่งดังนี้
คอมไพล์ lazarus svn

คอมไพล์ lazarus svn

  • ถ้าไฟล์ fpc.cfg แก้ไขได้ถูกต้อง lazarus หาไฟล์ไบนารีและ units ได้ กระบวนการคอมไพล์จะเริ่มจนจบ จากนั้นต่อด้วยคำสั่ง c:make install

คำสั่ง make install สำหรับ lazarus svn

คำสั่ง make install สำหรับ lazarus svn

  • คำสั่ง make install จะ copy file ที่ได้จากการ compile & build จาก lazarus svn (d:\sourcecodes\lazarus\lazarussvn) เข้าไปที่ c:\lazarus\ ใช้เวลานานพอสมควรเนื่องจากมีไฟล์เป็นจำนวนมาก

Start lazarus

  • หลังจากที่รอกันมานานถ้ามาถึงขั้นนี้ความสำเร็จประมาณ 90% แล้ว เปิด windows explorer ไปที่  c:\lazarus\ มองหาไฟล์ startlazarus.exe แล้ว double click เพื่อรันได้เลย คลิกที่เมนู Help > About เพื่อดูเลข revision
lazarus แรกบนวินโดส์

lazarus แรกบนวินโดส์

การตั้งค่าให้ lazarus

  • เราจะมาตั้งค่าให้กับ lazarus ก่อนจะติดตั้ง component เพิ่มเติม ที่เมนเมนูคลิกที่ Environment > Options… ที่พาเนลด้านซ้ายคลิกที่ Environment ดูเรื่อง path ให้ดี จะมีที่สำคัญคือที่ผมเน้นไว้ ตรง Compiler path ก็เลือกพาธของ fpc ที่ FPC source directory ก็เลือกที่ที่เราติดตั้ง lazarussvn ก็คือ d:\sourcecodes\lazarus\lazarussvn แล้วก็ Make path ดูภาพด้านล่างประกอบ
ตั้งค่า path ให้ lazarus

ตั้งค่า path ให้ lazarus

  • จากไดอะล็อกรูปข้างบน ที่พาเนลด้านซ้ายคลิกไปที่ Debugger ตั้งค่าดังรูปด้านล่าง
ตั้งค่า Debugger

ตั้งค่า Debugger

Rebuild lazarus

  • เราจะทำการ rebuild เพื่อทดสอบว่าเราตั้งค่าต่างๆที่กล่าวมาแล้วว่าถูกต้องแล้วหรือยัง ถ้าถูกต้องการ rebuild ก็จะผ่าน ที่เมนเมนูคลิกที่ Tools > Configure “Build Lazarus…”

Rebuild lazarus

Rebuild lazarus

  • จากรูปด้านบนผมเลือกจะไม่ติดตั้ง Examples เพราะโอกาส fail สูงมาก เมื่อเลือกได้ดังรูปก็คลิกที่ Build ต่อ จากนั้นรออีกไม่นานโปรแกรม Lazarus ก็จะเริ่มมาใหม่

การติดตั้ง Component เพิ่มเติม

  • การติดตั้ง component เพิ่มเติมจากโฟลเดอร์ c:\lazarus\components ผมจะทดสอบโดยการติดตั้ง component “tdbf” ที่เมนเมนูคลิกที่ Package > Open package file (.lpk)… ไปที่ c:\lazarus\components\tdbf เลือกไฟล์ dbflaz.lpk คลิก Open จะเห็นไดอะล็อกดังรูปข้างล่าง
การติดตั้ง component เพิ่มเติมให้ lazarus

การติดตั้ง component เพิ่มเติมให้ lazarus

  • อย่าเพิ่งรีบคลิกที่ Install ให้คลิกที่ Compiler Options ก่อนจะพบไดอะล็อก ดังรูปข้างล่าง ตั้งค่า
การตั้งค่า OS Platform และ CPU Platform

การตั้งค่า OS Platform และ CPU Platform

  • จากรูปข้างบนที่พาเนลด้านซ้ายคลิกที่ Code จากนั้นเลือก Target OS (-T) เป็น Win32 เลือก Target CPU family (-P) เป็น i386 ถ้าไม่ได้ตั้งค่าการ Install component จะ fail ทันที่เสร็จแล้วคลิกที่ OK และจากไดอะล็อกก่อนหน้านี้ให้คลิกที่ Install เพื่อ lazarus จะได้ compile & rebuild จะเห็น component “tdbf” โผล่มาที่ palette ของ Data Access จากรูปจะเห็น component ที่ผมติดตั้งเพิ่ม SQLdb และคอมโพเน็นท์ของผมเอง OpenGPSX

Additional components

Additional components

  • ถ้าการติดตั้ง component เพิ่มเติมไม่มีปัญหา ก็สามารถนำ lazarus ไปพัฒนาโปรแกรม cross-platform ได้เลย

การ Update fpc และ lazarus ด้วย TortoiseSVN

  • การใช้ lazarus แบบ subversion สามารถจะ update ได้ทุกๆวันด้วย TortoiseSVN แต่ไม่จำเป็นเพราะการคอมไพล์จะใช้เวลามาก ผมเองจะ update ประมาณอาทิตย์ละครั้ง การ update ทำได้โดยเปิด windows explorer ที่โฟลเดอร์เราเก็บ svn คือ d:\sourcecodes\lazarus\ คลิกขวาที่โฟลเดอร์ fpcsvn และ lazarussvn เลือกใช้เมนู SVN Update… เมื่อได้ source code ที่ update แล้วก็เริ่มกระบวนการคอมไพล์ ที่ผมกล่าวมาแล้วข้างต้นจาก make clean all => make install

ปัญหาของ make install เมื่อ write permission denied

  • เมื่อทำการ build ด้วยคำสั่ง make install จะพบปัญหาว่าไม่สามารถเขียนไฟล์ได้เนื่องจาก permission denied ปัญหานี้จะพบในวินโดส์ไม่พบใน linux เพราะสิทธิ์ของผู้ดูแลระบบใน linux จะชัดเจน
  • วิธีการแก้ปัญหาถ้า make install ของ fpc svn ไม่ผ่านเนื่อง write permission denied ดูว่า source code ของ fpc svn อยู่ไหนเช่นอยู่ drive d: ให้ใช้ windows explorer ลบโฟลเดอร์ d:\pp ทิ้ง แล้วรัน make install ใหม่
  • วิธีการแก้ไขถ้า make install ของ lazarus svn ไม่ผ่านเนื่องจาก write permission denied ใ้ห้ลบโฟลเดอร์ c:\lazarus ทิ้งไปเลย ทำการรัน make install ใหม่
  • ผมนึกว่าการติดตั้งแบบ Subversion บน Windows น่าจะสั้นกว่าใน Linux แต่มีอะไรให้เขียนมากพอสมควรก็เลยยาวเหยียดกว่า แต่อย่างไรก็ตาม ถ้ามีโปรแกรมเมอร์ใช้ lazarus พัฒนาโปรแกรมกันมากๆ ชุมชน lazarus ก็จะใหญ่ขึ้น สุดท้ายจะมีคนเข้าไปช่วยพัฒนากันมากขึ้น lazarus ก็จะแข็งแกร่งขึ้นด้วยร่วมด้วยช่วยกันนี่แหละครับ และก็น่าภูมิใจเพราะเป็นของฟรีตามหลัก Opensource

การคำนวณการรังวัด GPS ด้วยวิธี Online service (ตอนที่ 2)

  • ตอนที่ 1 ผมแนะนำการใช้บริการการคำนวณรังวัด GPS บนอินเทอร์เน็ต(GPS Online Service)  ของ Natural Resources Canada (NRCan) ต่อไปผมจะใช้บริการ website ของ Australia หน่วยงานที่ให้บริการเรียกว่า AUSLIG Online GPS Processing Service (AUSPOS) โดยที่ AUSLIG ย่อมาจาก Geoscience Australia, The Australian Government mapping agency
  • ผมจะใช้ข้อมูลเดียวกัน (Rinex file) ส่งไปให้ website ของ AusPos คำนวณแล้วมาเปรียบเทียบค่ากับของ Canada (NRCan) ว่าค่าที่ได้มีความแตกต่างกันเพียงใด

AUSPOS – Online GPS Processing Service

  • ข้อมูลรังวัด GPS ที่ทาง AUSPOS อ้างว่าถ้ารังวัดประมาณ 24 ชั่วโมง ค่าพิกัดทางราบจะมีความถูกต้องดีกว่า 1 เซนติเมตร(น้อยกว่า 1 ซม.)  ค่าระดับจะถูกต้องประมาณ 1-2 เซนติเมตร แต่มีข้อแม้ว่าความถูกต้องขนาดนี้ หมุด GPS ที่เราทำการรังวัดต้องห่างหมุดโครงข่าย GPS ของ THE INTERNATIONAL GPS SERVICE (IGS) ประมาณ 1000 กม.
  • เดิมทีหมุดโครงข่าย GPS ของ AUSPOS มีเพียง 15 หมุดครอบคลุมทวีปออสเตรเลียและแอนตาร์กติก ซึ่งจำกัดการให้บริการสำหรับ Surveyor ของออสเตรเลียเท่านั้น ตอนหลัง AUSPOS ได้เข้าร่วมกับหมุดโครงข่ายของ IGS ทำให้รวมๆกันมีหมุด 200 หมุดทั่วโลก ทำให้ไม่ว่าอยู่ภูมิภาคใดของโลกก็สามารถเข้ามาใช้บริการได้ แต่สังเกตว่าย่านบ้านเราคือเอเชียตะวันออกเฉียงใต้มีหมุดของ IGS น้อยมาก
IGS Network

ภาพแสดงหมุดโครงข่าย IGS Network

  • เมื่อเราส่งข้อมูลการรังวัดเป็นไฟล์ Rinex ไปให้ AUSPOS จะใช้โปรแกรม MicroCosm (ชอง Van Martin Systems, Inc.) ในการคำนวณโดยใช้ข้อมูลวงโคจร GPS ของ IGS เช่น Precise satellite orbit, Earth Orientation parameters และก็ใช้ค่าพิกัดหมุดโครงข่าย GPS ของ IGS ที่อยู่ใกล้หมุด GPS ของเราจำนวน 3 หมุด
  • ข้อสังเกตถ้าเป็นการคำนวณของ NRCan ที่ผมกล่าวไปในตอนที่ 1 จะไม่ใช้ค่าพิกัดชองหมุด GPS โครงข่ายใดๆเลยมาช่วยคำนวณ ซึ่งน่าทึ่งมาก และใช้ Algorithms “Precise Point Positioning” อย่างเดียว เพียงแต่ข้อมูลที่นำมาใช้เหมือนกันคือ ข้อมูลวงโคจร GPS ของ IGS ทาง NRCan เรียกการคำนวณนี้ว่า Un-differenced carrier phase observations ส่วนของ AUSPOS เรียกว่า Double-difference carrier phase observations

สิ่งที่ต้องเตรียม

  • เตรียมไฟล์ข้อมูลรังวัด GPS ในรุปแบบ Rinex ในตอนที่ 1 ผมกล่าวถึงไฟล์ไบนารีของ GPS ที่ผมใช้รังวัดคือ Trimble 5700 ซึ่งก่อนการ log ข้อมูลการรังวัดจะป้อนค่าความสูงของจานรับสัญญาณ (Antenna) และชนิดของจานรับสัญญาณ  เมื่อแปลงเป็น .Dat และก็แปลงเป็น Rinex ด้วยโปรแกรม teqc ความสูงของ Antenna จะติดมาด้วย เมื่อส่งข้อมูลเข้า website NRCan ความสูงจะถูกตรวจสอบได้อัตโนมัติ แต่ถ้าตอนรังวัดป้อนตัวเลขผิด ค่าระดับที่ได้จากการคำนวณจาก NRCan ก็จะผิดตามต้องทำการทอนความสูงใหม่ให้ถูกต้อง
  • แต่ของ AUSPOS จะให้ป้อนความสูงใหม่เมื่อตอนส่งข้อมูล Rinex ไปให้พร้อมระบุชนิดของ Antenna ซึ่งข้อดีถ้าป้อนความสูงขณะรังวัดผิดก็แก้ใหม่ได้เลย แต่ข้อเสียก็ไม่มีอะไรมากก็แค่เสียเวลาป้อนเพิ่มเท่านั้น
  • ถ้าสงสัยว่า Antenna ของเรารุ่นไหนเข้าไปตรวจสอบได้ที่ http://www.ngs.noaa.gov/ANTCAL/ ตัวอย่างเช่นผมใช้ของ Trimble เมื่อเลือกรุ่นที่ใช้จะเห็นดังรูปข้างล่าง
จานรับสํญญาณ Trimble Zephyr Geodetic

จานรับสํญญาณ Trimble Zephyr Geodetic

การใช้บริการ GPS Online Service ของ AUSPOS

เลือกไฟล์ Rinex, ป้อนรุ่น Antenna , ความสูง และ email address

เลือกไฟล์ Rinex, ป้อนรุ่น Antenna , ความสูง และ email address

  • คลิกที่ submit จากนั้น website จะแจ้งให้ทราบว่า upload ไฟล์ไปสำเร็จหรือไม่

AUSPOS แจ้งสถานะการ upload

AUSPOS แจ้งสถานะการ upload

ภาพแสดงการ Upload ไฟล์ได้สำเร็จ
  • จากรูปด้านบนรอประมาณ 24 นาที website ก็ส่งผลการคำนวณมาทาง Email address ของเรา เป็นไฟล์ pdf ดาวน์โหลดมาดูตรวจสอบได้
ภาพแสดงหมุด GPS ที่ทำการรังวัด

ภาพแสดงหมุด GPS ที่ทำการรังวัด

  • จากรูปด้านบน จะเห็นหมุดโครงข่าย GPS ของ IGS ที่ AUSPOS เลือกมาเป็น Base Station หมุด kunm อยู่ที่คุณหมิง จีน ระยะทางประมาณ 1300 กม. หมุด iisc อินเดีย ระยะทางประมาณ 2200 กม. และหมุด ntus ที่สิงคโปร์ ระยะทาง 1560 กม.
ผลลัพธ์ที่เป็นค่าพิกัดหมุด 2462

ผลลัพธ์ที่เป็นค่าพิกัดหมุด 2762

  • ความจริง website ที่ให้บริการ GPS Online Service ไม่ได้มีเพียง 2 site ตามที่ผมแนะนำไปในตอนที่ 1 เช่นของ NOAA (National Oceanic And Atmospheric Administration) ที่  http://www.ngs.noaa.gov/OPUS/ ค่าพิกัดที่คำนวณได้จาก website นี้ค่อนข้างกระโดดจากกลุ่มไปมาก เลยยังไม่ขอแนะนำ
  • ที่น่าสนใจมากๆก็คือของ NASA โปรแกรมที่ให้บริการดั้งเดิมเป็นที่รู้จักกันมากคือ Auto-GIPSY ตอนนีุ้ถูกแทนที่ด้วย APPS (Automatic Precise Position Service) ดูตามชื่อแล้วโปรแกรมน่าจะถูกปรับปรุงให้ใช้ Algorithms “Precise Point Positioning”  ถ้าจะใช้บริการไปที่ http://apps.gdgps.net/ ผมจะไม่กล่าวถึงในรายละเอียดเพราะแต่ละ website  ที่ให้บริการแบบนี้จะคล้ายๆกัน
  • ผมขอเสนอผลลัพธ์การคำนวณที่ได้จาก NRCan จากตอนที่ 1, AUSPOS ที่กล่าวไปข้างต้น, จากของ NOAA และของ NASA (JPL) ค่าพิกัดทางราบและทางดิ่งของผลการคำนวณจะใกล้เคียงกันมาก ยกเว้นของ NOAA จะแสดงให้เห็นดังตารางข้างล่าง


Website ผู้ให้บริการ
Northing (N)
Easting (E)
Ortho Height (m.)
Natural Resorces Canada (NRCan) 1579064.181 398040.484 5.475
AUSLIG Online GPS Processing Service (AUSPOS) 1579064.178 398040.484 5.530
National Oceanic And Atmospheric Administration (NOAA)
1579062.003 398839.367 14.072
Jet Propulsion Laboratory (JPL / NASA) 1579064.211 398840.497 5.819
  • จากตารางผมสรุปได้ว่า ค่าพิกัดทางราบและทางดิ่ง ทั้ง 3 website (ยกเว้น NOAA) ให้ค่าพิกัดที่น่าพอใจมากเกาะกลุ่มกันทั้งค่าพิกัดทางราบ (N,E) และค่าพิกัดทางดิ่ง (ค่าระดับ ที่เป็น Orthometric Height (MSL)) จะบอกว่าใครถูกต้องที่สุดคงไม่ได้ ขึ้นอยู่กับความสะดวกและความชอบมากกว่า
Categories: GPS, Surveying ป้ายกำกับ:, , , , , ,
ติดตาม

Get every new post delivered to your Inbox.

Join 51 other followers