- รูปด้านล่างแสดงตัว object ที่วางลงบนฟอร์มประกอบด้วย 2 combobox เก็บ datum ซึ่งผมจำกัดไว้แค่เท่าที่ประเทศไทยใช้กัน ส่วนตัวอื่นๆ ก็ไม่มีอะไรมาก มี Textbox จำนวน 4 ตัวสำหรับให้ผู้ใช้โปรแกรมกรอกค่าพิกัดที่ต้องการแปลง และ Label จำนวน 4 ตัวเพื่อแสดงค่าพิกัดที่แปลงแล้ว และสุดท้ายคือมี Button จำนวน 2 ตัวเอาไว้แสดงระบบพิกัดที่ดูแล้วเข้าใจง่าย และเรียกรูปแบบการแสดงผลนี้ว่า Well Known Text (WKT)
ระบบค่าพิกัดในรูปแบบ EPSG (Europian Petrolium Survey Group)
- ระบบพิกัด (Coordinate Reference System) ที่ GDAL สนับสนุนนั้นมีหลายรูปแบบเช่นรูปแบบของ ESRI, EPSG และรูปแบบ Proj.4 รวมทั้งสามารถกำหนดเองได้ด้วยรูปแบบที่ผมกล่าวไปข้างต้นคือ Well Known Text แต่ในความรู้สึกของผมแล้วของ EPSG นั้นง่ายที่สุด
- ก่อนจะไปต่อขอพูดถึง EPSG ว่าเป็นสมาคมที่ตั้งในปี 1985 นำโดย Jean-Patrick GIRBIG ด้วยบุคคลที่เกี่ยวข้องกับกลุ่มอุตสาหกรรมพลังงาน (oil & gas) แถวๆยุโรป ก่อนหน้านี้ระบบพิกัดที่ใช้กันทั่วโลกนั้นหลากหลายมาก ความแตกต่างระหว่างทรงรี ความแตกต่างของเส้นโครงแผนที่ (Map Projection) ทำให้ยากต่อการใช้งาน สมาคมนี้จึงได้กำหนดมาตรฐานขึ้นมาโดยใช้ตัวเลข (code) เป็นตัวอ้างอิง ด้านล่างเป็น EPSG ที่เกี่ยวข้องกับประเทศไทย ซึ่งผมจะนำมาใช้ในโปรแกรมนี้ด้วย ตัวที่สำคัญและใช้มากก็คือ 4326 ก็คือ WGS84 แและเป็นค่าพิกัด geographic (latitude/longitude)
4326 : WGS84 (Geographic)
32647 : WGS84 / UTM Zone 47 North
32648 : WGS84 / UTM Zone 48 North
4240 : INDIAN 1975 (Geographic)
24047 : INDIAN 1975 / UTM Zone 47 North
24048 : INDIAN 1975 / UTM Zone 48 North
- แต่การใช้ EPSG ก็มีปัญหาเล็กน้อยถ้าเป็น Indian 1975 คือค่าพารามิเตอร์ที่ใช้ในการ transformation จะไม่ตรงกับที่เราใช้งานกันแต่ GDAL เปิดโอกาสให้เราแก้ไขได้ด้วยฟังก์ชั่นง่ายๆ
Sourcecode ของเมนฟอร์ม
- มาดูโค๊ดของฟอร์ม (mainform)
unit main;
{$mode objfpc}{$H+}
interface
uses
Classes, SysUtils, LResources, Forms, Controls, Graphics, Dialogs,
StdCtrls, gdal, gdalcore, ogr, ogrcore;
type
{ TfrmMain }
TCharSet = set of char;
TfrmMain = class(TForm)
cmdChkDatumFrom: TButton;
cmdChkDatumTo: TButton;
cmdCompute: TButton;
cmbDatumFrom: TComboBox;
cmbDatumTo: TComboBox;
GroupBox1: TGroupBox;
GroupBox2: TGroupBox;
Label1: TLabel;
Label10: TLabel;
labNorthingTo: TLabel;
Label2: TLabel;
Label3: TLabel;
Label4: TLabel;
Label5: TLabel;
Label6: TLabel;
Label7: TLabel;
Label8: TLabel;
Label9: TLabel;
labEastingTo: TLabel;
labLatitudeTo: TLabel;
labLongitudeTo: TLabel;
txtEastingFrom: TEdit;
txtLatitudeFrom: TEdit;
txtLongitudeFrom: TEdit;
txtNorthingFrom: TEdit;
procedure cmdChkDatumFromClick(Sender: TObject);
procedure cmdChkDatumToClick(Sender: TObject);
procedure cmbDatumFromChange(Sender: TObject);
procedure cmbDatumToChange(Sender: TObject);
procedure cmdComputeClick(Sender: TObject);
procedure FormActivate(Sender: TObject);
procedure FormCreate(Sender: TObject);
private
{ private declarations }
FEPSGFrom, FEPSGTo : longint;
procedure Tokenize(tokens: TStrings; const sztarget: string; szdelims: TCharSet);
procedure CheckDatum (epsgcode : longint);
public
{ public declarations }
end;
var
frmMain: TfrmMain;
implementation
procedure TfrmMain.Tokenize(tokens: TStrings; const sztarget: string; szdelims: TCharSet);
var
Start: Integer;
i: Integer;
StrLength: Integer;
begin
Start := 1;
StrLength := Length(sztarget);
for i := 1 to StrLength do
begin
if sztarget[i] in szdelims then
begin
if Start <= i - 1 then
tokens.Add(Copy(sztarget, Start, i - start));
Start := i + 1;
end;
end;
if not (sztarget[strLength] in szdelims) then
begin
tokens.Add(Copy(sztarget, Start, strLength - start + 1));
end;
end;
procedure TfrmMain.CheckDatum (epsgcode : longint);
var
pszWKT : PCHAR;
hSRS : TOGRSpatialReferenceH;
err : TOGRErr;
i : longint;
t1 : string;
tokens : TStrings;
begin
pszWKT := NIL;
hSRS := ogrcore.OSRNewSpatialReference(NIL);
err := ogrcore.OSRImportFromEPSG(hSRS, epsgcode);
if (err = OGRERR_NONE) then
begin
err := ogrcore.OSRExportToPrettyWkt(hSRS, @pszWKT, 0);
try
tokens := TStringList.Create;
tokenize(tokens, pszWKT, [#10]);
t1 := '';
for i := 0 to tokens.Count - 1 do
t1 := t1 + tokens.Strings[i] + #10;
ShowMessage (t1);
finally
tokens.Free;
end;
OGRFree(pszWKT);
end;
ogrcore.OSRDestroySpatialReference(hSRS);
end;
procedure TfrmMain.cmdChkDatumFromClick(Sender: TObject);
begin
CheckDatum(FEPSGFrom);
end;
procedure TfrmMain.cmdChkDatumToClick(Sender: TObject);
begin
CheckDatum(FEPSGTo);
end;
procedure TfrmMain.cmbDatumFromChange(Sender: TObject);
begin
if (cmbDatumFrom.ItemIndex = 0) then
begin
txtNorthingFrom.Text := '';
txtEastingFrom.Text := '';
labLatitudeTo.Caption := '';
labLongitudeTo.Caption := '';
end;
//Set EPSG Code to match with Combobox of datums.
case cmbDatumFrom.ItemIndex of
0 : FEPSGFrom := 4326; //WGS84 (Geographic)
1 : FEPSGFrom := 32647; //WGS84 / UTM Zone 47 North
2 : FEPSGFrom := 32648; //WGS84 / UTM Zone 48 North
3 : FEPSGTo := 4240; //INDIAN 1975 (Geographic)
4 : FEPSGFrom := 24047; //INDIAN 1975 / UTM Zone 47 North
5 : FEPSGFrom := 24048; //INDIAN 1975 / UTM Zone 48 North
end;
end;
procedure TfrmMain.cmbDatumToChange(Sender: TObject);
begin
if (cmbDatumTo.ItemIndex = 0) then
begin
labNorthingTo.Caption := '';
labEastingTo.Caption := '';
txtLatitudeFrom.Text := '';
txtLongitudeFrom.Text := '';
end;
case cmbDatumTo.ItemIndex of
0 : FEPSGTo := 4326;
1 : FEPSGTo := 32647;
2 : FEPSGTo := 32648;
3 : FEPSGTo := 4240;
4 : FEPSGTo := 24047;
5 : FEPSGTo := 24048;
end;
end;
procedure TfrmMain.cmdComputeClick(Sender: TObject);
var
hSrc, hTarget : TOGRSpatialReferenceH;
hCT : TOGRCoordinateTransformationH;
dX, dY, dZ : double;
err : TOGRErr;
//EPSGFrom, EPSGTo : longint;
success : longint;
retval : longint;
wkt : PCHAR;
//adCoeffs : TDoubleArray7;
begin
if (FEPSGFrom = FEPSGTo) then
begin
Showmessage('Can not convert coordinates at the same datum.');
exit;
end;
try
hSrc := ogrcore.OSRNewSpatialReference(NIL);
try
hTarget := ogrcore.OSRNewSpatialReference(NIL);
err := ogrcore.OSRImportFromEPSG(hSrc, FEPSGFrom);
err := ogrcore.OSRImportFromEPSG(hTarget, FEPSGTo);
//For Thailand we change transformation parameters to Royai Thai Survey Department (RTSD).
//Indian 1975 to WGS84 => TX=206, TY=837, TZ=295
if (ogrcore.OSRIsProjected(hSrc) = 1) then
err := ogrcore.OSRSetTOWGS84(hSrc, 206, 837, 295, 0, 0, 0, 0);
if (ogrcore.OSRIsProjected(hTarget) = 1) then
err := ogrcore.OSRSetTOWGS84(hTarget, 206, 837, 295, 0, 0, 0, 0);
//We can check the transformation parameters by this function.
//err := ogrcore.OSRGetTOWGS84(hSrc, @adCoeffs, 7);
success := 0;
try
hCT := ogrcore.OCTNewCoordinateTransformation(hSrc, hTarget);
if (ogrcore.OSRIsGeographic(hSrc) = 1) then
begin
dY := strtofloat(txtLatitudeFrom.Text);
dX := strtofloat(txtLongitudeFrom.Text);
dZ := 0;
end
else
begin
dY := strtofloat(txtNorthingFrom.Text);
dX := strtofloat(txtEastingFrom.Text);
dZ := 0;
end;
retval := ogrcore.OCTTransformEx(hCT, 1, @dX, @dy, @dZ, @success );
if( hCT = NIL) or ( retval = 0) then
Showmessage('Computation failed.\n')
else
begin
if (ogrcore.OSRIsGeographic(hTarget) = 1) then
begin
labLatitudeTo.Caption := format('%10.6f', [dY]);
labLongitudeTo.Caption := format('%10.6f', [dX]);
end
else
begin
labNorthingTo.Caption := format('%10.3f', [dY]);
labEastingTo.Caption := format('%10.3f', [dX]);
end;
end;
finally
ogrcore.OCTDestroyCoordinateTransformation(hCT);
end;
finally
ogrcore.OSRDestroySpatialReference(hTarget);
end;
finally
ogrcore.OSRDestroySpatialReference(hSrc);
end;
end;
procedure TfrmMain.FormActivate(Sender: TObject);
begin
cmbDatumFromChange(Sender);
cmbDatumToChange(Sender);
end;
procedure TfrmMain.FormCreate(Sender: TObject);
begin
gdalcore.AllRegister;
//Set the path to GDAL Library.
//If the path is not correct. We can't use EPSG database.
gdalcore.SetConfigOption('GDAL_DATA', 'C:\gdalwin32-1.6\data');
end;
initialization
{$I main.lrs}
end.- ผมจะอธิบายโค๊ดพอสังเขป เริ่มจากตรง uses clause ในโปรแกรมนี้จะใช้ไฟล์ pascal 4 ไฟล์คือ gdal.pas, gdalcore.pas, ogr.pas และ ogrcore.pas และประกาศตัวแปรตรง private section ได้แก่ตัวแปร FEPSGFrom และ FEPSGTo ไว้เก็บตัวเลข EPSG ส่วนโปรซีเจอร์ Tokenize เอาไว้แยกสตริงที่ใช้ตัวคั่น (ตัวอย่างไฟล์ csv ที่ใช้เครื่องหมายคอมม่าคั่น) ส่วนฟังก์ชั่น CheckDatum เอาไว้แสดงระบบพิกัดเมื่อส่งตัวเลข EPSG มาให้
private
{ private declarations }
FEPSGFrom, FEPSGTo : longint;
procedure Tokenize(tokens: TStrings; const sztarget: string; szdelims: TCharSet);
procedure CheckDatum (epsgcode : longint);- เมื่อเริ่มสร้างฟอร์มเราจะตั้งค่าให้ GDAL ได้ทราบ path ที่ติดตั้งไว้ โดยไปดักที่อีเวนต์ FromCreate ผมสังเกตว่าในวินโดส์ต้องกำหนดให้ผ่านฟังก์ชั่น SetConfigOption ส่วนบรรทัดแรกนั้นสั่งให้ไลบรารีทำการ register ทุกฟอร์แม็ตที่สนับสนุน
procedure TfrmMain.FormCreate(Sender: TObject);
begin
gdalcore.AllRegister;
//Set the path to GDAL Library.
//If the path is not correct. We can't use EPSG database.
gdalcore.SetConfigOption('GDAL_DATA', 'C:\gdalwin32-1.6\data');
end;ฟังก์ชั่นคำนวณ
- การคำนวณจะเรียกใช้ไฟล์ฟังก์ชั่นที่อยู่ใน ogrcore.pas (ส่วน gdalcore.pas จะเกีียวข้องกับ raster) มาดูโค๊ด
procedure TfrmMain.cmdComputeClick(Sender: TObject);
var
hSrc, hTarget : TOGRSpatialReferenceH;
hCT : TOGRCoordinateTransformationH;
dX, dY, dZ : double;
err : TOGRErr;
//EPSGFrom, EPSGTo : longint;
success : longint;
retval : longint;
wkt : PCHAR;
//adCoeffs : TDoubleArray7;
begin
if (FEPSGFrom = FEPSGTo) then
begin
Showmessage('Can not convert coordinates at the same datum.');
exit;
end;
try
hSrc := ogrcore.OSRNewSpatialReference(NIL);
try
hTarget := ogrcore.OSRNewSpatialReference(NIL);
err := ogrcore.OSRImportFromEPSG(hSrc, FEPSGFrom);
err := ogrcore.OSRImportFromEPSG(hTarget, FEPSGTo);
//For Thailand we change transformation parameters to Royai Thai Survey Department (RTSD).
//Indian 1975 to WGS84 => TX=206, TY=837, TZ=295
if (ogrcore.OSRIsProjected(hSrc) = 1) then
err := ogrcore.OSRSetTOWGS84(hSrc, 206, 837, 295, 0, 0, 0, 0);
if (ogrcore.OSRIsProjected(hTarget) = 1) then
err := ogrcore.OSRSetTOWGS84(hTarget, 206, 837, 295, 0, 0, 0, 0);
//We can check the transformation parameters by this function.
//err := ogrcore.OSRGetTOWGS84(hSrc, @adCoeffs, 7);
success := 0;
try
hCT := ogrcore.OCTNewCoordinateTransformation(hSrc, hTarget);
if (ogrcore.OSRIsGeographic(hSrc) = 1) then
begin
dY := strtofloat(txtLatitudeFrom.Text);
dX := strtofloat(txtLongitudeFrom.Text);
dZ := 0;
end
else
begin
dY := strtofloat(txtNorthingFrom.Text);
dX := strtofloat(txtEastingFrom.Text);
dZ := 0;
end;
retval := ogrcore.OCTTransformEx(hCT, 1, @dX, @dy, @dZ, @success );
if( hCT = NIL) or ( retval = 0) then
Showmessage('Computation failed.\n')
else
begin
if (ogrcore.OSRIsGeographic(hTarget) = 1) then
begin
labLatitudeTo.Caption := format('%10.6f', [dY]);
labLongitudeTo.Caption := format('%10.6f', [dX]);
end
else
begin
labNorthingTo.Caption := format('%10.3f', [dY]);
labEastingTo.Caption := format('%10.3f', [dX]);
end;
end;
finally
ogrcore.OCTDestroyCoordinateTransformation(hCT);
end;
finally
ogrcore.OSRDestroySpatialReference(hTarget);
end;
finally
ogrcore.OSRDestroySpatialReference(hSrc);
end;
end;- จะเห็นว่าเมื่อผู้ใช้โปรแกรมคลิกที่ปุ่มคำนวณจะเรียกอีเวนต์นี้ ต้นฟังก์ชั่นผมประกาศตัวแปรสำหรับการแปลงค่าพิกัดคือ ซึ่งทั้งหมดเป็น pointer หรือเรียกอีกนัยหนึ่งว่า handle
hSrc, hTarget : TOGRSpatialReferenceH; hCT : TOGRCoordinateTransformationH;- ซึ่ง hSrc และ hTarget จะได้ค่า handle จากการเรียกฟังก์ชั่น OGR คือ OSRNewSpatialReference(NIL) เมื่อได้ค่า handle แล้วระบบพิกัดต้นทางและปลายทางมีตัวตนขึ้นมาแล้ว ต่อไปจะกำหนดค่าตัวเลข EPSG ให้ผ่านฟังก์ชั่น OSRImportFromEPSG(handle, EPSGcode) ซึ่งตัวแปร FEPSGFrom และ FEPSGTo จะได้ตัวเลขจากการ map ค่ากับ combobox ที่แสดง datum
- ที่ผมกล่าวไปแล้วว่าค่าพารามิเตอร์ที่ใช้ในการแปลงค่าพิกัดไม่ตรงกับที่เราใช้กัน ผม assign ค่าให้ใหม่ด้วย OSRSetTOWGS84(hSrc, 206, 837, 295, 0, 0, 0, 0) และกำหนดให้ hTarget เช่นเดียวกัน
- handle ตัว hSrc และ hTarget เก็บระบบพิกัดเท่านั้นจะคำนวณแปลงค่าพิกัดด้วยตัวเองไม่ได้ การคำนวณจะใช้ handle อีกตัวคือ hCT โดยการเรียกฟังก์ชั่น hCT := ogrcore.OCTNewCoordinateTransformation(hSrc, hTarget) เพื่อสร้าง handle ขึ้นมา ต่อไปใช้ฟังก์ชั่น OSRIsGeographic(hSrc) เพื่อตรวจสอบว่าระบบพิกัดที่ผู้ใช้เลือกว่าเป็น geographic หรือเป็น projection (UTM) จะได้เลือกเก็บค่าพิกัดที่ผู้ใช้ป้อนได้ถูกที่ (ผมไม่ได้เขียนดัก error ที่อาจจะเกิดจากการป้อนตัวอักษรแทนตัวเลข ถ้าจะเขียนโปรแกรมไปใช้เป็นเรื่องเป็นราวต้องเขียนดัก)
- ทำการคำนวณการแปลงพิกัดด้วยฟังก์ชั่น retval := ogrcore.OCTTransformEx(hCT, 1, @dX, @dy, @dZ, @success ) ด้วยการส่ง pointer ค่าพิกัด ตรวจสอบผลลัพธ์การคำนวณด้วย retval หรือ success ว่าเท่ากับ 1 หรือไม่ ฟังก์ชั่นตัวนี้สามารถคำนวณค่าพิกัดได้มากกว่า 1 จุด ด้วยการส่ง array of double ไปแทน dX, dY, dZ ถ้าดูในฟังก์ชั่นจะเห็นพารามิเตอร์ที่สอง จะเป็น 1 คือจำนวนจุดนั่นเอง ถ้ามีสัก 10 จุดที่ต้องการแปลงต้องส่งเลข 10 ไปให้ และอาร์เรย์จะต้องประกาศเป็น array[0..9] of double ส่งไปให้
ทอสอบโปรแกรมบน windows
- พร้อมแล้วก็กดคีย์ F9 เพื่อรันโปรแกรมได้เลย ผมลองแปลงพิกัดจาก WGS84 /UTM Zone 47N ไปยัง Indian 1975/UTM Zone 47N ดังรูปด้านล่าง
ทดสอบโปรแกรมบน Ubuntu linux
- รัน lazarus แล้วเปิดโครงการที่จำเพาะกับ linux มาลองรันดูแล้วป้อนค่าพิกัดเข้าไป
- เพื่อยืนยันว่าการแปลงพิกัดได้ถูกต้อง ถ้าผู้อ่านมี Geocalc บนวินโดส์ก็สามารถทดสอบได้ค่าจะได้เท่ากันครับ ต่อไปลองคลิกที่ปุ่ม Coordinate System Information ดูเพื่อดูระบบพิกัดในรูปแบบ Well Known Text
- เป็นอย่างไรบ้างครับ ความสามารถของ GDAL ไม่ธรรมดาเลยและยัง cross-platform ด้วย ผมดูแล้วถ้ามานั่งเขียนให้ support ทุกๆ datum ก็ไม่น่ายากเกินไปเพราะไลบรารีเตรียมมาให้หมดแล้ว ลองศึกษาที่ website ของ GDAL ดู ผมเองยังดูได้ไม่หมดเลยเพราะไลบรารีค่อนข้างใหญ่ แต่ก็ไม่ยากจนเกินไปครับ