การเขียนโปรแกรมคำนวณแปลงพิกัดระหว่าง Datum ด้วยไลบรารี GDAL/OGR (ตอนที่ 2)

17/01/2010 Leave a comment
  • รูปด้านล่างแสดงตัว 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 ดังรูปด้านล่าง
แปลงพิกัดจาก WGS84 ไปยัง Indian 1975

ทดสอบโปรแกรมบน Ubuntu linux

  • รัน lazarus แล้วเปิดโครงการที่จำเพาะกับ linux มาลองรันดูแล้วป้อนค่าพิกัดเข้าไป
แปลงค่าพิกัดทดสอบใน linux
  • เพื่อยืนยันว่าการแปลงพิกัดได้ถูกต้อง ถ้าผู้อ่านมี Geocalc บนวินโดส์ก็สามารถทดสอบได้ค่าจะได้เท่ากันครับ ต่อไปลองคลิกที่ปุ่ม Coordinate System Information ดูเพื่อดูระบบพิกัดในรูปแบบ Well Known Text
แสดงระบบพิกัดแบบ WKT
  • เป็นอย่างไรบ้างครับ ความสามารถของ GDAL ไม่ธรรมดาเลยและยัง cross-platform ด้วย ผมดูแล้วถ้ามานั่งเขียนให้ support ทุกๆ datum ก็ไม่น่ายากเกินไปเพราะไลบรารีเตรียมมาให้หมดแล้ว ลองศึกษาที่ website ของ GDAL ดู ผมเองยังดูได้ไม่หมดเลยเพราะไลบรารีค่อนข้างใหญ่ แต่ก็ไม่ยากจนเกินไปครับ

Leave a Reply

Your email address will not be published. Required fields are marked *