ตุลาคม | 2009 | Prajuab Riabroy's Blog

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

30/10/2009 prajuab riabroy 5 ของความคิดเห็น

ข้อดีของ 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 แล้วเรียกใช้คำสั่งดังรูปด้านล่างเพื่อเรียกโปรแกรม

ผมเปิดฐานข้อมูลด้วยคำสั่ง sqlite3 e:\sourcecodes\lazarus\sqliteproj1\datums.s3db ถ้าสงสัยหรือจำคำสั่งไม่ได้ก็ให้พิมพ์ .help เข้าไป
ผมใช้คำสั่ง .databases เพื่อดูว่า database อะไรที่เราใช้งานอยู่ ตามด้วย .header ON เพื่อให้แสดงคอลัมภ์ชื่อฟิลด์ และ .show เพื่อแสดงพารามิเตอร์ แล้วลอง query ดูด้วย syntax “select” ดังรูปด้านล่าง

จากรูปด้านบนจะเห็นผลลัพธ์ของ sql ด้วย syntax “select” ที่จริงยังมีคำสั่งอีกมากที่ทำได้ในสภาวะ command line แต่ตอนนี้ขอพอแค่นี้ก่อน ตอนหน้าจะลงลึกไปอีกนิดว่าจะนำ SQLite มาใช้กับ Lazarus อย่างไร

Categories: Linux, Programming, Windows ป้ายกำกับ:database, Lazarus, Programming, sqlite

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

25/10/2009 prajuab riabroy ใส่ความเห็น

ตอนที่แล้วได้แนะนำสูตรที่จะใช้ในการคำนวณและแสดงยูนิต 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 ดังรูปด้านล่าง

ที่ code คลิกเลือก Target OS เป็น win32/win64 และ Target CPU Family เป็น i386 ส่วน Target processor เลือกเป็น Default ดังรูปด้านล่าง

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

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

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

คลิกที่เมนู Project > Save Project As… ตั้งชื่อเป็น Geodesics จะเห็นหน้าตาของ Lazarus และ project Geodesics ดังรูปด้านล่าง

ต่อไปจะแปะพวก object ทั้งหลายเฃ่น Label, Edit, ComboBox, Button ลงไปและ ตั้งฃื่อ (Name) ดังรูปด้านล่าง

ส่วน 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

Categories: Lazarus, Linux, Programming, Surveying, Windows ป้ายกำกับ:คำนวณ, ระะทาง, อะซิมัท, เขียนโปรแกรม, Compute, direct, geodesic, inverse, Lazarus, vincenty

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

20/10/2009 prajuab riabroy ใส่ความเห็น

ใน 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 ให้ถูกต้อง

จากรูปด้านบนที่ 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

โปรแกรม Geocalc เรียกว่า Map System ท่านผู้อ่านอาจจะพบคำว่า Coordinate Systems มากกว่า มีความหมายคือระบบพิกัดนั่นเอง

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

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

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

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

ที่เมนูหลักคลิก Process > Manual Input

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

โจทย์ของผมวันนี้คือเรามีค่าพิกัด UTM อยู่ในระบบ Indian 1975 มีค่า Easting(E) = 391024.838 ค่า Northing(N) = 1576384.958 ต้องการแปลงเป็นค่า Lat/Long ในระบบพิกัด 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

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

ผลลัพธ์ของการจุดค่าพิกัดบน 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 ก็คำนวณได้

Categories: GIS, GPS, Surveying, Windows ป้ายกำกับ:ลองจิจูด, แปลงพิกัด, แลตติจูด, Calculate, Comput, Compute, Coordinate, Datum, GeoCalc, Google Earth, Latitude, Longitude, Transform

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

17/10/2009 prajuab riabroy ใส่ความเห็น

ตอนก่อนหน้านี้ ผมเขียนโปรแกรมแปลงพิกัดระหว่าง 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 อ้างอิงจาก 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 (a−b)/a
φ₁, φ₂ = geodetic latitude
L = difference in longitude
U₁ = atan((1−f).tanφ₁) (U is ‘reduced latitude’)
U₂ = atan((1−f).tanφ₂)
λ = L (first approximation)
iterate until change in λ is negligible (e.g. 10^-12 ≈ 0.06mm) {
	sinσ = √[ (cosU₂.sinλ)² + (cosU₁.sinU₂ − sinU₁.cosU₂.cosλ)² ]	(14)
	cosσ = sinU₁.sinU₂ + cosU₁.cosU₂.cosλ	(15)
	σ = atan2(sinσ, cosσ)	(16)
	sinα = cosU₁.cosU₂.sinλ / sinσ	(17)
	cos²α = 1 − sin²α (trig identity; §6)
	cos2σ_m = cosσ − 2.sinU₁.sinU₂/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)
α₁ = atan2(cosU₂.sinλ, cosU₁.sinU₂ − sinU₁.cosU₂.cosλ)		(20)
α₂ = atan2(cosU₁.sinλ, −sinU₁.cosU₂ + cosU₁.sinU₂.cosλ)		(21)

where :

- s is the distance (in the same units as a & b)
- α₁ is the initial bearing, or forward azimuth
- α₂ is the final bearing (in direction p₁→p₂)

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

ยูนิต GeoEllipsoids (geoellipsoids.pas) เก็บสัณฐานทรงรีที่ใช้กันทั่วโลก ที่ผมเขียนไว้เรื่อง การเขียนโปรแกรมคำนวณการแปลงพิกัดระหว่าง UTM และ Geographic (ตอนที่ 2)

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

ผมจะสร้างยูนิตที่ทำหน้าที่คำนวณชื่อ 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 ดูรูปผลลัพธ์ของโปรแกรมที่ผมเขียนในตอนนี้

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

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

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

ทิ้งท้ายกันก่อนในตอนนี้ ติดตามกันตอนหน้า

Categories: Linux, Programming, Windows ป้ายกำกับ:cross-platform, ellipsoid, geodesic, Latitude, Lazarus, Longitude

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

12/10/2009 prajuab riabroy ใส่ความเห็น

ตัว 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)

จากรูปด้านบนจะเห็นจุดที่เรียงเป็นกริดนั้น 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.

จะมีไดอะล็อก ให้ตั้งเงื่อนไขของ 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

จากรูปด้านบนเมื่อกด 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

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

ตั้งค่า GeoTiff

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

Viewer by Erdas Imagine

Categories: GIS, Windows ป้ายกำกับ:3D, Create, DEM, DTM, Generate, Global Mapper, Points

Older Entries

ตุลาคม 2009
จ	อ	พ	พฤ	ศ	ส	อา
« ก.ย.				พ.ย. »
			1	2	3	4
5	6	7	8	9	10	11
12	13	14	15	16	17	18
19	20	21	22	23	24	25
26	27	28	29	30	31

Prajuab Riabroy's Blog

Archive