Month: January 2010

การคำนวณความสูง Geoid โดยใช้ EGM96 และ EGM2008

  • ขอพักเรื่องโปรแกรมมิ่งสักตอนเพราะว่า เขียนเกี่ยวกับเรื่องโปรแกรมมิ่งต้องใช้พลังความคิดมาก นอกจาก library “GDAL” ที่ผมกำลังเขียนถึง ซึ่งมีแง่มุมให้เขียนเกี่ยวกับการใช้งานได้เป็นร้อยๆตอนเลยละครับ ถ้ามีแรงกายและใจขนาดนั้น ถ้ายังจำกันได้คือ GeographicLib ของ Charles Karney ที่เคยอ้างอิงถึงไปแล้ว มีเรื่องอีกเรื่องที่น่าสนใจคือ การคำนวณ Geoid Height ซึ่งดูโค๊ดของ Charles Karney โมดูลนี้ต้องเป็นงงเรื่องอัลกอริทึ่ม มากไปกว่านั้นมีรุ่นน้องส่งโค๊ดมาให้คำนวณเรื่องเดียวกันแต่เป็นโค๊ดของภาษา Fortran ยิ่งงงเข้าไปอีก ทั้งที่เคยร่ำเรียนตอนอยู่มหาวิทยาลัย แต่ก็ลืม syntax ไปหมด
  • EGM (Earth Gravitational Model) ที่คนที่รังวัด GPS คงต้องทราบกันดี รุ่นก่อนหน้านี้คือ EGM96 ใช้งานแพร่หลายที่สุด โปรแกรมคำนวณ GPS ทั้งหลายเช่น Leica Geo Office (LGO) หรือ Trimble Geo Office (TGO) ผมยังใช้ LGO5 ยังเป็นรุ่นเก่าอยู่ไม่ทราบว่าในรุ่นใหม่ ทั้ง TGO และ LGO ปรับมาใช้ EGM2008 แล้วหรือยัง
  • ที่อยากเขียนเรื่องคำนวณ Geoid Height เพราะในการใช้ GPS ในปัจจุบันการอ้างอิงความสูงจะเป็นความสูงบน Ellipsoid ก็คือ WGS84 ถ้าบอกว่าความสูงบน WGS84 เท่ากับ 35 เมตร จะเท่าไหร่เมื่อเมื่อเทียบกับ รทก. (ระดับน้ำทะเลปานกลาง หรือ MSL – Mean sea level) ที่เราคุ้นเคยกันดี การคำนวณ Geoid height จึงเป็นเรื่องที่คำนวณการทอนความสูงจาก WGS84 มาบน MSL (หรือเรียกอีกอย่างว่า Orthometric Height) แต่คนใช้ระดับชาวบ้าน บางครั้งไม่ต้องสนใจเพราะในเครื่อง GPS มือถือรุ่นใหม่ทั้งหลายได้คำนวณความสูงจากทรงรี WGS84 เป็น ให้เป็น MSL ให้เป็นที่เรียบร้อยแล้ว (คือภายในเครื่องได้บรรจุ EGM96 เข้าไปเรียบร้อย)

รูปทรง Geoid

  • รูปทรงรีที่ใช้แทนสัณฐานโลก จะมีรูปทรงที่แน่นอน สามารถใช้สูตรทางคณิตศาสตร์คำนวณได้ แต่ Geoid เป็นจะเป็นทรงที่ไม่แน่นอน บุบๆบู้บี้ ดูรูปข้างล่างประกอบ
รูปทรง Geoid (ภาพจาก http://en.wikipedia.org/wiki/Geoid)

ความสัมพันธ์ระหว่าง Geoid และ ทรงรี

  • เนื่องจากค่าความสูงที่เราต้องการคือต้องอยู่บน Geoid (เทียบเท่าระดับน้ำทะเลปานกลาง) แต่ความสูงที่ได้จากรังวัด GPS (ตัวอย่างเช่นการรังวัดแบบ Fast Static และ Static) จะอยู่บนทรงรี ซึ่งแต่ละสถานที่ความต่างเนื่องรูปทรงที่ไม่แน่นอนของ Geoid จึงทำให้แต่ละสถานทีค่าความสูงต่างจากทรงรีและบน Geoid จะไม่เท่ากัน
ความสัมพันธ์ Geoid และทรงรี (ภาพจาก http://principles.ou.edu/earth_figure_gravity/geoid/index.html)
  • จากรูปด้านบน(ผมชอบรูปนี้มากแสดงได้ชัดเจน) ค่าความสูง H (รทก.) = h (ความสูงบนทรงรี ได้จาก GPS) – N (ค่าความสูง Geoid) จากสูตรความสูงของจุดที่เราต้องการจริงๆแล้วอยู่บนภูมิประเทศ (Topography หรือ Terrain) ถ้าเรารังวัด GPS จะได้ค่า h ออกมา ถ้าสามารถคำนวณหา N ได้จากค่าพิกัด geographic (lat/long) เราก็สามารถหาค่า H ที่เราต้องการได้ จากรูปด้านบนจะสังเกตเห็นว่าแนวคำนวณหา H จะเป็นแนวตั้งฉากกับ Geoid ซึ่งจะเป็นแนวแรงไปตามแรงโน้มถ่วงเสมอ แต่แนว h จะเป็นแนวตั้งฉากกับทรงรี ลองดูรูปอีกรูปหนึ่งด้านล่าง
ความสัมพันธ์ระหว่าง Geoid และทรงรี (ภาพจาก http://en.wikipedia.org/wiki/Geoid)

ความเป็นมาของ EGM (Earth Gravitational Model)

  • EGM96 เป็นผลงานของความร่วมมือของ National Imagery and Mapping Agency (NIMA), NASA Goddard Space Flight Center (GSFC) และ Ohio State University โครงการนี้ได้รวบรวมข้อมูลความโน้มถ่วงจากหลายๆแหล่งของโลก ด้วยวิธีที่แตกต่างกันบ้าง เช่นในมหาสมุทรใช้ดาวเทียม GEOSAT และ ERS-1 ข้อมูลเหล่านี้จะนำมาหาค่าที่เรียกว่า coefficients ของ EGM96
  • แรงโน้มถ่วงของโลกมีลักษณะเป็นฟังก์ชั่น Harmonic  ซึ่งฟังก์ชั่นฮาร์โมนิค จะคำนวณได้ก็ต่อเมื่อทราบค่า coefficients ที่กล่าวไว้ข้างต้นนั่นเอง โมเดล EGM96 จะประกอบไปด้วยกริดแต่ละช่องขนาด 15 ลิปดา x 15 ลิปดา และแต่ละช่องจะเก็บค่า coefficients ที่แตกต่างกันไปตามค่าพิกัด การคำนวณหาความสูง geoid (Geoid Undulation) จะรับ Input จากผู้ใช้คือค่าพิกัด แล้วนำค่าพิกัดไปดึงเอาค่า coefficient แล้วนำค่า coefficient ไปคำนวณหาค่าความสูง Geoid (N – เรียกอีกอย่างว่า Geoid separation)

EGM2008

  • จัดทำโดย U.S. National Geospatial-Intelligence Agency (NGA) ชื่อเดิมก็คือ NIMA นั่นเอง โดยได้ปรับปรุงจาก EGM96 โดยเพิ่มข้ิอมูลจากดาวเทียม Gravity Recovery and Climate Experiment(GRACE) เป็นดาวเทียมของ NASA ที่วัดสนามความโน้มถ่วงของโลกปล่อยสู่วงโคจรในปี 2002 โดยที่โมเดลนี้ได้เผยแผ่ในปี 2008 จึงเรียกว่า EGM2008 ความละเอียดของ EGM2008 ประกอบด้วยกริดที่แต่ละช่องมีขนาด 1′ x 1′  หรือประมาณบรรจุค่า coefficients ประมาณ 4 ล้านค่า

EGM96 vs. EGM2008

  • EGM2008 ควรจะมี accuracy สูงกว่า EGM96 แต่บางรายงานก็ยังก้ำกึ่ง ตัวอย่างรายงานและการทดสอบที่แสดงว่า EGM2008 นั้นดีกว่า EGM96 ได้แก่ Wuhan University หรือที่ โดยศาสตราจารย์  Charles Merry, University of Cape Town ท้ายข้อสรุประบุว่ายากที่จะบอกว่า EGM2008 ละเอียดกว่า EGM96 มากเท่าไหร่ แต่ที่ทวีปอเมริกาเหนือ ทวีปยุโรป ค่า Accuracy ของ EGM2008 ดีถึง ±10 ซม.

การคำนวณความสูง Geoid บน EGM96

  • หาโปรแกรมที่เป็นโปรแกรม standalone ไม่ได้ส่วนใหญ่แล้วจะแฝงอยู่ในโปรแกรมด้านการแปลงพิกัด เช่น Trimble Coordinate Calculator หรือโปรแกรมด้านการคำนวณรังวัด GPS เช่น TGO และ LGO นอกเหนือจากนั้นเป็นโปรแกรม online ซะมากกว่า ที่แรกก็คือที่ NGA โดยตรงไปที่ http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/intpt.html ดูตัวอย่างด้านล่าง
ป้อนข้อมูลคำนวณ
  • ผลลัพธ์ที่ได้
ผลลัพธ์ที่ได้คือค่าความสูง Geoid
  • สมมติว่าตรงที่ผมป้อนค่าพิกัดไป รังวัด GPS ได้ความสูงบนทรงรีเท่ากับ (h) 2.504 เมตร และค่า N จากรูปด้านบนเท่ากับ -34.39 เมตร จากสูตร H = h – N จะได้ค่าความสูงที่เป็นรทก. H = 2.504 – – 34.39 = 2.504 + 34.39 = 36.894 เมตร

ดาวน์โหลดและติดตั้งโปรแกรม Alltrans EGM2008 Calculator

  • เป็นโปรแกรมฟรี สามารถดาวน์โหลดได้ที่ alltransegm2008.zip พัฒนาโดย Hans-Gerd Duenck-Kerst เป็นโปรแกรมเมอร์ชาวเยอรมัน ในตัวโปรแกรมจะมีค่า coefficients ของ  EGM2008 มาด้วยขนาด 10′ x 10′ ซึ่งเวลาดาวน์โหลด zip file ของโปรแกรมนี้จะมีขนาดเกือบๆ 7 MB แต่ถ้าติดตั้งแล้วไปดูขนาดของโปรแกรม ขนาดเพียงแค่ 80 KB เขียนด้วย MS C++ จะเห็นอีกไฟล์ชื่อ Und_min10x10_egm2008_isw=82_WGS84_TideFree_SE เป็นไฟล์ข้อมูลค่า coefficients ของ EGM2008 ขนาดเมื่อ unzip แล้วประมาณ 9 MB
  • ในหน้า website ของ NGA ที่ http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008/egm08_wgs84.html จะเห็นว่ามีหมวดหมู่ของโปรแกรมเช่นการคำนวณความสูง Geoid ด้วยวิธี Interpolation หมวดหมู่แรกคือ

The 1 x 1 minute Grid Interpolation

  • จะมีโปรแกรมที่เขียนด้วย Fortran ให้ดาวน์โหลด ซึ่งผมจะไม่แนะนำเพราะเป็น command line เกรงว่าจะยากไป ประกอบไปด้วยไฟล์ coefficients ของ EGM2008 ขนาดของกริด 1′ x 1′ ชื่อ Und_min1x1_egm2008_isw=82_WGS84_TideFree_SE.gz ไฟล์มีขนาดประมาณ 825 MB ซึ่งใหญ่มาก ซึ่งถ้าใช้ internet ที่ไม่เร็วเช่นผ่านมือถือซึ่งผมต้องใช้บ่อยเพราะนอกพื้นที่แล้วจะไม่มี high speed ให้ใช้ ดังนั้นผมยังไม่แนะนำตอนนี้ และใน website เองก็เตือนว่าถ้าใช้ไฟล์ข้อมูลตัวนี้จะกินแรมมาก หมวดหมู่ที่สองคือ

The 2.5 x 2.5 Minute Grid Interpolation

  • เช่นเดียวกันจะมีโปรแกรมที่เขียนด้วย Fortran เพื่อคำนวณไฟล์ค่า coefficients EGM2008 ขนาด 2.5 ลิปดา x 2.5 ลิปดา ชื่อไฟล์คือ Und_min2.5×2.5_egm2008_isw=82_WGS84_TideFree_SE.gz ผมแนะนำให้ดาวน์โหลดไฟล์นี้ ซึ่งมีขนาดประมาณ 135 MB ซึ่งพอรับได้ เราจะใช้ไฟล์นี้ไปใช้กับโปรแกรม Alltrans EGM2008 Calculator ที่แนะนำไปข้างต้น เมื่อดาวน์โหลดไฟล์นี้แล้วให้ unzip ไว้

AllTrans EGM2008 Calculator

  • ลองรันโปรแกรมดู จะเห็นโปรแกรมมีหน้าตาที่เรียบง่ายมาก
โปรแกรม Alltrans EGM2008 Calculator

 

 

    • จากรูปด้านบนจะมีเมนูที่เป็นแท็ปให้เลือกอยู่ 3 ทางเลือกหลัก คือ EGM2008 Manual Calc, EGM2008 GridMaker และ EGM2008 File Calc

 

    • ก่อนจะทำการคำนวณ เราจะตั้งไฟล์ข้อมูลของ EGM2008 ขนาด 2.5′ x 2.5′ ดังรูปด้านล่าง โดยการ browse ไปที่ไฟล์ที่ดาวน์โหลดมาและ unzip แล้วที่ผมกล่าวไปแล้วข้างต้น อย่าลืมคลิกเพื่อเลือก 2.5′ x 2.5′  Grid ให้สอดคล้องกับไฟล์ข้อมูลด้วย

 

 

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

  • ลองคำนวณแบบ manual ดู โดยใช้ค่าพิกัด latitude 12°30’10” (12.5027778) longitude 98°45’15” (98.7541667) ดูผลลัพธ์ แต่โปรแกรมจะแสดงผลลัพธ์ให้ด้วยอัลกอริทึ่มหลายๆวิธ๊
ผลลัพธ์การคำนวณแบบ manual

การสร้างกริดค่าความสูง Geoid

  • ลองป้อนค่าพิกัด ตั้งค่า spacing แล้วอย่าลืมคลิกที่ปุ่ม Calc จากนั้นคลิกที่ View ดูผลลัพธ์
สร้างกริดความสูง Geoid
  • ผลลัพธ์จะถูกเปิดด้วยโปรแกรม Notepad ของวินโดส์
Nodepad แสดงค่าความสูง Geoid เป็นกริด
  • ส่วนการคำนวณด้วยการ input เป็นไฟล์ ก็ไม่ยากลองดูโปรแกรม ผมไม่ได้แสดงตัวอย่างไว้
การคำนวณในโหมด input file
  • สรุปแล้ว accuracy ของ EGM2008 ก็ยังไม่ชัดเจนว่าจะดีกว่า EGM96 แต่เราได้ศึกษาการใช้ EGM2008 ก็ดีครับ ส่วนนักวิทยาศาสตร์ก็มีหน้าที่ต้องปรับปรุงโมเดลให้ดีไปกว่านี้ แต่ขนาด EGM96 ที่ผมใช้ในประเทศไทยเราพบว่าได้ค่าระดับเมื่อเทียบกับค่าระดับที่ได้จากงานระดับ (Levelling) อยู่ในเกณฑ์ประมาณเกือบๆงานชั้น 3 ของงานระดับ แต่ทั้งนี้ก็ไม่มีใครกล้าการันตี (ส่วนค่าพิกัดที่ได้จากรังวัด GPS แบบ fast static หรือ static นั้นยอดเยี่ยมอยู่แล้วสามารถการันตีได้) ส่วนผมแล้วก็เริ่มนำ EGM 2008 มาใช้แล้ว

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

  • รูปด้านล่างแสดงตัว 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 ดู ผมเองยังดูได้ไม่หมดเลยเพราะไลบรารีค่อนข้างใหญ่ แต่ก็ไม่ยากจนเกินไปครับ

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

  • ตอนก่อนนี้ ผมเขียนโปรแกรมทดสอบไลบรารี GDAL/OGR ด้วยการเปิดไฟล์รูปแล้วอ่าน Metadata, ระบบพิกัด ตลอดจนแปลงฟอร์แม็ตของไฟล์รูป จะเห็นถึงความสามารถของไลบรารี ที่เตรียมฟังก์ชั่นทุกสิ่งทุกอย่างครอบคลุมด้าน Geospatial ไว้พร้อมสรรพ และไม่ต้องแปลกใจที่โปรแกรมดังๆ เช่น Google Earth, ArcGIS, Quantum GIS ต่างก็นำไปใช้ ดูชื่อโปรแกรมที่นำไลบรารีไปใช้  ที่นี่ และการดาวน์โหลดและการติดตั้งอ่านที่ blog ของผมได้ ที่นี่
  • ครั้งก่อนผมเขียนโปรแกรมสำหรับคำนวณแปลงค่าพิกัดระหว่าง UTM และ Geographic (Latitude/Longitude) ด้วยการดัดแปลงโค๊ดภาษา c++ เป็น lazarus จาก GeographicLib ของคุณ Charles Karney ซึ่งมีข้อจำกัดแค่แปลงพิกัด UTM และ Geographic และต้องอยู่ใน datum เดียวกันเท่านั้น ความจริงในไลบรารี GeographicLib ได้เขียนไว้ครอบคลุมการแปลงค่าพิกัดระหว่าง datum ทั้งหมด แต่ด้วยข้อจำกัดภาษา c++ ของผมจึงได้แค่นั้น
  • เมื่อมาศึกษาไลบรารีของ GDAL/OGR พบว่ามีฟังก์ชั่นด้าน Geospatial ครบครันสมกับชื่อ GDAL(Geospatial Data Abstract Library) ที่ทางผู้พัฒนาตั้งชื่อไว้ ดูฟังก์ชั่นด้านการแปลงพิกัดของ GDAL พบว่าครอบคลุมทุกๆ datum ซึงข้อดีของการใช้ไลบรารี ก็คือเขียนโค๊ดน้อยมาก แต่ข้อเสียคือต้องศึกษาฟังก์ชั่นต่างๆในเชิงลึกพอสมควร คู่มือต้องอ่านจาก website ซึ่งก็พอเข้าใจถึงแม้จะไม่ละเอียดมากนักก็ตาม

ปัญหาการแปลงโค๊ดจาก c/c++ เป็น lazarus

  • เรื่องปัญหาการแปลงโค๊ดผมเคยเขียนไปบ้างแล้ว แต่ที่เสียเวลามากที่สุดก็คือเรื่องธรรมเนียมการเรียกใช้ฟังก์ชั่น stdcall กับ cdecl ตอนแรกไม่ได้ระวังผมใช้ stdcall ทุกฟังก์ชั่น แต่ปัญหาก็คือตอนเรียกใช้ฟังก์ชั่นใน lazarus พบว่าเกิด error ลองดูโค๊ดของ c
OGRSpatialReferenceH CPL_DLL CPL_STDCALL
OSRNewSpatialReference( const char * /* = NULL */);
OGRSpatialReferenceH CPL_DLL CPL_STDCALL OSRCloneGeogCS( OGRSpatialReferenceH );
OGRSpatialReferenceH CPL_DLL CPL_STDCALL OSRClone( OGRSpatialReferenceH );
void CPL_DLL CPL_STDCALL OSRDestroySpatialReference( OGRSpatialReferenceH );

int CPL_DLL OSRReference( OGRSpatialReferenceH );
int CPL_DLL OSRDereference( OGRSpatialReferenceH );
void CPL_DLL OSRRelease( OGRSpatialReferenceH );

OGRErr CPL_DLL OSRValidate( OGRSpatialReferenceH );
OGRErr CPL_DLL OSRFixupOrdering( OGRSpatialReferenceH );
OGRErr CPL_DLL OSRFixup( OGRSpatialReferenceH );
OGRErr CPL_DLL OSRStripCTParms( OGRSpatialReferenceH );
  • จะเห็นคีย์เวิร์ดของภาษาซีคือ CPL_STDCALL ถ้าฟังก์ชั่นมีคีย์เวิร์ดนี้ใน lazarus ต้องมีคีย์ stdcall ข้างท้ายฟังก์ชั่น ถ้าไม่มีคีย์ตัวนี้ต้องเขียน cdecl แทน ดูโค๊ดของ lazarus
function OSRNewSpatialReference(const WKT : PCHAR
) : TOGRSpatialReferenceH;
stdcall; external External_Lib name ‘_OSRNewSpatialReference@4’;
function OSRCloneGeogCS(const Handle : TOGRSpatialReferenceH
) : TOGRSpatialReferenceH;
stdcall; external External_Lib name ‘_OSRCloneGeogCS@4’;

function OSRClone(const Handle : TOGRSpatialReferenceH ) : TOGRSpatialReferenceH ;
stdcall; external External_Lib name ‘_OSRClone@4’;
procedure OSRDestroySpatialReference (const Handle : TOGRSpatialReferenceH); stdcall; external External_Lib name ‘_OSRDestroySpatialReference@4’;

function OSRReference(hSRS : TOGRSpatialReferenceH) : longint; cdecl; external External_Lib name ‘OSRReference’;
function OSRDereference(hSRS : TOGRSpatialReferenceH) : longint; cdecl; external External_Lib name ‘OSRDereference’;
procedure OSRRelease(hSRS : TOGRSpatialReferenceH ); cdecl; external External_Lib name ‘OSRRelease’;
function OSRValidate(hSRS : TOGRSpatialReferenceH ) : TOGRErr; cdecl; external External_Lib name ‘OSRValidate’;
function OSRFixupOrdering(hSRS : TOGRSpatialReferenceH) : TOGRErr; cdecl; external External_Lib name ‘OSRFixupOrdering’;
function OSRFixup(hSRS : TOGRSpatialReferenceH ) : TOGRErr; cdecl; external External_Lib name ‘OSRFixup’;
function OSRStripCTParms(hSRS : TOGRSpatialReferenceH) : TOGRErr; cdecl; external External_Lib name ‘OSRStripCTParms’;
  • ความจริงแล้ว stdcall ก็คือการเรียกฟังก์ชั่นตามแบบ Win32 นั่นเอง โธ่เอ๊ย ผมอ่านตำราของ pascal นานมาแล้วก็พูดถึงเรียกแบบนี้ มันนานจนลืม การใช้ cdecl เป็นการเรียกใช้ฟังก์ชั่นตามรูปแบบภาษาซี ซึ่งทั้งสองวิธีมีการจัดการ stack ของอาร์กิวเมนต์ของฟังก์ชั่นแตกต่างกัน cdecl ยอมให้จำนวนอาร์กิวเมนต์แปรเปลี่ยนได้เช่นฟังก์ชั่นมี 4 อาร์กิวเมนต์ แต่การเรียกใช้ฟังก์ชั่นอาจจะส่งไปแค่ 3 ตัว แต่ใน pascal จะไม่ยอม
  • สนใจเรื่องนี้อ่านวีกิได้ที่ http://en.wikibooks.org/wiki/X86_Disassembly/Calling_Conventions

ความต่างของไลบรารี GDAL/OGR ระหว่างสองโลก Linux และ Windows

  • เมื่อเข้าใจเรื่องการเรียกฟังก์ชั่น (Calling convention) ระหว่าง cdecl และ stdcall ผมก็เข้าใจได้ทันทีว่าทำไมไลบรารี GDAL ใน linux ถึงไม่เวิร์คในตอนแรกเมื่อคอมไพล์ด้วย lazarus ใน linux เองผมใช้ subversion ทำการติดตั้งและคอมไพล์ GDAL จะได้ไฟล์ไลบรารีที่เป็น .so อยู่หลายไฟล์ ตัวที่สำคัญที่สุดตอนนี้เป็น libgdal.so.1.14.0 ตัวเลขตามหลังคือรุ่นหรือเวอร์ชั่นนั่นเอง
  • และทีสำคัญที่สุดใน linux คือจะต้องเรียกฟังก์ชั่นทุกตัวด้วยคีย์ cdecl เท่านั้น ส่วนในวินโดส์ต้องดูในคีย์ภาษาซีว่าเป็น stdcall หรือ cdecl ตามที่ผมแสดงโค๊ดตัวอย่างให้ดูข้างต้น ตอนนี้ในโค๊ดผมยังไม่ได้ใช้ directive เพื่อตรวจสอบว่ากำลังอยู่ในวินโดส์หรือลีนุกซ์ โค๊ดจึงต้องแยกกันก่อน

Download sourcecode

  • สนใจก็ดาวน์โหลดโค๊ด lazarus สำหรับวินโดส์ได้ที่ gdalutm2latlong_win32.zip
  • หรือถ้าเป็น linux ก็ดาวน์โหลดได้ที่ gdalutm2latlong_linux.zip

Sourcecode ของ GDAL/OGR

  • มีไฟล์หลักๆอยู่แค่ 4 ไฟล์คือ gdal.pas, gdalcore.pas, ogr.pas และ ogrcore.pas ไฟล์ gdal.pas และ ogr.pas ผม declare ตัวแปรต่างๆที่แปลงมาจากภาษาซี ส่วน gdalcore.pas และ ogrcore.pas เป็นไฟล์ที่รวมฟังก์ชั่นที่เป็น export function จาก GDAL ไฟล์ที่จะแสดงให้ดูโค๊ดด้านล่างจะเป็นเวอร์ชั่นสำหรับวินโดส์
  • มาลองดูโค๊ดของ gdal.pas
//
//*****************************************************************************
// Build 0135
//
// Project: GDAL Lazarus Bindings
// Purpose: Types and Constants for GDAL
// Author: Frank Warmerdam, warmerdam@pobox.com
// Lazarus Wrapper : Prajuab Riabroy
//*****************************************************************************
// Copyright (c) 2005, Frank Warmerdam <warmerdam@pobox.com>
//
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//*****************************************************************************

unit gdal;

{$mode objfpc}{$H+}

interface

uses
Classes, SysUtils;

const
//External_lib = ‘gdal’;
External_lib = ‘gdal16.dll’;

type
TGInt32 = longint;
TGUInt32 = longword;
TGInt16 = shortInt;
TGUInt16 = word;
TGByte = byte;
TGBool = boolean;
TGFloat32 = single;
TGFloat64 = double;

PGInt32 = ^TGInt32;
PGUInt32 = ^TGUInt32;
PGInt16 = ^TGInt16;
PGUInt16 = ^TGUInt16;
PGByte = ^TGByte;
PGBool = ^TGBool;
PGFloat32 = ^TGFloat32;
PGFloat64 = ^TGFloat64;

PPINTEGER = ^PINTEGER;

TDynArrayDouble = array of double;
TDynArrayInt16 = array of shortint;
TDoubleArray6 = array[0..5] of double;
TDoubleArray2 = array[0..1] of double;
TDoubleArray7 = array[0..6] of double;

//* ——————————————————————– */
//* Significant constants. */
//* ——————————————————————– */
TCPLErr = (CE_None = 0, CE_Debug = 1, CE_Warning = 2, CE_Failure = 3,
CE_Fatal = 4);
//*! Pixel data types */
TGDALDataType = (GDT_Unknown = 0, GDT_Byte = 1, GDT_UInt16 = 2, GDT_Int16 = 3,
GDT_UInt32 = 4, GDT_Int32 = 5, GDT_Float32 = 6, GDT_Float64 = 7,
GDT_CInt16 = 8, GDT_CInt32 = 9, GDT_CFloat32 = 10, GDT_CFloat64 = 11,
GDT_TypeCount = 12);
//*! Flag indicating read/write, or read-only access to data. */
TGDALAccess = (GA_ReadOnly = 0, GA_Update = 1);
//*! Read/Write flag for RasterIO() method */
TGDALRWFlag = (GF_Read = 0, GF_Write = 1);
//*! Types of color interpretation for raster bands. */
TGDALColorInterp = (GCI_Undefined = 0, GCI_GrayIndex = 1, GCI_PaletteIndex = 2, GCI_RedBand = 3,
GCI_GreenBand = 4, GCI_BlueBand = 5, GCI_AlphaBand = 6, GCI_HueBand = 7,
GCI_SaturationBand = 8, GCI_LightnessBand = 9, GCI_CyanBand = 10, GCI_MagentaBand = 11,
GCI_YellowBand = 12, GCI_BlackBand = 13, GCI_YCbCr_YBand = 14, GCI_YCbCr_CbBand = 15,
GCI_YCbCr_CrBand = 16, GCI_Max = 16);
//*! Types of color interpretations for a GDALColorTable. */
TGDALPaletteInterp = (GPI_Gray = 0, GPI_RGB = 1, GPI_CMYK = 2, GPI_HLS = 3);
TGDALRATFieldType = (GFT_Integer, GFT_Real, GFT_String);
TGDALRATFieldUsage = (GFU_Generic = 0, GFU_PixelCount = 1, GFU_Name = 2, GFU_Min = 3,
GFU_Max = 4, GFU_MinMax = 5, GFU_Red = 6, GFU_Green = 7,
GFU_Blue = 8, GFU_Alpha = 9, GFU_RedMin = 10, GFU_GreenMin = 11,
GFU_BlueMin = 12, GFU_AlphaMin = 13, GFU_RedMax = 14, GFU_GreenMax = 15,
GFU_BlueMax = 16, GFU_AlphaMax = 17, GFU_MaxCount);

//* ——————————————————————– */
//* Callback "progress" function. */
//* ——————————————————————– */
TGDALProgressFunc = function(dfComplete : double; const pszMessage : PCHAR;
pProgressArg : POINTER) : longint;
TGDALDerivedPixelFunc = function(papoSources : PPOINTER; nSources : longint;
pData : POINTER; nBufXSize, nBufYSize : longint;
eSrcType : TGDALDataType;
eBufType : TGDALDataType;
nPixelSpace, nLineSpace : longint): TCPLErr;

//* ——————————————————————– */
//* Define handle types related to various internal classes. */
//* ——————————————————————– */
TGDALMajorObjectH = POINTER;
TGDALDatasetH = POINTER;
TGDALRasterBandH = POINTER;
TGDALDriverH = POINTER;
TGDALColorTableH = POINTER;
TGDALRasterAttributeTableH = POINTER;

PGDALRasterBandH = ^TGDALRasterBandH;
PGDALDatasetH = ^TGDALDatasetH;
PPGDALDatasetH = ^PGDALDatasetH;

const
//* ==================================================================== */
//* Registration/driver related. */
//* ==================================================================== */
DMD_SHORTNAME : PCHAR = ‘DMD_SHORTNAME’;
DMD_LONGNAME : PCHAR = ‘DMD_LONGNAME’;
DMD_HELPTOPIC : PCHAR = ‘DMD_HELPTOPIC’;
DMD_MIMETYPE : PCHAR = ‘DMD_MIMETYPE’;
DMD_EXTENSION : PCHAR = ‘DMD_EXTENSION’;
DMD_CREATIONOPTIONLIST : PCHAR = ‘DMD_CREATIONOPTIONLIST’;
DMD_CREATIONDATATYPES : PCHAR = ‘DMD_CREATIONDATATYPES’;

DCAP_CREATE : PCHAR = ‘DCAP_CREATE’;
DCAP_CREATECOPY : PCHAR = ‘DCAP_CREATECOPY’;

GDALMD_AREA_OR_POINT = ‘AREA_OR_POINT’;
GDALMD_AOP_AREA = ‘Area’;
GDALMD_AOP_POINT = ‘Point’;
CPLE_WrongFormat = 200;
GDAL_DMD_LONGNAME = ‘DMD_LONGNAME’;
GDAL_DMD_HELPTOPIC = ‘DMD_HELPTOPIC’;
GDAL_DMD_MIMETYPE = ‘DMD_MIMETYPE’;
GDAL_DMD_EXTENSION = ‘DMD_EXTENSION’;
GDAL_DMD_CREATIONOPTIONLIST = ‘DMD_CREATIONOPTIONLIST’;
GDAL_DMD_CREATIONDATATYPES = ‘DMD_CREATIONDATATYPES’;
GDAL_DCAP_CREATE = ‘DCAP_CREATE’;
GDAL_DCAP_CREATECOPY = ‘DCAP_CREATECOPY’;
GDAL_DCAP_VIRTUALIO = ‘DCAP_VIRTUALIO’;
//SRCVAL(papoSource, eSrcType, ii)
//SRCVAL – Macro which may be used by pixel functions to obtain a pixel from a source buffer.
GMF_ALL_VALID = $01;
GMF_PER_DATASET = $02;
GMF_ALPHA = $04;
GMF_NODATA = $08;

type
    TGDAL_GCP = record
    pszId : PCHAR;//Unique identifier, often numeric.
    pszInfo : PCHAR;//Informational message or "".
    dfGCPPixel : double;//Pixel (x) location of GCP on raster.
    dfGCPLine : double;//Line (y) location of GCP on raster.
    dfGCPX : double;//X position of GCP in georeferenced space.
    dfGCPY : double;//Y position of GCP in georeferenced space.
    dfGCPZ : double;//Elevation of GCP, or zero if not known.
end;
PGDAL_GCP = ^TGDAL_GCP;

TGDALColorEntry = record
    c1 : shortint;
    c2 : shortint;
    c3 : shortint;
    c4 : shortint;
end;
PGDALColorEntry = ^TGDALColorEntry;

TGDALRPCInfo = record
    dfLINE_OFF : double;
    dfSAMP_OFF : double;
    dfLAT_OFF : double;
    dfLONG_OFF : double;
    dfHEIGHT_OFF : double;
    dfLINE_SCALE : double;
    dfSAMP_SCALE : double;
    dfLAT_SCALE : double;
    dfLONG_SCALE : double;
    dfHEIGHT_SCALE : double;
    adfLINE_NUM_COEFF : array[0..20] of double;
    adfLINE_DEN_COEFF : array[0..20] of double;
    adfSAMP_NUM_COEFF : array[0..20] of double;
    adfSAMP_DEN_COEFF : array[0..20] of double;
    dfMIN_LONG : double;
    dfMIN_LAT : double;
    dfMAX_LONG : double;
    dfMAX_LAT : double;
end;
PGDALRPCInfo = ^TGDALRPCInfo;

implementation

end.
  • ต่อไปเป็น gdal.pas ด้านท้ายสุดเป็นคลาส TGDALMajorobject ที่ได้จากการ wrapper
//
//*****************************************************************************
// Build 0135
//
// Project:  GDAL Lazarus Bindings
// Purpose:  Types and Constants for GDAL
// Author:   Frank Warmerdam, warmerdam@pobox.com
// Lazarus Wrapper : Prajuab Riabroy
//*****************************************************************************
// Copyright (c) 2005, Frank Warmerdam <warmerdam@pobox.com>
//
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//*****************************************************************************

unit gdal;

{$mode objfpc}{$H+}

interface

uses
  Classes, SysUtils;

const
{$IFDEF MSWINDOWS}
  GDALLib = 'gdal16.dll';
{$ENDIF}
{$IFDEF LINUX}
  GDALLib = 'libgdal.so';//.1.14.0';
{$ENDIF}



type
  TGInt32   = longint;
  TGUInt32  = longword;
  TGInt16   = shortInt;
  TGUInt16  = word;
  TGByte    = byte;
  TGBool    = boolean;
  TGFloat32 = single;
  TGFloat64 = double;

  PGInt32   = ^TGInt32;
  PGUInt32  = ^TGUInt32;
  PGInt16   = ^TGInt16;
  PGUInt16  = ^TGUInt16;
  PGByte    = ^TGByte;
  PGBool    = ^TGBool;
  PGFloat32 = ^TGFloat32;
  PGFloat64 = ^TGFloat64;

  PPINTEGER = ^PINTEGER;

  TDynArrayDouble = array of double;
  TDynArrayInt16  = array of shortint;
  TDoubleArray6 = array[0..5] of double;
  TDoubleArray2 = array[0..1] of double;
  TDoubleArray7 = array[0..6] of double;

  //* -------------------------------------------------------------------- */
  //*      Significant constants.                                          */
  //* -------------------------------------------------------------------- */
  TCPLErr   =  (CE_None = 0, CE_Debug = 1, CE_Warning = 2, CE_Failure = 3,
                       CE_Fatal = 4);
  //*! Pixel data types */
  TGDALDataType = (GDT_Unknown = 0, GDT_Byte = 1, GDT_UInt16 = 2, GDT_Int16 = 3,
                   GDT_UInt32 = 4, GDT_Int32 = 5, GDT_Float32 = 6, GDT_Float64 = 7,
                   GDT_CInt16 = 8, GDT_CInt32 = 9, GDT_CFloat32 = 10, GDT_CFloat64 = 11,
                   GDT_TypeCount = 12);
  //*! Flag indicating read/write, or read-only access to data. */
  TGDALAccess = (GA_ReadOnly = 0, GA_Update = 1);
  //*! Read/Write flag for RasterIO() method */
  TGDALRWFlag = (GF_Read = 0, GF_Write = 1);
  //*! Types of color interpretation for raster bands. */
  TGDALColorInterp = (GCI_Undefined = 0, GCI_GrayIndex = 1, GCI_PaletteIndex = 2, GCI_RedBand = 3,
                      GCI_GreenBand = 4, GCI_BlueBand = 5, GCI_AlphaBand = 6, GCI_HueBand = 7,
                      GCI_SaturationBand = 8, GCI_LightnessBand = 9, GCI_CyanBand = 10, GCI_MagentaBand = 11,
                      GCI_YellowBand = 12, GCI_BlackBand = 13, GCI_YCbCr_YBand = 14, GCI_YCbCr_CbBand = 15,
                      GCI_YCbCr_CrBand = 16, GCI_Max = 16);
  //*! Types of color interpretations for a GDALColorTable. */
  TGDALPaletteInterp = (GPI_Gray = 0, GPI_RGB = 1, GPI_CMYK = 2, GPI_HLS = 3);
  TGDALRATFieldType = (GFT_Integer, GFT_Real, GFT_String);
  TGDALRATFieldUsage = (GFU_Generic = 0, GFU_PixelCount = 1, GFU_Name = 2, GFU_Min = 3,
                        GFU_Max = 4, GFU_MinMax = 5, GFU_Red = 6, GFU_Green = 7,
                        GFU_Blue = 8, GFU_Alpha = 9, GFU_RedMin = 10, GFU_GreenMin = 11,
                        GFU_BlueMin = 12, GFU_AlphaMin = 13, GFU_RedMax = 14, GFU_GreenMax = 15,
                        GFU_BlueMax = 16, GFU_AlphaMax = 17, GFU_MaxCount);

  //* -------------------------------------------------------------------- */
  //*      Callback "progress" function.                                   */
  //* -------------------------------------------------------------------- */
  TGDALProgressFunc = function(dfComplete : double; const pszMessage : PCHAR;
                              pProgressArg : POINTER) : longint;
  TGDALDerivedPixelFunc = function(papoSources : PPOINTER; nSources : longint;
                                  pData : POINTER; nBufXSize, nBufYSize : longint;
                                  eSrcType : TGDALDataType;
                                  eBufType : TGDALDataType;
                                  nPixelSpace, nLineSpace : longint): TCPLErr;

  //* -------------------------------------------------------------------- */
  //*      Define handle types related to various internal classes.        */
  //* -------------------------------------------------------------------- */
  TGDALMajorObjectH = POINTER;
  TGDALDatasetH = POINTER;
  TGDALRasterBandH = POINTER;
  TGDALDriverH = POINTER;
  TGDALColorTableH = POINTER;
  TGDALRasterAttributeTableH = POINTER;

  PGDALRasterBandH = ^TGDALRasterBandH;
  PGDALDatasetH = ^TGDALDatasetH;
  PPGDALDatasetH = ^PGDALDatasetH;

const
//* ==================================================================== */
//*      Registration/driver related.                                    */
//* ==================================================================== */
  DMD_SHORTNAME : PCHAR = 'DMD_SHORTNAME';
  DMD_LONGNAME : PCHAR = 'DMD_LONGNAME';
  DMD_HELPTOPIC : PCHAR = 'DMD_HELPTOPIC';
  DMD_MIMETYPE : PCHAR = 'DMD_MIMETYPE';
  DMD_EXTENSION : PCHAR = 'DMD_EXTENSION';
  DMD_CREATIONOPTIONLIST : PCHAR = 'DMD_CREATIONOPTIONLIST';
  DMD_CREATIONDATATYPES : PCHAR = 'DMD_CREATIONDATATYPES';

  DCAP_CREATE : PCHAR = 'DCAP_CREATE';
  DCAP_CREATECOPY : PCHAR = 'DCAP_CREATECOPY';

  GDALMD_AREA_OR_POINT = 'AREA_OR_POINT';
  GDALMD_AOP_AREA = 'Area';
  GDALMD_AOP_POINT = 'Point';
  CPLE_WrongFormat  = 200;
  GDAL_DMD_LONGNAME = 'DMD_LONGNAME';
  GDAL_DMD_HELPTOPIC = 'DMD_HELPTOPIC';
  GDAL_DMD_MIMETYPE = 'DMD_MIMETYPE';
  GDAL_DMD_EXTENSION = 'DMD_EXTENSION';
  GDAL_DMD_CREATIONOPTIONLIST = 'DMD_CREATIONOPTIONLIST';
  GDAL_DMD_CREATIONDATATYPES = 'DMD_CREATIONDATATYPES';
  GDAL_DCAP_CREATE = 'DCAP_CREATE';
  GDAL_DCAP_CREATECOPY = 'DCAP_CREATECOPY';
  GDAL_DCAP_VIRTUALIO = 'DCAP_VIRTUALIO';
  //SRCVAL(papoSource, eSrcType, ii)
  //SRCVAL - Macro which may be used by pixel functions to obtain a pixel from a source buffer.
  GMF_ALL_VALID  = $01;
  GMF_PER_DATASET =  $02;
  GMF_ALPHA  = $04;
  GMF_NODATA  = $08;

type
  TGDAL_GCP = record
    pszId : PCHAR;//Unique identifier, often numeric.
    pszInfo : PCHAR;//Informational message or "".
    dfGCPPixel : double;//Pixel (x) location of GCP on raster.
    dfGCPLine : double;//Line (y) location of GCP on raster.
    dfGCPX : double;//X position of GCP in georeferenced space.
    dfGCPY : double;//Y position of GCP in georeferenced space.
    dfGCPZ : double;//Elevation of GCP, or zero if not known.
  end;
  PGDAL_GCP = ^TGDAL_GCP;

  TGDALColorEntry = record
    c1 : shortint;
    c2 : shortint;
    c3 : shortint;
    c4 : shortint;
  end;
  PGDALColorEntry = ^TGDALColorEntry;

  TGDALRPCInfo = record
    dfLINE_OFF : double;
    dfSAMP_OFF : double;
    dfLAT_OFF  : double;
    dfLONG_OFF : double;
    dfHEIGHT_OFF : double;
    dfLINE_SCALE : double;
    dfSAMP_SCALE : double;
    dfLAT_SCALE  : double;
    dfLONG_SCALE : double;
    dfHEIGHT_SCALE : double;
    adfLINE_NUM_COEFF : array[0..20] of double;
    adfLINE_DEN_COEFF : array[0..20] of double;
    adfSAMP_NUM_COEFF : array[0..20] of double;
    adfSAMP_DEN_COEFF : array[0..20] of double;
    dfMIN_LONG : double;
    dfMIN_LAT  : double;
    dfMAX_LONG : double;
    dfMAX_LAT  : double;
  end;
  PGDALRPCInfo = ^TGDALRPCInfo;

implementation

end.
  • ต่อไปมาดูโค๊ด ogr.pas
//
//'*****************************************************************************
//' Build 0135
//'
//' Project:  GDAL Lazarus Bindings
//' Purpose:  Types and Constants for OGR
//' Author:   Frank Warmerdam, warmerdam@pobox.com
//' Lazarus Wrapper : Prajuab Riabroy
//'*****************************************************************************
//' Copyright (c) 2005, Frank Warmerdam <warmerdam@pobox.com>
//'
//' Permission is hereby granted, free of charge, to any person obtaining a
//' copy of this software and associated documentation files (the "Software"),
//' to deal in the Software without restriction, including without limitation
//' the rights to use, copy, modify, merge, publish, distribute, sublicense,
//' and/or sell copies of the Software, and to permit persons to whom the
//' Software is furnished to do so, subject to the following conditions:
//'
//' The above copyright notice and this permission notice shall be included
//' in all copies or substantial portions of the Software.
//'
//' THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
//' OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
//' FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
//' THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
//' LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
//' FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
//' DEALINGS IN THE SOFTWARE.
//'*****************************************************************************

unit ogr;

{$mode objfpc}{$H+}

interface

uses
  Classes, SysUtils, gdal;

type
  TOGRSpatialReferenceH = POINTER;
  TOGRCoordinateTransformationH = POINTER;

  TOGREnvelope = record
    MinX : double;
    MaxX : double;
    MinY : double;
    MaxY : double;
  end;

  TOGRErr = (OGRERR_NONE = 0,
             OGRERR_NOT_ENOUGH_DATA = 1,    //* not enough data to deserialize */
             OGRERR_NOT_ENOUGH_MEMORY = 2,
             OGRERR_UNSUPPORTED_GEOMETRY_TYPE = 3,
             OGRERR_UNSUPPORTED_OPERATION = 4,
             OGRERR_CORRUPT_DATA = 5,
             OGRERR_FAILURE = 6,
             OGRERR_UNSUPPORTED_SRS = 7,
             OGRERR_INVALID_HANDLE = 8);

  POGRErr = ^ TOGRErr;

  TOGRwkbGeometryType = (wkbUnknown = 0,
                         wkbPoint = 1,
                         wkbLineString = 2,
                         wkbPolygon = 3,
                         wkbMultiPoint = 4,
                         wkbMultiLineString = 5,
                         wkbMultiPolygon = 6,
                         wkbGeometryCollection = 7,
                         wkbNone = 100,
                         wkbLinearRing = 101,
                         wkbPoint25D = $80000001,
                         wkbLineString25D = $80000002,
                         wkbPolygon25D = $80000003,
                         wkbMultiPoint25D = $80000004,
                         wkbMultiLineString25D = $80000005,
                         wkbMultiPolygon25D = $80000006,
                         wkbGeometryCollection25D = $80000007) ;
const
  wkb25DBit = $80000000;
  ogrZMarker = $21125711;
type
  TOGRwkbByteOrder = (wkbXDR = 0,     //* MSB/Sun/Motoroloa: Most Significant Byte First   */
                      wkbNDR = 1);    //* LSB/Intel/Vax: Least Significant Byte First      */
  TOGRFieldType = (OFTInteger = 0,
                   OFTIntegerList = 1,
                   OFTReal = 2,
                   OFTRealList = 3,
                   OFTString = 4,
                   OFTStringList = 5,
                   OFTWideString = 6,
                   OFTWideStringList = 7,
                   OFTBinary = 8,
                   OFTDate = 9,
                   OFTTime = 10,
                   OFTDateTime = 11,
                   OFTMaxType = 11);
  TOGRJustification = (OJUndefined = 0,
                       OJLeft = 1,
                       OJRight = 2);
const
  OGRNullFID = -1;
  OGRUnsetMarker = -21121;

  function OGRGeometryTypeToName(eType : TOGRwkbGeometryType) : PCHAR; stdcall; external GDALLib name 'OGRGeometryTypeToName';
  function OGRMergeGeometryTypes(eMain, eExtra : TOGRwkbGeometryType) : TOGRwkbGeometryType; stdcall; external GDALLib name 'OGRMergeGeometryTypes';

type
  POGRField = ^TOGRField;
  TOGRField = record
    case longint of
      0 : ( Integer : longint );
      1 : ( Real : double);
      2 : ( _String : Pchar );
      3 : ( IntegerList : record
            nCount : longint;
            paList : Plongint;
           end );
      4 : ( TRealList : record
              nCount : longint;
              paList : PDOUBLE;
           end );
      5 : ( TStringList : record
            nCount : longint;
            paList : PPCHAR;
           end );
      6 : ( TBinary : record
            nCount : longint;
            paData : PGByte;
          end );
      7 : ( TSet : record
            nMarker1 : longint;
            nMarker2 : longint;
           end );
      8 : ( TDate : record
            Year : TGInt16;
            Month : TGByte;
            Day : TGByte;
            Hour : TGByte;
            Minute : TGByte;
            Second : TGByte;
            TZFlag : TGByte;
          end );
  end;

  function OGRParseDate(pszInput : PCHAR; psOutput : POGRField;
                        nOptions : longint) : longint; stdcall; external GDALLib name 'OGRParseDate';
  //* -------------------------------------------------------------------- */
  //*      Constants from ogrsf_frmts.h for capabilities.                  */
  //* -------------------------------------------------------------------- */
const
  OLCRandomRead         = 'RandomRead';
  OLCSequentialWrite    = 'SequentialWrite';
  OLCRandomWrite        = 'RandomWrite';
  OLCFastSpatialFilter  = 'FastSpatialFilter';
  OLCFastFeatureCount   = 'FastFeatureCount';
  OLCFastGetExtent      = 'FastGetExtent';
  OLCCreateField        = 'CreateField';
  OLCTransactions       = 'Transactions';
  OLCDeleteFeature      = 'DeleteFeature';
  OLCFastSetNextByIndex = 'FastSetNextByIndex';
  OLCStringsAsUTF8      = 'StringsAsUTF8';
  ODsCCreateLayer       = 'CreateLayer';
  ODsCDeleteLayer       = 'DeleteLayer';
  ODrCCreateDataSource  = 'CreateDataSource';
  ODrCDeleteDataSource  = 'DeleteDataSource';
type
  Togr_style_tool_class_id = (OGRSTCNone = 0,
                              OGRSTCPen = 1,
                              OGRSTCBrush = 2,
                              OGRSTCSymbol = 3,
                              OGRSTCLabel = 4,
                              OGRSTCVector = 5);

  Togr_style_tool_param_pen_id = (OGRSTPenColor = 0,
                                  OGRSTPenWidth = 1,
                                  OGRSTPenPattern = 2,
                                  OGRSTPenId = 3,
                                  OGRSTPenPerOffset = 4,
                                  OGRSTPenCap = 5,
                                  OGRSTPenJoin = 6,
                                  OGRSTPenPriority = 7,
                                  OGRSTPenLast = 8);

  Togr_style_tool_param_symbol_id = (OGRSTSymbolId = 0,
                                     OGRSTSymbolAngle = 1,
                                     OGRSTSymbolColor = 2,
                                     OGRSTSymbolSize = 3,
                                     OGRSTSymbolDx = 4,
                                     OGRSTSymbolDy = 5,
                                     OGRSTSymbolStep = 6,
                                     OGRSTSymbolPerp = 7,
                                     OGRSTSymbolOffset = 8,
                                     OGRSTSymbolPriority = 9,
                                     OGRSTSymbolFontName = 10,
                                     OGRSTSymbolOColor = 11,
                                     OGRSTSymbolLast = 12);

  Togr_style_tool_units_id = (OGRSTUGround = 0,
                              OGRSTUPixel = 1,
                              OGRSTUPoints = 2,
                              OGRSTUMM = 3,
                              OGRSTUCM = 4,
                              OGRSTUInches = 5);
  Togr_style_tool_param_brush_id = (OGRSTBrushFColor = 0,
                                    OGRSTBrushBColor = 1,
                                    OGRSTBrushId = 2,
                                    OGRSTBrushAngle = 3,
                                    OGRSTBrushSize = 4,
                                    OGRSTBrushDx = 5,
                                    OGRSTBrushDy = 6,
                                    OGRSTBrushPriority = 7,
                                    OGRSTBrushLast = 8);
  Togr_style_tool_param_label_id  = (OGRSTLabelFontName = 0,
                                     OGRSTLabelSize = 1,
                                     OGRSTLabelTextString = 2,
                                     OGRSTLabelAngle = 3,
                                     OGRSTLabelFColor = 4,
                                     OGRSTLabelBColor = 5,
                                     OGRSTLabelPlacement = 6,
                                     OGRSTLabelAnchor = 7,
                                     OGRSTLabelDx = 8,
                                     OGRSTLabelDy = 9,
                                     OGRSTLabelPerp = 10,
                                     OGRSTLabelBold = 11,
                                     OGRSTLabelItalic = 12,
                                     OGRSTLabelUnderline = 13,
                                     OGRSTLabelPriority = 14,
                                     OGRSTLabelStrikeout = 15,
                                     OGRSTLabelStretch = 16,
                                     OGRSTLabelAdjHor = 17,
                                     OGRSTLabelAdjVert = 18,
                                     OGRSTLabelHColor = 19,
                                     OGRSTLabelOColor = 20,
                                     OGRSTLabelLast = 21);


  TOGRSTUnitId = Togr_style_tool_units_id;
  TOGRSTPenParam = Togr_style_tool_param_pen_id;
  TOGRSTBrushParam = Togr_style_tool_param_brush_id;
  TOGRSTSymbolParam = Togr_style_tool_param_symbol_id;
  TOGRSTClassId = Togr_style_tool_class_id;
  TOGRSTLabelParam = Togr_style_tool_param_label_id;

  Pogr_style_tool_class_id  = ^Togr_style_tool_class_id;
  Pogr_style_tool_param_brush_id  = ^Togr_style_tool_param_brush_id;
  Pogr_style_tool_param_label_id  = ^Togr_style_tool_param_label_id;
  Pogr_style_tool_param_pen_id  = ^Togr_style_tool_param_pen_id;
  Pogr_style_tool_param_symbol_id  = ^Togr_style_tool_param_symbol_id;
  Pogr_style_tool_units_id  = ^Togr_style_tool_units_id;
  POGRSTBrushParam  = ^TOGRSTBrushParam;
  POGRSTClassId  = ^TOGRSTClassId;
  POGRSTLabelParam  = ^TOGRSTLabelParam;
  POGRSTPenParam  = ^TOGRSTPenParam;
  POGRSTSymbolParam  = ^TOGRSTSymbolParam;
  POGRSTUnitId  = ^TOGRSTUnitId;

  TOGRAxisOrientation = (OAO_Other = 0,
                         OAO_North = 1,
                         OAO_South = 2,
                         OAO_East = 3,
                         OAO_West = 4,
                         OAO_Up = 5,
                         OAO_Down = 6) ;

  function OSRAxisEnumToName(eOrientation : TOGRAxisOrientation
                        ) : PCHAR; cdecl; external GDALLib name 'OSRAxisEnumToName';

  //* -------------------------------------------------------------------- */
  //*      Datum types (corresponds to CS_DatumType).                      */
  //* -------------------------------------------------------------------- */
type
  TOGRDatumType = (
      ODT_HD_Min = 1000,
      ODT_HD_Other = 1000,
      ODT_HD_Classic = 1001,
      ODT_HD_Geocentric = 1002,
      ODT_HD_Max = 1999,
      ODT_VD_Min = 2000,
      ODT_VD_Other = 2000,
      ODT_VD_Orthometric = 2001,
      ODT_VD_Ellipsoidal = 2002,
      ODT_VD_AltitudeBarometric = 2003,
      ODT_VD_Normal = 2004,
      ODT_VD_GeoidModelDerived = 2005,
      ODT_VD_Depth = 2006,
      ODT_VD_Max = 2999,
      ODT_LD_Min = 10000,
      ODT_LD_Max = 32767);


const
  //SRS_WKT_WGS84 : PCHAR = 'GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9108\"]],AUTHORITY[\"EPSG\",\"4326\"]]';

  //* ==================================================================== */
  //*      Some "standard" strings.                                        */
  //* ==================================================================== */

  SRS_PT_ALBERS_CONIC_EQUAL_AREA = 'Albers_Conic_Equal_Area';
  SRS_PT_AZIMUTHAL_EQUIDISTANT = 'Azimuthal_Equidistant';
  SRS_PT_CASSINI_SOLDNER = 'Cassini_Soldner';
  SRS_PT_CYLINDRICAL_EQUAL_AREA = 'Cylindrical_Equal_Area';
  SRS_PT_BONNE = 'Bonne';
  SRS_PT_ECKERT_I = 'Eckert_I';
  SRS_PT_ECKERT_II = 'Eckert_II';
  SRS_PT_ECKERT_III = 'Eckert_III';
  SRS_PT_ECKERT_IV = 'Eckert_IV';
  SRS_PT_ECKERT_V = 'Eckert_V';
  SRS_PT_ECKERT_VI = 'Eckert_VI';
  SRS_PT_EQUIDISTANT_CONIC = 'Equidistant_Conic';
  SRS_PT_EQUIRECTANGULAR = 'Equirectangular';
  SRS_PT_GALL_STEREOGRAPHIC = 'Gall_Stereographic';
  SRS_PT_GAUSSSCHREIBERTMERCATOR = 'Gauss_Schreiber_Transverse_Mercator';
  SRS_PT_GEOSTATIONARY_SATELLITE = 'Geostationary_Satellite';
  SRS_PT_GOODE_HOMOLOSINE = 'Goode_Homolosine';
  SRS_PT_GNOMONIC = 'Gnomonic';
  SRS_PT_HOTINE_OBLIQUE_MERCATOR = 'Hotine_Oblique_Mercator';
  SRS_PT_HOTINE_OBLIQUE_MERCATOR_TWO_POINT_NATURAL_ORIGIN=
                'Hotine_Oblique_Mercator_Two_Point_Natural_Origin';
  SRS_PT_LABORDE_OBLIQUE_MERCATOR = 'Laborde_Oblique_Mercator';
  SRS_PT_LAMBERT_CONFORMAL_CONIC_1SP = 'Lambert_Conformal_Conic_1SP';
  SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP = 'Lambert_Conformal_Conic_2SP';
  SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP_BELGIUM =
                                  'Lambert_Conformal_Conic_2SP_Belgium';
  SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA =
                                  'Lambert_Azimuthal_Equal_Area';
  SRS_PT_MERCATOR_1SP = 'Mercator_1SP';
  SRS_PT_MERCATOR_2SP = 'Mercator_2SP';
  SRS_PT_MILLER_CYLINDRICAL = 'Miller_Cylindrical';
  SRS_PT_MOLLWEIDE = 'Mollweide';
  SRS_PT_NEW_ZEALAND_MAP_GRID = 'New_Zealand_Map_Grid';
  SRS_PT_OBLIQUE_STEREOGRAPHIC = 'Oblique_Stereographic';
  SRS_PT_ORTHOGRAPHIC = 'Orthographic';
  SRS_PT_POLAR_STEREOGRAPHIC = 'Polar_Stereographic';
  SRS_PT_POLYCONIC='Polyconic';
  SRS_PT_ROBINSON = 'Robinson';
  SRS_PT_SINUSOIDAL = 'Sinusoidal';
  SRS_PT_STEREOGRAPHIC = 'Stereographic';
  SRS_PT_SWISS_OBLIQUE_CYLINDRICAL = 'Swiss_Oblique_Cylindrical';
  SRS_PT_TRANSVERSE_MERCATOR = 'Transverse_Mercator';
  SRS_PT_TRANSVERSE_MERCATOR_SOUTH_ORIENTED
                                 = 'Transverse_Mercator_South_Orientated';

  //* special mapinfo variants on Transverse Mercator */
  SRS_PT_TRANSVERSE_MERCATOR_MI_21 =
                                  'Transverse_Mercator_MapInfo_21';
  SRS_PT_TRANSVERSE_MERCATOR_MI_22 =
                                  'Transverse_Mercator_MapInfo_22';
  SRS_PT_TRANSVERSE_MERCATOR_MI_23 =
                                  'Transverse_Mercator_MapInfo_23';
  SRS_PT_TRANSVERSE_MERCATOR_MI_24 =
                                  'Transverse_Mercator_MapInfo_24';
  SRS_PT_TRANSVERSE_MERCATOR_MI_25 =
                                  'Transverse_Mercator_MapInfo_25';

  SRS_PT_TUNISIA_MINING_GRID =
                                  'Tunisia_Mining_Grid';
  SRS_PT_TWO_POINT_EQUIDISTANT =
                                  'Two_Point_Equidistant';
  SRS_PT_VANDERGRINTEN =  'VanDerGrinten';
  SRS_PT_KROVAK =          'Krovak';
  SRS_PT_IMW_POLYCONIC =   'International_Map_of_the_World_Polyconic';
  SRS_PT_WAGNER_I  =       'Wagner_I';
  SRS_PT_WAGNER_II =       'Wagner_II';
  SRS_PT_WAGNER_III =      'Wagner_III';
  SRS_PT_WAGNER_IV =       'Wagner_IV';
  SRS_PT_WAGNER_V =        'Wagner_V';
  SRS_PT_WAGNER_VI =       'Wagner_VI';
  SRS_PT_WAGNER_VII =      'Wagner_VII';


  SRS_PP_CENTRAL_MERIDIAN =        'central_meridian';
  SRS_PP_SCALE_FACTOR     =        'scale_factor';
  SRS_PP_STANDARD_PARALLEL_1 =     'standard_parallel_1';
  SRS_PP_STANDARD_PARALLEL_2 =     'standard_parallel_2';
  SRS_PP_PSEUDO_STD_PARALLEL_1 =   'pseudo_standard_parallel_1';
  SRS_PP_LONGITUDE_OF_CENTER =     'longitude_of_center';
  SRS_PP_LATITUDE_OF_CENTER =      'latitude_of_center';
  SRS_PP_LONGITUDE_OF_ORIGIN =     'longitude_of_origin';
  SRS_PP_LATITUDE_OF_ORIGIN =      'latitude_of_origin';
  SRS_PP_FALSE_EASTING  =          'false_easting';
  SRS_PP_FALSE_NORTHING =          'false_northing';
  SRS_PP_AZIMUTH =                 'azimuth';
  SRS_PP_LONGITUDE_OF_POINT_1 =    'longitude_of_point_1';
  SRS_PP_LATITUDE_OF_POINT_1  =    'latitude_of_point_1';
  SRS_PP_LONGITUDE_OF_POINT_2 =    'longitude_of_point_2';
  SRS_PP_LATITUDE_OF_POINT_2  =    'latitude_of_point_2';
  SRS_PP_LONGITUDE_OF_POINT_3 =    'longitude_of_point_3';
  SRS_PP_LATITUDE_OF_POINT_3  =    'latitude_of_point_3';
  SRS_PP_RECTIFIED_GRID_ANGLE =    'rectified_grid_angle';
  SRS_PP_LANDSAT_NUMBER =          'landsat_number';
  SRS_PP_PATH_NUMBER =             'path_number';
  SRS_PP_PERSPECTIVE_POINT_HEIGHT = 'perspective_point_height';
  SRS_PP_SATELLITE_HEIGHT =        'satellite_height';
  SRS_PP_FIPSZONE =                'fipszone';
  SRS_PP_ZONE =                    'zone';
  SRS_PP_LATITUDE_OF_1ST_POINT =   'Latitude_Of_1st_Point';
  SRS_PP_LONGITUDE_OF_1ST_POINT =  'Longitude_Of_1st_Point';
  SRS_PP_LATITUDE_OF_2ND_POINT =   'Latitude_Of_2nd_Point';
  SRS_PP_LONGITUDE_OF_2ND_POINT =  'Longitude_Of_2nd_Point';

  SRS_UL_METER =           'Meter';
  SRS_UL_FOOT  =           'Foot (International)'; //* or just 'FOOT'? */
  SRS_UL_FOOT_CONV =                   '0.3048';
  SRS_UL_US_FOOT =         'Foot_US'; //* or 'US survey foot' from EPSG */
  SRS_UL_US_FOOT_CONV =                '0.3048006096012192';
  SRS_UL_NAUTICAL_MILE =   'Nautical Mile';
  SRS_UL_NAUTICAL_MILE_CONV =          '1852.0';
  SRS_UL_LINK =            'Link';          //* Based on US Foot */
  SRS_UL_LINK_CONV =                   '0.20116684023368047';
  SRS_UL_CHAIN =           'Chain';         //* based on US Foot */
  SRS_UL_CHAIN_CONV =                  '20.116684023368047';
  SRS_UL_ROD =             'Rod';           //* based on US Foot */
  SRS_UL_ROD_CONV =                    '5.02921005842012';

  SRS_UA_DEGREE =          'degree';
  SRS_UA_DEGREE_CONV =                 '0.0174532925199433';
  SRS_UA_RADIAN =          'radian';

  SRS_PM_GREENWICH =       'Greenwich';

  SRS_DN_NAD27 =           'North_American_Datum_1927';
  SRS_DN_NAD83 =           'North_American_Datum_1983';
  SRS_DN_WGS72 =           'WGS_1972';
  SRS_DN_WGS84 =           'WGS_1984';

  SRS_WGS84_SEMIMAJOR =    6378137.0;
  SRS_WGS84_INVFLATTENING = 298.257223563;

implementation

end.
  • สุดท้ายเป็นไฟล์ ogrcore.pas
{'*****************************************************************************
' Build 0135
'
' Project:  GDAL Lazarus Bindings
' Purpose:  OGR Lazarus internal support functions.
' Author:   Frank Warmerdam, warmerdam@pobox.com
' Lazarus Wrapper : Prajuab Riabroy
'*****************************************************************************
' Copyright (c) 2005, Frank Warmerdam <warmerdam@pobox.com>
'
' Permission is hereby granted, free of charge, to any person obtaining a
' copy of this software and associated documentation files (the "Software"),
' to deal in the Software without restriction, including without limitation
' the rights to use, copy, modify, merge, publish, distribute, sublicense,
' and/or sell copies of the Software, and to permit persons to whom the
' Software is furnished to do so, subject to the following conditions:
'
' The above copyright notice and this permission notice shall be included
' in all copies or substantial portions of the Software.
'
' THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
' OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
' FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
' THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
' LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
' FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
' DEALINGS IN THE SOFTWARE.
'*****************************************************************************
}

unit ogrcore;

{$mode objfpc}{$H+}

interface

uses
  Classes, SysUtils, gdal, ogr;
  procedure OGRFree(p : POINTER); cdecl; external GDALLib name 'OGRFree';
{$IFDEF MSWINDOWS}
  function OSRNewSpatialReference(const WKT : PCHAR
           ) : TOGRSpatialReferenceH; stdcall; external GDALLib name '_OSRNewSpatialReference@4';
  function OSRCloneGeogCS(const Handle : TOGRSpatialReferenceH
           ) : TOGRSpatialReferenceH; stdcall; external GDALLib name '_OSRCloneGeogCS@4';
  function OSRClone(const Handle : TOGRSpatialReferenceH
           ) : TOGRSpatialReferenceH ; stdcall; external GDALLib name '_OSRClone@4';
  procedure OSRDestroySpatialReference (const Handle : TOGRSpatialReferenceH
            ); stdcall; external GDALLib name '_OSRDestroySpatialReference@4';
  {$ENDIF}
  {$IFDEF LINUX}
  function OSRNewSpatialReference(const WKT : PCHAR
           ) : TOGRSpatialReferenceH; cdecl; external GDALLib name 'OSRNewSpatialReference';
  function OSRCloneGeogCS(const Handle : TOGRSpatialReferenceH
           ) : TOGRSpatialReferenceH; cdecl; external GDALLib name 'OSRCloneGeogCS';
  function OSRClone(const Handle : TOGRSpatialReferenceH
           ) : TOGRSpatialReferenceH ; cdecl; external GDALLib name 'OSRClone';
  procedure OSRDestroySpatialReference (const Handle : TOGRSpatialReferenceH
            ); cdecl; external GDALLib name 'OSRDestroySpatialReference';
  {$ENDIF}
  function  OSRReference(hSRS : TOGRSpatialReferenceH
            ) : longint; cdecl; external GDALLib name 'OSRReference';
  function  OSRDereference(hSRS : TOGRSpatialReferenceH
            ) : longint; cdecl; external GDALLib name 'OSRDereference';
  procedure OSRRelease(hSRS : TOGRSpatialReferenceH
            ); cdecl; external GDALLib name 'OSRRelease';
  function  OSRValidate(hSRS : TOGRSpatialReferenceH
            ) : TOGRErr; cdecl; external GDALLib name 'OSRValidate';
  function  OSRFixupOrdering(hSRS : TOGRSpatialReferenceH
            ) : TOGRErr; cdecl; external GDALLib name 'OSRFixupOrdering';
  function  OSRFixup(hSRS : TOGRSpatialReferenceH
            ) : TOGRErr; cdecl; external GDALLib name 'OSRFixup';
  function  OSRStripCTParms(hSRS : TOGRSpatialReferenceH
            ) : TOGRErr; cdecl; external GDALLib name 'OSRStripCTParms';
  {$IFDEF MSWINDOWS}
  function  OSRImportFromEPSG(hSRS : TOGRSpatialReferenceH;
            EPSGCode : longint
            ) : TOGRErr; stdcall; external GDALLib name '_OSRImportFromEPSG@8';
  {$ENDIF}
  {$IFDEF LINUX}
  function  OSRImportFromEPSG(hSRS : TOGRSpatialReferenceH;
            EPSGCode : longint
            ) : TOGRErr; cdecl; external GDALLib name 'OSRImportFromEPSG';
  {$ENDIF}
  function  OSRImportFromWkt(hSRS : TOGRSpatialReferenceH;
            ppszInput : PPCHAR
            ) : TOGRErr; cdecl; external GDALLib name 'OSRImportFromWkt';
  function  OSRImportFromProj4(hSRS : TOGRSpatialReferenceH;
            pszProj4 : PCHAR
            ) : TOGRErr; cdecl; external GDALLib name 'OSRImportFromProj4';
  function  OSRImportFromESRI(hSRS : TOGRSpatialReferenceH;
            papszPrj : PPCHAR
            ) : TOGRErr; cdecl; external GDALLib name 'OSRImportFromESRI';
  function  OSRImportFromPCI(hSRS : TOGRSpatialReferenceH; pszProj : PCHAR;
	    pszUnits : PCHAR; padfPrjParams : PDOUBLE
            ) : TOGRErr; cdecl; external GDALLib name 'OSRImportFromPCI';
  function  OSRImportFromUSGS(hSRS : TOGRSpatialReferenceH; iProjsys : longint;
                              iZone : longint; padfPrjParams : PDOUBLE;
                              iDatum : longint
                              ) : TOGRErr; cdecl; external GDALLib name 'OSRImportFromUSGS';
  function  OSRImportFromXML(hSRS : TOGRSpatialReferenceH;
                              pszXML : PCHAR) : TOGRErr; cdecl; external GDALLib name 'OSRImportFromXML';
  function  OSRImportFromMICoordSys(hSRS : TOGRSpatialReferenceH;
                                    pszCoordSys : PCHAR
                                   ) : TOGRErr; cdecl; external GDALLib name 'OSRImportFromMICoordSys';
  function  OSRImportFromUrl(hSRS : TOGRSpatialReferenceH;
                             pszUrl : PCHAR
                            ) : TOGRErr; cdecl; external GDALLib name 'OSRImportFromUrl';

  function  OSRExportToWkt(hSRS : TOGRSpatialReferenceH;
                           ppszReturn : PPCHAR) : TOGRErr; stdcall; external GDALLib name '_OSRExportToWkt@8';
  function  OSRExportToPrettyWkt(hSRS : TOGRSpatialReferenceH;
                                 ppszReturn : PPCHAR;
                                 bSimplify: longint) : TOGRErr; stdcall; external GDALLib name '_OSRExportToPrettyWkt@12';

  function  OSRExportToProj4(hSRS : TOGRSpatialReferenceH;
                             ppszReturn : PPCHAR) : TOGRErr; stdcall; external GDALLib name '_OSRExportToProj4@8';
  function  OSRExportToPCI(hSRS : TOGRSpatialReferenceH; ppszProj : PPCHAR;
                           ppszUnits : PPCHAR;
                           ppadfPrjParams : PPDOUBLE
                           ) : TOGRErr; cdecl; external GDALLib name 'OSRExportToPCI';
  function  OSRExportToUSGS(hSRS : TOGRSpatialReferenceH; piProjSys : PINTEGER;
                            piZone : PLONGINT;
                            ppadfPrjParams : PPDOUBLE;
                            piDatum : PLONGINT
                            ) : TOGRErr; cdecl; external GDALLib name 'OSRExportToUSGS';
  function  OSRExportToXML(hSRS : TOGRSpatialReferenceH; ppszRawXML : PPCHAR;
                           pszDialect : PCHAR
                           ) : TOGRErr; cdecl; external GDALLib name 'OSRExportToXML';
  function  OSRExportToMICoordSys(hSRS : TOGRSpatialReferenceH;
            ppszReturn : PPCHAR
            ) : TOGRErr; cdecl; external GDALLib name 'OSRExportToMICoordSys';
  function  OSRMorphToESRI(hSRS : TOGRSpatialReferenceH
            ) : TOGRErr; cdecl; external GDALLib name 'OSRMorphToESRI';
  function  OSRMorphFromESRI(hSRS : TOGRSpatialReferenceH
            ) : TOGRErr; cdecl; external GDALLib name 'OSRMorphFromESRI';
  {$IFDEF MSWINDOWS}
  function  OSRSetAttrValue(hSRS : TOGRSpatialReferenceH;
            pszNodePath : PCHAR;
            pszNewNodeValue : PCHAR
            ) : TOGRErr; stdcall; external GDALLib name '_OSRSetAttrValue@12';
  function  OSRGetAttrValue(hSRS : TOGRSpatialReferenceH;
            pszName : PCHAR;
            iChild  : longint(* = 0 *)
            ) : PCHAR; stdcall; external GDALLib name '_OSRGetAttrValue@12';
  {$ENDIF}
  {$IFDEF LINUX}
  function  OSRSetAttrValue(hSRS : TOGRSpatialReferenceH;
            pszNodePath : PCHAR;
            pszNewNodeValue : PCHAR
            ) : TOGRErr; cdecl; external GDALLib name 'OSRSetAttrValue';
  function  OSRGetAttrValue(hSRS : TOGRSpatialReferenceH;
            pszName : PCHAR;
            iChild  : longint(* = 0 *)
            ) : PCHAR; cdecl; external GDALLib name 'OSRGetAttrValue';
  {$ENDIF}

  function  OSRSetAngularUnits(hSRS : TOGRSpatialReferenceH; pszUnits : PCHAR;
            dfInRadians : double
            ) : TOGRErr; cdecl; external GDALLib name 'OSRSetAngularUnits';
  function  OSRGetAngularUnits(hSRS : TOGRSpatialReferenceH;
            ppszName : PPCHAR
            ) : double; cdecl; external GDALLib name 'OSRGetAngularUnits';
  function  OSRSetLinearUnits(hSRS : TOGRSpatialReferenceH; pszUnits : PCHAR;
            dfInMeters : double
            ) : TOGRErr; cdecl; external GDALLib name 'OSRSetLinearUnits';
  function  OSRSetLinearUnitsAndUpdateParameters(hSRS : TOGRSpatialReferenceH;
            pszUnits : PCHAR;
            dfInMeters : double
            ) : TOGRErr; cdecl; external GDALLib name 'OSRSetLinearUnitsAndUpdateParameters';
  function  OSRGetLinearUnits(hSRS : TOGRSpatialReferenceH;
            ppszName : PPCHAR
            ) : double; cdecl; external GDALLib name 'OSRGetLinearUnits';
  function  OSRGetPrimeMeridian(hSRS : TOGRSpatialReferenceH;
            ppszName : PPCHAR
            ) : double; cdecl; external GDALLib name 'OSRGetPrimeMeridian';
  function  OSRIsGeographic(hSRS : TOGRSpatialReferenceH
            ) : longint; cdecl; external GDALLib name 'OSRIsGeographic';
  function  OSRIsLocal(hSRS : TOGRSpatialReferenceH
            ) : longint; cdecl; external GDALLib name 'OSRIsLocal';
  function  OSRIsProjected(hSRS : TOGRSpatialReferenceH
            ) : longint; cdecl; external GDALLib name 'OSRIsProjected';
  function  OSRIsSameGeogCS(hSRS1 : TOGRSpatialReferenceH;
            hSRS2 : TOGRSpatialReferenceH
            ) : longint; cdecl; external GDALLib name 'OSRIsSameGeogCS';
  function  OSRIsSame(hSRS1 : TOGRSpatialReferenceH;
            hSRS2 : TOGRSpatialReferenceH
            ) : longint; cdecl; external GDALLib name 'OSRIsSame';
  function  OSRSetLocalCS(hSRS : TOGRSpatialReferenceH;
            pszName : PCHAR
            ) : TOGRErr; cdecl; external GDALLib name 'OSRSetLocalCS';
  function  OSRSetProjCS(hSRS : TOGRSpatialReferenceH;
            pszName : PCHAR
            ) : TOGRErr; cdecl; external GDALLib name 'OSRSetProjCS';
  function  OSRSetWellKnownGeogCS(hSRS : TOGRSpatialReferenceH;
            pszName : PCHAR
            ) : TOGRErr; cdecl; external GDALLib name 'OSRSetWellKnownGeogCS';
  function  OSRSetFromUserInput(hSRS : TOGRSpatialReferenceH;
            pszDef : PCHAR
            ): TOGRErr; stdcall; external GDALLib name '_OSRSetFromUserInput@8';
  function  OSRCopyGeogCSFrom(hSRS : TOGRSpatialReferenceH;
            hSrcSRS : TOGRSpatialReferenceH
            ) : TOGRErr; cdecl; external GDALLib name 'OSRCopyGeogCSFrom';
  function  OSRSetTOWGS84(hSRS :  TOGRSpatialReferenceH;
            dfDX, dfDY,
            dfDZ, dfEX,
            dfEY, dfEZ,
            dfPPM : double
            ) : TOGRErr; cdecl; external GDALLib name 'OSRSetTOWGS84';
  function  OSRGetTOWGS84(hSRS : TOGRSpatialReferenceH; padfCoeff : PDOUBLE;
            nCoeffCount : longint) : TOGRErr; cdecl; external GDALLib name 'OSRGetTOWGS84';
  function  OSRSetGeogCS(hSRS : TOGRSpatialReferenceH;
            pszGeogName, pszDatumName, pszEllipsoidName : PCHAR;
            dfSemiMajor, dfInvFlattening : double;
            pszPMName : PCHAR; (* = NULL *)
            dfPMOffset : double; (* = 0.0 *)
            pszUnits : PCHAR;(* = NULL *)
            dfConvertToRadians : double(* = 0.0 *)
            ) : TOGRErr; cdecl; external GDALLib name 'OSRSetGeogCS';
  function  OSRGetSemiMajor(hSRS : TOGRSpatialReferenceH;
            pnErr : POGRErr
            ) : double; cdecl; external GDALLib name 'OSRGetSemiMajor';
  function  OSRGetSemiMinor(hSRS : TOGRSpatialReferenceH;
            pnErr : POGRErr
            ) : double; cdecl; external GDALLib name 'OSRGetSemiMinor';
  function  OSRGetInvFlattening(hSRS : TOGRSpatialReferenceH;
            pnErr : POGRErr
            ) : double; cdecl; external GDALLib name 'OSRGetInvFlattening';
  function  OSRSetAuthority(hSRS :  TOGRSpatialReferenceH;
            pszTargetKey : PCHAR;
            pszAuthority : PCHAR;
            nCode : longint
            ) : TOGRErr; cdecl; external GDALLib name 'OSRSetAuthority';
  function  OSRGetAuthorityCode(hSRS : TOGRSpatialReferenceH;
            pszTargetKey : PCHAR
            ) : PCHAR; cdecl; external GDALLib name 'OSRGetAuthorityCode';
  function  OSRGetAuthorityName(hSRS : TOGRSpatialReferenceH;
            pszTargetKey : PCHAR
            ) : PCHAR; cdecl; external GDALLib name 'OSRGetAuthorityName';
  function  OSRSetProjection(hSRS : TOGRSpatialReferenceH;
            pszProjection : PCHAR
            ) : TOGRErr; cdecl; external GDALLib name 'OSRSetProjection';
  function  OSRSetProjParm(hSRS : TOGRSpatialReferenceH;
            pszParmName : PCHAR;
            dfValue : double
            ) : TOGRErr; cdecl; external GDALLib name 'OSRSetProjParm';
  function  OSRGetProjParm(hSRS :  TOGRSpatialReferenceH;
            pszParmName : PCHAR;
            dfDefault : double;
            pnErr : POGRErr
            ) : double; cdecl; external GDALLib name 'OSRGetProjParm';
  function  OSRSetNormProjParm(hSRS : TOGRSpatialReferenceH;
            pszParmName : PCHAR;
            dfValue : double
            ) : TOGRErr; cdecl; external GDALLib name 'OSRSetNormProjParm';
  function  OSRGetNormProjParm(hSRS :  TOGRSpatialReferenceH;
            pszParmName : PCHAR;
            dfDefault : double;
            pnErr : POGRErr
            ) : double; cdecl; external GDALLib name 'OSRGetNormProjParm';
  function  OSRSetUTM(hSRS :  TOGRSpatialReferenceH;
            nZone : longint;
            bNorth : longint
            ) : TOGRErr; cdecl; external GDALLib name 'OSRSetUTM';
  function  OSRGetUTMZone(hSRS :  TOGRSpatialReferenceH;
            pbNorth : PINTEGER
            ) : longint; cdecl; external GDALLib name 'OSRGetUTMZone';
  //** Albers Conic Equal Area */
  function   OSRSetACEA(hSRS :  TOGRSpatialReferenceH;
                        dfStdP1, dfStdP2,
                        dfCenterLat, dfCenterLong,
                        dfFalseEasting, dfFalseNorthing : double
                        ) : TOGRErr; cdecl; external GDALLib name 'OSRSetACEA';

  //** Azimuthal Equidistant */
  function    OSRSetAE(hSRS :  TOGRSpatialReferenceH;
                       dfCenterLat, dfCenterLong,
                       dfFalseEasting, dfFalseNorthing : double
                       ) : TOGRErr; cdecl; external GDALLib name 'OSRSetACEA';

  //** Bonne */
  function   OSRSetBonne(hSRS : TOGRSpatialReferenceH;
                         dfStandardParallel, dfCentralMeridian,
                         dfFalseEasting, dfFalseNorthing : double
                         ) : TOGRErr; cdecl; external GDALLib name 'OSRSetBonne';

  //** Cylindrical Equal Area */
  function   OSRSetCEA(hSRS :  TOGRSpatialReferenceH;
                       dfStdP1, dfCentralMeridian,
                       dfFalseEasting, dfFalseNorthing : double
                       ) : TOGRErr; cdecl; external GDALLib name 'OSRSetCEA';

  //** Cassini-Soldner */
  function   OSRSetCS(hSRS :  TOGRSpatialReferenceH;
                      dfCenterLat, dfCenterLong,
                      dfFalseEasting, dfFalseNorthing : double
                      ) : TOGRErr; cdecl; external GDALLib name 'OSRSetCS';

  //** Equidistant Conic */
  function   OSRSetEC(hSRS :  TOGRSpatialReferenceH;
                      dfStdP1, dfStdP2,
                      dfCenterLat, dCenterLong,
                      dfFalseEasting, dfFalseNorthing : double
                      ) : TOGRErr; cdecl; external GDALLib name 'OSRSetEC';

  ///** Eckert I-VI */
  function   OSRSetEckert(hSRS :  TOGRSpatialReferenceH;
                          nVariation : longint;
                          dfCentralMeridian,
                          dfFalseEasting,
                          dfFalseNorthing : double
                          ) : TOGRErr; cdecl; external GDALLib name 'OSRSetEckert';

  //** Eckert IV */
  function   OSRSetEckertIV(hSRS :  TOGRSpatialReferenceH;
                            dfCentralMeridian,
                            dfFalseEasting,
                            dfFalseNorthing : double
                            ) : TOGRErr; cdecl; external GDALLib name 'OSRSetEckertIV';

  //** Eckert VI */
  function   OSRSetEckertVI(hSRS :  TOGRSpatialReferenceH;
                            dfCentralMeridian,
                            dfFalseEasting,
                            dfFalseNorthing : double
                            ) : TOGRErr; cdecl; external GDALLib name 'OSRSetEckertVI';

  //** Equirectangular */
  function   OSRSetEquirectangular(hSRS : TOGRSpatialReferenceH;
                                   dfCenterLat, dfCenterLong,
                                   dfFalseEasting, dfFalseNorthing : double
                                   ) : TOGRErr; cdecl; external GDALLib name 'OSRSetEquirectangular';
  //** Equirectangular generalized form */
  function   OSRSetEquirectangular2(hSRS : TOGRSpatialReferenceH;
                                    dfCenterLat, dfCenterLong,
                                    dfPseudoStdParallel1,
                                    dfFalseEasting,
                                    dfFalseNorthing : double
                                    ) : TOGRErr; cdecl; external GDALLib name 'OSRSetEquirectangular2';

  //** Gall Stereograpic */
  function   OSRSetGS(hSRS :  TOGRSpatialReferenceH;
                      dfCentralMeridian,
                      dfFalseEasting,
                      dfFalseNorthing : double
                      ) : TOGRErr; cdecl; external GDALLib name 'OSRSetGS';

  //** Goode Homolosine */
  function   OSRSetGH(hSRS :  TOGRSpatialReferenceH;
                      dfCentralMeridian,
                      dfFalseEasting,
                      dfFalseNorthing : double
                      ) : TOGRErr; cdecl; external GDALLib name 'OSRSetGH';

  //** GEOS - Geostationary Satellite View */
  function   OSRSetGEOS(hSRS : TOGRSpatialReferenceH;
                        dfCentralMeridian, dfSatelliteHeight,
                        dfFalseEasting, dfFalseNorthing : double
                        ) : TOGRErr; cdecl; external GDALLib name 'OSRSetGEOS';

  //** Gauss Schreiber Transverse Mercator */
  function   OSRSetGaussSchreiberTMercator(hSRS :  TOGRSpatialReferenceH;
                                           dfCenterLat, dfCenterLong,
                                           dfScale,
                                           dfFalseEasting,
                                           dfFalseNorthing : double
                                           ) : TOGRErr; cdecl; external GDALLib name 'OSRSetGaussSchreiberTMercator';
  //** Gnomonic */
  function   OSRSetGnomonic(hSRS : TOGRSpatialReferenceH;
                            dfCenterLat, dfCenterLong,
                            dfFalseEasting, dfFalseNorthing : double
                            ) : TOGRErr; cdecl; external GDALLib name 'OSRSetGnomonic';
  //** Hotine Oblique Mercator using azimuth angle */
  function   OSRSetHOM(hSRS :  TOGRSpatialReferenceH;
                       dfCenterLat, dfCenterLong,
                       dfAzimuth, dfRectToSkew,
                       dfScale,
                       dfFalseEasting, dfFalseNorthing : double
                       ) : TOGRErr; cdecl; external GDALLib name 'OSRSetHOM';
  //** Hotine Oblique Mercator using two points on centerline */
  function   OSRSetHOM2PNO(hSRS :  TOGRSpatialReferenceH;
                           dfCenterLat,
                           dfLat1, dfLong1,
                           dfLat2, dfLong2,
                           dfScale,
                           dfFalseEasting, dfFalseNorthing : double
                           ) : TOGRErr; cdecl; external GDALLib name 'OSRSetHOM2PNO';
  //** International Map of the World Polyconic */
  function   OSRSetIWMPolyconic(hSRS :  TOGRSpatialReferenceH;
                                dfLat1, dfLat2,
                                dfCenterLong,
                                dfFalseEasting,
                                dfFalseNorthing : double
                                ) : TOGRErr; cdecl; external GDALLib name 'OSRSetIWMPolyconic';

  ///** Krovak Oblique Conic Conformal */
  function   OSRSetKrovak(hSRS :  TOGRSpatialReferenceH;
                          dfCenterLat, dfCenterLong,
                          dfAzimuth, dfPseudoStdParallelLat,
                          dfScale,
                          dfFalseEasting,
                          dfFalseNorthing : double
                          ) : TOGRErr; cdecl; external GDALLib name 'OSRSetKrovak';
  //** Lambert Azimuthal Equal-Area */
  function   OSRSetLAEA(hSRS :  TOGRSpatialReferenceH;
                        dfCenterLat, dfCenterLong,
                        dfFalseEasting, dfFalseNorthing : double
                        ) : TOGRErr; cdecl; external GDALLib name 'OSRSetLAEA';
  //** Lambert Conformal Conic */
  function   OSRSetLCC(hSRS :  TOGRSpatialReferenceH;
                       dfStdP1, dfStdP2,
                       dfCenterLat, dfCenterLong,
                       dfFalseEasting, dfFalseNorthing : double
                       ) : TOGRErr; cdecl; external GDALLib name 'OSRSetLCC';

  //** Lambert Conformal Conic 1SP */
  function   OSRSetLCC1SP(hSRS :  TOGRSpatialReferenceH;
                          dfCenterLat, dfCenterLong,
                          dfScale,
                          dfFalseEasting,
                          dfFalseNorthing : double
                          ) : TOGRErr; cdecl; external GDALLib name 'OSRSetLCC1SP';

  //** Lambert Conformal Conic (Belgium) */
  function   OSRSetLCCB(hSRS :  TOGRSpatialReferenceH;
                        dfStdP1, dfStdP2,
                        dfCenterLat, dfCenterLong,
                        dfFalseEasting, dfFalseNorthing : double
                        ) : TOGRErr; cdecl; external GDALLib name 'OSRSetLCCB';
  //** Miller Cylindrical */
  function   OSRSetMC(hSRS :  TOGRSpatialReferenceH;
                      dfCenterLat, dfCenterLong,
                      dfFalseEasting, dfFalseNorthing : double
                      ) : TOGRErr; cdecl; external GDALLib name 'OSRSetMC';

  //** Mercator */
  function   OSRSetMercator(hSRS :  TOGRSpatialReferenceH;
                            dfCenterLat, dfCenterLong,
                            dfScale,
                            dfFalseEasting, dfFalseNorthing : double
                            ) : TOGRErr; cdecl; external GDALLib name 'OSRSetMercator';

  //** Mollweide */
  function    OSRSetMollweide(hSRS :  TOGRSpatialReferenceH;
                              dfCentralMeridian,
                              dfFalseEasting,
                              dfFalseNorthing : double
                              ) : TOGRErr; cdecl; external GDALLib name 'OSRSetMollweide';

  //** New Zealand Map Grid */
  function   OSRSetNZMG(hSRS :  TOGRSpatialReferenceH;
                        dfCenterLat, dfCenterLong,
                        dfFalseEasting, dfFalseNorthing : double
                        ) : TOGRErr; cdecl; external GDALLib name 'OSRSetNZMG';

  //** Oblique Stereographic */
  function   OSRSetOS(hSRS :  TOGRSpatialReferenceH;
                      dfOriginLat, dfCMeridian,
                      dfScale,
                      dfFalseEasting, dfFalseNorthing : double
                      ) : TOGRErr; cdecl; external GDALLib name 'OSRSetOS';
  //** Orthographic */
  function   OSRSetOrthographic(hSRS :  TOGRSpatialReferenceH;
                                dfCenterLat, dfCenterLong,
                                dfFalseEasting,
                                dfFalseNorthing : double
                                ) : TOGRErr; cdecl; external GDALLib name 'OSRSetOrthographic';
  //** Polyconic */
  function   OSRSetPolyconic(hSRS :  TOGRSpatialReferenceH;
                             dfCenterLat, dfCenterLong,
                             dfFalseEasting, dfFalseNorthing : double
                             ) : TOGRErr; cdecl; external GDALLib name 'OSRSetPolyconic';
  //** Polar Stereographic */
  function   OSRSetPS(hSRS :  TOGRSpatialReferenceH;
                      dfCenterLat, dfCenterLong,
                      dfScale,
                      dfFalseEasting, dfFalseNorthing : double
                      ) : TOGRErr; cdecl; external GDALLib name 'OSRSetPS';
  //** Robinson */
  function   OSRSetRobinson(hSRS :  TOGRSpatialReferenceH;
                            dfCenterLong,
                            dfFalseEasting, dfFalseNorthing : double
                            ) : TOGRErr; cdecl; external GDALLib name 'OSRSetRobinson';
  //** Sinusoidal */
  function   OSRSetSinusoidal(hSRS :  TOGRSpatialReferenceH;
                              dfCenterLong,
                              dfFalseEasting,
                              dfFalseNorthing : double
                              ) : TOGRErr; cdecl; external GDALLib name 'OSRSetSinusoidal';

  //** Stereographic */
  function   OSRSetStereographic(hSRS :  TOGRSpatialReferenceH;
                                 dfCenterLat, dfCenterLong,
                                 dfScale,
                                 dfFalseEasting,
                                 dfFalseNorthing : double
                                 ) : TOGRErr; cdecl; external GDALLib name 'OSRSetStereographic';

  //** Swiss Oblique Cylindrical */
  function   OSRSetSOC(hSRS :  TOGRSpatialReferenceH;
                       dfLatitudeOfOrigin, dfCentralMeridian,
                       dfFalseEasting, dfFalseNorthing : double
                       ) : TOGRErr; cdecl; external GDALLib name 'OSRSetSOC';

  //** Transverse Mercator */
  function   OSRSetTM(hSRS :  TOGRSpatialReferenceH;
                      dfCenterLat, dfCenterLong,
                      dfScale,
                      dfFalseEasting, dfFalseNorthing : double
                      ) : TOGRErr; cdecl; external GDALLib name 'OSRSetTM';

  //** Transverse Mercator variant */
  function   OSRSetTMVariant(hSRS :  TOGRSpatialReferenceH;
                             pszVariantName : PCHAR;
                             dfCenterLat, dfCenterLong,
                             dfScale,
                             dfFalseEasting, dfFalseNorthing : double
                             ) : TOGRErr; cdecl; external GDALLib name 'OSRSetTMVariant';

  //** Tunesia Mining Grid  */
  function   OSRSetTMG(hSRS :  TOGRSpatialReferenceH;
                       dfCenterLat, dfCenterLong,
                       dfFalseEasting, dfFalseNorthing : double
                       ) : TOGRErr; cdecl; external GDALLib name 'OSRSetTMG';

  //** Transverse Mercator (South Oriented) */
  function   OSRSetTMSO(hSRS :  TOGRSpatialReferenceH;
                        dfCenterLat, dfCenterLong,
                        dfScale,
                        dfFalseEasting, dfFalseNorthing : double
                        ) : TOGRErr; cdecl; external GDALLib name 'OSRSetTMSO';

  //** VanDerGrinten */
  function   OSRSetVDG(hSRS :  TOGRSpatialReferenceH;
                       dfCenterLong,
                       dfFalseEasting, dfFalseNorthing : double
                       ) : TOGRErr; cdecl; external GDALLib name 'OSRSetVDG';

  //** Wagner I -- VII */
  function   OSRSetWagner(hSRS :  TOGRSpatialReferenceH; nVariation : longint;
                          dfFalseEasting,
                          dfFalseNorthing : double
                          ) : TOGRErr; cdecl; external GDALLib name 'OSRSetWagner';

procedure OSRCleanup; cdecl; external GDALLib name 'OSRCleanup';


//* -------------------------------------------------------------------- */
//*      Projection transform dictionary query.                          */
//* -------------------------------------------------------------------- */

function OPTGetProjectionMethods : PPCHAR; cdecl; external GDALLib name 'OPTGetProjectionMethods';
function OPTGetParameterList(pszProjectionMethod : PCHAR;
                             ppszUserName : PPCHAR
                             ) : PPCHAR; cdecl; external GDALLib name 'OPTGetParameterList';
function OPTGetParameterInfo(pszProjectionMethod : PCHAR;
                             pszParameterName : PCHAR;
                             ppszUserName : PPCHAR;
                             ppszType : PPCHAR;
                             pdfDefaultValue : PDOUBLE
                             ) : longint; cdecl; external GDALLib name 'OPTGetParameterInfo';

{$IFDEF MSWINDOWS}
  function  OCTNewCoordinateTransformation (hSourceSRS : TOGRCoordinateTransformationH;
            hTargetSRS : TOGRCoordinateTransformationH
            ) : TOGRCoordinateTransformationH ;
            stdcall; external GDALLib name '_OCTNewCoordinateTransformation@8';
  procedure OCTDestroyCoordinateTransformation (const Handle : TOGRCoordinateTransformationH
            ); stdcall; external GDALLib name '_OCTDestroyCoordinateTransformation@4';
  function  OCTTransformEx (const Handle : TOGRCoordinateTransformationH; const PointCount : longint;
            X, Y, Z : PDOUBLE;
            pSuccess : PINTEGER
            ) : longint; stdcall; external GDALLib name '_OCTTransformEx@24';
  function  OCTTransform (const Handle : TOGRCoordinateTransformationH; const PointCount : longint;
            X, Y, Z : PDOUBLE
            ) : longint; stdcall; external GDALLib name '_OCTTransform@20';
{$ENDIF}
{$IFDEF LINUX}
function  OCTNewCoordinateTransformation (hSourceSRS : TOGRCoordinateTransformationH;
          hTargetSRS : TOGRCoordinateTransformationH
          ) : TOGRCoordinateTransformationH ;
          cdecl; external GDALLib name 'OCTNewCoordinateTransformation';
procedure OCTDestroyCoordinateTransformation (const Handle : TOGRCoordinateTransformationH
          ); stdcall; external GDALLib name '_OCTDestroyCoordinateTransformation@4';
function  OCTTransformEx (const Handle : TOGRCoordinateTransformationH; const PointCount : longint;
          X, Y, Z : PDOUBLE;
          pSuccess : PINTEGER
          ) : longint; cdecl; external GDALLib name 'OCTTransformEx';
function  OCTTransform (const Handle : TOGRCoordinateTransformationH; const PointCount : longint;
          X, Y, Z : PDOUBLE
          ) : longint; cdecl; external GDALLib name 'OCTTransform';
{$ENDIF}
implementation

end.

เริ่มโครงการโปรแกรมคำนวณค่าพิกัดด้วย Lazarus

  • เมื่อรัน Lazarus ลองเปิด Project ชื่อ UTM2LatLong.lpi จะเห็นหน้าตาของโปรแกรมที่ผมเขียนไว้ดังรูปด้านล่าง
เปิดโครงการโปรแกรมแปลงพิกัด
  • โปรแกรมนี้มีฟอร์มอยู่ 1 ฟอร์มเท่านั้น ผมวาง combobox ซึ่งเก็บชื่อ datum ซึ่งในโปรแกรมตัวอย่างนี้ผมเขียนไว้แค่ 2 datum คือ WGS84 และ Indian 1975 เท่านั้นโดยแยกย่อยเป็น WGS84 (geographic), Indian 1975 (geographic), WGS84/UTM Zone 47N, WGS84/UTM Zone 48N, Indian 1975/UTM Zone 47N และ Indian 1975/UTM Zone 48N
ฟอร์มของโปรแกรม
  • รู้สึกว่า post ในตอนนี้จะยาวมากไป ติดตามต่อตอนที่ 2

เปรียบเทียบ ASTER GDEM (30m) กับ SRTM DEM (90m)

SRTM DEM

  • SRTM (Shuttle Radar Topography Mission) DEM ของ NASA จัดทำโดย NASA เปิดให้บริการฟรีตั้งแต่ปี 2003 ครอบคลุมพื้นที่ประมาณ 80% ของพื้นที่โลก ขนาด pixel ของ DEM คือ 3 ฟิลิปดา(1 ฟิลิปดา ประมาณ 30 เมตร) หรือขนาดประมาณ 90 เมตร ข้อบกพร่อง ในเบื้องต้นของ SRTM DEM จะมี void หรือรูที่ไม่มีข้อมูลอยู่มากพอสมควร แต่ในรุ่นที่ 4 ปรากฎว่าได้แก้ไขข้อบกพร่องนี้ไปพอสมควร DEM ของ SRTM ผมได้นำไปใช้งาน ช่วยได้มากทีเดียว

ASTER GDEM

  • เป็นโครงการร่วมมือระหว่าง NASA และ METI (The Ministry of Economy, Trade, and Industry ของ ญี่ปุ่น) ตอนนี้ยังเป็น version 1 คาดว่าจะปรับปรุงข้อมูลไปอีกพอสมควร เนื่องจากมีข้อบกพร่องอยู่บ้าง ด้วยความที่ DEM มีขนาด pixel ที่ 1 ฟิลิปดา (ประมาณ 30 เมตร) ทำให้คนคาดหวังว่าจะได้ใช้ DEM ที่มีความละเอียดมากกว่า SRTM DEM แต่หลายๆคน ค่อนข้างผิดหวัง แต่อย่างไรก็ตามของใหม่ก็ต้องดีกว่าของเก่าอย่างแน่นอน เพียงแต่ความคาดหวังมันมากไปหน่อยแค่นั้นเอง
  • กระบวนการจัดทำ DEM เริ่มจากญี่ปุ่นสร้างดาวเทียม Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) และปล่อยสู่วงโคจรโดยกระสวยอวกาศของ NASA ในปี 1999
  • หลักของการทำ DEM คล้ายๆกับงาน Photogrammetry  คือถ่ายภาพเป็นภาพ stereo-pair คือในวงโคจรเดียวกัน จะถ่ายภาพในมุมจุดต่ำสุด(nadir) และถ่ายอีกครั้งเมื่อวงโคจรเลื่อนไปด้านหน้าและมุมถ่ายชี้กลับมาด้านหลัง (backward angle) ทำให้ได้ภาพคู่สามารถนำมาฟอร์มเป็นภาพ 3D และนำมาวัดความสูงได้ด้วยวิธี Automated processing ด้วยโปรแกรม Erdas Imagine software (Othobase Pro)

ASTER Satelite

  • หลักการจัดสร้าง DEM ดูได้จากไดอะแกรมด้านล่าง

ไดอะแกรมการจัดสร้าง DEM

  • ภาพถ่ายดาวเทียม 1 scene กินพื้นที่ประมาณ  60 กม. x 60 กม. การจัดทำ DEM ต้องโคจรประมาณ 4200 แนวเพื่อให้ภาพมี overlap ที่เพียงพอ และแต่ละแนวถ่ายภาพประมาณ 4100 ภาพ สรุปแล้วในการจัดทำ DEM ทั้งโลกนี้ ต้องใช้ภาพประมาณ 4100 x 4200 แต่ที่เป็นภาพคู่นำมาคำนวณความสูงมีประมาณ 1.5 ล้านภาพ ข้อมูลของ GDEM ครอบคลุมตั้งแต่แลตติจูด 83S จนถึง 83N การจัดสร้าง DEM จะแบ่ง GDEM ออกเป็น tile (เหมือนกระเบื้องปูพื้น) 1 tile กินพื้นที่ 1 องศา X 1 องศา (ประมาณ 108 กม. X 108 กม.) แล้วเปิดบริการให้ดาวน์โหลดไปใช้งานโดยไม่ต้องเสียค่าใช้จ่าย เปิดให้ดาวน์โหลดไปใช้งานได้ประมาณ ครึ่งปีแล้วครับ (ตั้งแต่กลางปี 2009)

ตารางเปรียบเทียบระหว่าง ASTER DEM และ SRTM DEM

ASTER GDEM SRTM3*
Data source ASTER Space shuttle radar
Generation and distribution METI/NASA NASA/USGS
Release year 2009 ~ 2003 ~
Data acquisition period 2000 ~ ongoing 10  days (in 2000)
Posting interval 30m 90m
DEM accuracy (stdev.) 7~14m 10m
DEM coverage 83 degrees north ~ 83 degrees south 60 degrees north ~ 56 degrees south
Area of missing data Areas with no ASTER data due to constant cloud cover (supplied by other DEM) Topographically steep area (due to radar characteristics)
  • จากตารางด้านบนนำมาจาก website ของ ASTER GDEM ดูเหมือนว่า ASTER GDEM จะดีกว่า SRTM DEM ทุกๆด้าน

การ Download GDEM

ลงทะเบียนเพื่อดาวน์โหลด GDEM

วิธีการดาวน์โหลด GDEM

  • การเลือกดาวน์โหลดมีวิธีอยู่ 4 วิธีคือ เลือกโดยตรง(ด้วยการคลิกเลือก tile), เลือกจากการวาด polygon, เลือกจากกำหนด shape file และเลือกจากการกำหนดค่าพิกัด สำหรับผมเองเลือกจากการคลิกที่ tile สะดวกที่สุด

เลือกวิธีการดาวน์โหลด GDEM

  • เมื่อขยายรูปได้ตามที่ต้องการแล้ว คลิกที่ปุ่มกริด (Grid #) เพื่อแสดงเส้นกริดจะได้แบ่งรูป tile ได้ชัดเจน ดูรูปด้านล่าง

คลิกที่ tile เพื่อเลือก GDEM ที่จะดาวน์โหลด

  • การจะเลือกพื้นที่ที่ต้องการ แรกให้ขยายรูปก่อนด้วยการเลื่อนสกรอลล์ของเมาส์ (ปุ่มกลาง) เพื่อขยายรูป แล้วคลิกที่ tile ได้ขอบเขตตามที่ต้องการก็คลิกที่ปุ่ม Next ไปต่อได้เลย หน้าต่อไปจะเป็น list ของชื่อไฟล์ของ DEM ถ้าไม่มีอะไรก็คลิก Next ไปต่อได้เลย

แสดง list ของไฟล์ GDEM ที่ต้องเลือกเพื่อจะดาวน์โหลด

  • หน้าต่อไปผมเข้าใจว่าทาง website คงต้องการเก็บสถิติ ว่าผู้ใช้ดาวน์โหลดไปใช้ในกิจการด้านไหน พยายามเลือกให้ใกล้เคียงครับ

เลือกตอบว่านำ GDEM ไปใช้ในด้านใด

  • สุดท้ายมาถึงหน้าที่ต้องการจริงๆคือหน้าดาวน์โหลด จะเห็นชื่อไฟล์ของ GDEM ลิสต์เรียงกันอยู่ด้านซ้าย (เลข 1) และปุ่มดาวน์โหลดแต่ละไฟล์ด้านขวา (เลข 2) แต่ถ้าต้องการดาวน์โหลดทั้งหมดให้คลิกที่เลข 3 ถ้าค้นหาใน google จะพบว่าการดาวน์โหลด GDEM จะช้า แต่ที่ผมลองมาก็ OK เร็วพอใช้ได้ (ช่วงที่ผมใช้  server อาจจะ่ว่างพอดี)

ดาวน์โหลดไฟล์ GDEM

รูปแบบไฟล์ของ GDEM

  • ไฟล์ของ GDEM เมื่อดาวน์โหลดจะได้ไฟล์ zip เช่นที่ latitude 13N และ longitude 97E จะได้ไฟล์ชื่อ ASTGTM_N13E097.zip ชื่อไฟล์ค่อนข้างง่ายบอกพิกัดของ DEM มาให้เสร็จสรรพ เมื่อ unzip จะได้ 2 ไฟล์ที่เป็น GeoTiff (signed 16 bit ระบบค่าพิกัดเป็น geographic บน WGS84) คือ ASTGTM_N13E097_dem.tif และ ASTGTM_N13E097_num.tif ไฟล์หลังเรียกว่า QA Plane file ซึ่งผมยังไม่เข้าใจนักว่าใช้ทำอะไร

เปรียบเทียบ ASTER GDEM แและ SRTM DEM

  • ผมจะเปรียบเทียบพื้นที่ตัวอย่างคือดาวน์โหลดไฟล์ ASTGTM_N13E098_dem.tif ของ ASTER GDEM ในพื้นที่เดียวกันผมดาวน์โหลดไฟล์ SRTM_56_10.tif ของ SRTM DEM (การดาวน์โหลด SRTM DEM เคยกล่าวไปตอนก่อนแล้ว) ผมเปิดไฟล์ของ DEM ด้วย MicroDEM
  • จากรูปด้านบนเป็นการแสดงผลของ SRTM DEM จากโปรแกรม MicroDEM บริเวณปากแม่น้ำทวาย พม่า จะเห็นว่า SRTM จัดเป็นข้อมูล void คือเป็นแม่น้ำและทะเล ซึ่งก็ถือว่าถูกต้อง (แต่การแสดงข้อมูลเป็น void เช่นบนภูเขาจะไม่ถูกต้อง)

  • รูปด้านบนเป็น ASTER GDEM ที่บริเวณเดียวกัน ดูความสูงของภูเขาตรงแหลมทวายก็ดูคมชัดเพราะละเอียดกว่าและก็ OK ดี แต่ที่ผิดร้ายแรงคือตรงแม่น้ำทวายที่ผมวงให้ ร่องน้ำหรือแม่น้ำที่สามารถใช้เดินเรือประมงเข้าไปในเมืองทวายได้ จุดนี้กลับมีความสูงหรือประมาณว่าเป็นเนินเขามางอกอยู่กลางแม่น้ำวัดได้สูงประมาณ 80 เมตรซึ่งไม่ถูกต้องอย่างแน่นอน แต่อย่างที่ผมกล่าวไปแล้ว DEM ชุดนี้ยังเพิ่งให้บริการเท่านั้นยังมีเวลาแก้ไขและปรับปรุงให้ถูกต้องได้อีกมาก

การเปรียบเทียบ ASTER GDEM และ STRM DEM ในทัศนะอื่นๆ

  • ผมไม่ใช่นักวิชาการจึงเปรียบเทียบด้วยการดูด้วยสายตาเปล่า ลองไปดูทัศนะคนอื่นๆที่น่าสนใจได้แก่ http://www.gisdevelopment.net/technology/tm/tm001a.htm และอีกท่านหนึ่ง Arrowsmith blog

การเขียนตัวเลขค่าระดับเข้า Autocad ด้วยโปรแกรม Spot Fire (ฟรี)

  • ช่วงเดือนธันวาคม 2009 ส่งท้ายปีต้องติดภารกิจของบริษัทฯ เข้าไปพม่าที่เมืองมะริดอยู่ 2 ครั้ง ทำให้ไม่มีเวลาเขียน blog เลย ปีนี้ 2010 คิดว่าเกือบครึ่งปีที่จะต้องไปอยู่ที่มะริดและทวาย ของประเทศพม่า ซึ่งเป็นเมืองปิดมานาน ไม่รับนักท่องเที่ยว คนที่มีโอกาสเข้าไปโดยเฉพาะคนไทย มีน้อยมาก สองเมืองนี้ตอนนี้น่าอยู่มาก เงียบ สงบ คล้ายๆกับบ้านเกิดของผมเมื่อสัก 40 ปีที่แล้ว มองไปทางไหนยังเห็นไร่นาอยู่มาก แต่สิ่งที่เปลี่ยนแปลงไปก็คือความเจริญด้านวัตถุ รถยนต์ก็ยังเห็นบ้างแต่ไม่มาก รถมอเตอร์ไซค์ ก็หนาตา แต่ที่ประทับใจคือคนที่นั่นยังศรัทธาต่อศาสนาพุทธอย่างล้ำลึก
  • แต่ที่เงียบกว่าคือเมืองทวายความเจริญยังน้อยกว่ามะริดมาก ภาพของขี่จักรยานยังเห็นอยู่หลากหลาย ผมมีโอกาสได้ไปนมัสการรอยพระพุทธบาททั้งสองรอยคือรอยเท้าซ้ายและขวา เป็นสิ่งที่ประทับใจมากๆ

โปรแกรม Spot Fire

  • มาเข้าเรื่องโปรแกรม Spot Fire เป็นโปรแกรมเล็กๆ ที่เขียนด้วย Delphi เขียนไว้หลายปีแล้ว ถ้าเรามีจุดค่าพิกัดและระดับ (X,Y,Z) เก็บอยู่ในรูป Text file และบางทีอาจจะคั่นด้วยเครื่องหมายคอมมาหรือคั่นด้วยช่องว่างก็ตาม สามารถใช้โปรแกรม Land Desktop หรือ Terramodel เขียนขึ้นไม่ยาก แต่จะให้สวยและแยกสีตามค่าระดับนั้นยาก
  • ที่เอาไปใช้บ่อยการเขียนตัวหนังสือแบบนี้จะเป็นงาน Hydrographic Survey มากกว่า ผมจำได้ว่าโปรแกรม HydroNav ของ Trimble รุ่นนั้นยังเป็น DOS อยู่ เขียนตัวหนังสือ (Text) ได้สวยงามมาก แต่ตอนหลังเป็นยุคของวินโดส์  Trimble เลยทิ้ง HydroNav หันไปซื้อ Terramodel สำหรับ windows มาจาก Spectra ถ้าจำไม่ผิด เพื่อมาแทน HydroNav

ตัวอย่างแผนที่ภูมิประเทศ(Topographic Map) ที่เขียนตัวเลขแสดงความสูงเป็นลัำกษณะกริด

  • จากรูปด้านบนจะเห็นค่าระดับเป็นลักษณะกริด จะมีค่าระดับเป็นหน่วยเมตรด้านหน้าและตัวทศนิยมเป็นตัวห้อยขนาดเล็กกว่าอยู่ด้านหลัง และส่วนใหญ่จะเอียงหมุนไปด้านซ้ายประมาณ 45 องศา เพื่อเลี่ยงตัวเลขทับกัน

แสดงตัวเลขค่าระดับบนแผนที่ Hydrographic Map

  • แผนที่แสดงใน Hydrographic Map ทิศทางค่าระดับ(แกน Z) จะตรงข้ามกับงานบนบก คือชี้ลงไปในน้ำเป็นบวก ขึ้นฝั่งเป็นลบ

ดาวน์โหลดโปรแกรมและติดตั้ง

  • สนใจดาวน์โหลดโปรแกรมได้ที่นี่ http://www.4shared.com/file/190515733/31c9f3b5/SpotFireSetup.html เมื่อดาวน์โหลดแล้วำทำการติดตั้งโปรแกรมจะถูกติดตั้งอยู่ที่โฟลเดอร์ c:\program files\survey suite\spot fire และจะมีโฟลเดอร์ย่อยชื่อ examples เก็บไฟล์ตัวอย่างไว้ ซึ่งผู้ใช้สามารถลองดูได้ เมื่อเปิดโปรแกรมมาจะเห็นดังรูปด้านล่าง

Spot Fire หลังจากเปิดโปรแกรม

  • หน้าตาของโปรแกรมเรียบง่ายมี 2 แท็ปคือ File และ Options เมื่อเปิด Text file ที่เก็บจุดที่ต้องการจะแสดงที่แท็ป File ส่วนแท็ป Options จะตั้งค่าตัวเลือกเช่นขนาดตัวหนังสือ มาตราส่วน แยกสีตามแบนด์

การใช้โปรแกรม

  • คลิกที่ไอคอนเปิดไฟล์ browse ไปที่ c:\program files\survey suite\spot fire\examples เปิดไฟล์ชื่อ STOCKPILE-PNEZD.csv อย่าลืมคลิก file type เป็น P,N,E,Z,D ตามรูปแบบของไฟล์ข้อมูล ดังรูปด้านล่าง

เปิด Text file เก็บค่าพิกัด

  • โปรแกรมจะอ่านข้อมูลจาก Text file แล้วเติมลงในตารางกริดดังรูปด้านล่าง

แสดงข้อมูลในตารางกริด

ตั้งค่า Options

  • คลิกที่แท็ป Options เพื่อตั้งค่า รูปด้านล่างเป็นค่าปริยาย

การตั้งค่า Options

  • จากรูปด้านบน ตรงชี้หมายเลข 1 เป็นการตั้งค่าขนาดตัวอักษร ตัวอย่างความสูงตัวอักษร 1.8 มิลลิเมตร (ตัวห้อยหรือทศนิยม ขนาด 1.2 มิลลิเมตร) ที่มาตราส่วน 1:1000 แต่ถ้าเปลี่ยนเป็นมาตราส่วน 1:500 ขนาดตัวอักษรที่ปรากฎใน Model ของ Autocad จะวัดได้ 1.8/2 = 0.9 ถ้าเป็นมาตราส่วน 1:2000 จะวัดขนาดตัวอักษรได้ 2 x 1.8 = 3.6 ที่สำคัญคือเมื่อเรา วัดที่ layout จะวัดได้เท่ากับ 1.8 มม. เสมอ ไม่ว่าจะใช้อัตราส่วนเท่าไหร่
  • ส่วนอื่นๆเช่นมาตราส่วน Horizontal scale, Depth Annotation layer ตั้งชื่อ layer ที่จะเก็บตัวหนังสือ, Points layer ต่อไปคือรูปแบบตัวหนังสือ Text Font ซึ่งใช้ font ของ Autocad คือ isoct.shx หรืออาจจะเป็น font อื่นๆก็นได้เช่น Romans.shx และที่กล่าวไปแล้วคือ Text height และที่สำคัญอีกตัวคือ Rotation Angle มุมที่จะหมุนกระดกตัวหนังสือจัดไว้ประมาณ 45 องศากำลังสวย ส่วน Minus Style คือลักษณะของเครื่องหมายลบ ถ้าเป็นงาน Hydrographic Survey เครื่องหมายลบจะขีดเส้นตรงสั้นๆใต้ตัวเลขความสูง และตัวเลือกตัวสุดท้ายคือ Select point every เลือกจุด 1 จุดทุกๆ 2 จุด (คำพูดง่ายๆคือเลือกจุดเว้นจุด ถ้าเลือกมา 1 จุด เว้นไปสองจุด ก็ป้อนเลข 3)

การเขียนตัวเลขค่าระดับแบบแยกสีตามความลึก

  • การเขียนแยกสีอักษรตามความลึกของจุด (Color Banding) สามารถทำได้โดยการคลิกที่ Automatic color banding interval แล้วตั้งค่า Band เช่นทุกๆ 2 เมตร ดูรูปด้านล่างประกอบ ครั้งแรกเลือกมาตราส่วนเป็น 1:500 เมื่อป้อนค่า interval เสร็จ ให้คลิกที่ปุ่ม Apply

ตั้งค่า options

การ Pump ตัวหนังสือแสดงค่าระดับเข้า Autocad

  • เมื่อทุกสิ่งทุกอย่างพร้อม ก็พร้อมจะปั๊มตัวหนังสือเข้า Autocad ได้คลิกที่ Icon แล้ว confirm ดังรูปด้านล่าง

  • โปรแกรม Spot Fire พยายามเปิด Autocad ให้ แล้ว pump ข้อมูลเข้าไปดังรูป

แสดงภาพเมื่อ pump ข้อมูลเข้า Autocad

  • ที่ Autocad  เรา Zoom เข้าไปดูตัวอักษร

ภาพขยาย Autocad แสดง Text ที่ได้จากการ pump จากโปรแกรม SpotFire

  • ที่ Autocad แล้ว save เป็น Block แล้วนำไป insert เข้ากับแผนที่่ Hydrographic Map จะได้แผนที่ที่แสดงความลึกตามที่ต้องการ
  • โปรแกรม Spot Fire เป็น tools ตัวเล็กๆ ที่ทีมงานผมยังใช้งานอยู่ประจำ ถึงแม้ไม่มีอะไรหวือหวา ก็ขอฝากไว้ด้วยนะครับ

การเขียนโปรแกรมทดสอบการใช้ไลบรารี GDAL/OGR

  • ตอนก่อนผมแนะนำไลบรารี GDAL/OGR ไปพอสมควร ตอนนี้มาเริ่มลองโปรแกรมมิ่งดูกัน โปรแกรมทดสอบผมดัดแปลงจากโค๊ดภาษาซี เป็น Lazarus ดูรายละเอียดโค๊ดภาษาซีได้ที่นี่ http://www.gdal.org/gdal_tutorial.html ส่วนไลบรารีส่วนมากแปลงจาก VB6

Download sourcecode

  • สนใจก็ดาวน์โหลดได้ที่นี่ GDALTest1.zip
  • ไลบราีรีที่ผมเขียน wrapper มามีทั้งหมด 11 ไฟล์ มีไฟล์ gdal.pas, gdalcore.pas, ogr.pas, ogrcore.pas, gdaldriver.pas, gdaldrivermanager.pas, gdalrasterband.pas, gdaldataset.pas, ogrspatialreference.pas, ogrcoordinatetramsformation.pas, gdalcolortable.pas ทั้งหมดนี้รวมอยูในไฟล์ zip ที่ดาวน์โหลดได้จากที่ผมกล่าวไปแล้ว

เปิดโครงการ GdalTest1

  • เมื่อ unzip มาแล้วให้ copy ไปลงโฟลเดอร์ เปิด Lazarus แล้วใช้เมนู Project > Open Project… แล้ว browse ไปที่ไฟล์ gdaltest1.lpi เปิดโครงการมาจะเห็นดังรูปด้านล่าง

เปิดโครงการ gdaltest1

ดู Mainform

  • ที่ Mainform ชื่อไฟล์คือ main.pas และชื่อของฟอร์มคือ frmmain ผมวาง object ลงบนฟอร์มมี MainMenu1, OpenDialog1, SaveDialog1 และ Memo1 ที่ MainMenu1 ผมสร้างเมนูมา 2 เมนูคือ File และ Tools ดูรูปด้านล่าง

เมนู file

เมนู Tools

จุดมุ่งหมายของโปรแกรม

  • ใน sourcecode ที่ให้ดาวน์โหลดจะมีไฟล์เพื่อทดสอบ 2 ไฟล์คือ Test_dem.tif และ Test_image.tif ไฟล์แรกเป็น dem ไฟล์ที่สองเป็นไฟล์ที่ได้จาก Landsat7 เราจะรันโปรแกรมแล้วใช้เมนู Files > Open เปิดไฟล์มาทำการทดสอบ แล้วใช้เมนู Tools ทำการทดสอบเช่น Read Test, Coordinate System Information, List Drivers, CreateCopy Test และ Create Test

ไฟล์รูปที่ใช้ในการทดสอบ

  • เมื่อเปิดดูด้วย ER Viewer จะเห็นไฟล์รูปทั้งสองดังรูปด้านล่าง

ไฟล์รูปแสดงด้วย ER Viewer

Sourcecode ของ main.pas

  • ต่อไปมาดู sourcecode ของฟอร์ม main.pas กัน

[sourcecode language=”delphi”]
unit main;

{$mode objfpc}{$H+}

interface

uses
Classes, SysUtils, FileUtil, LResources, Forms, Controls, Graphics, Dialogs,
StdCtrls, Menus, gdal, gdalcore, gdaldrivermanager, gdaldriver,
gdaldataset, ogr, ogrcore, ogrspatialreference, ogrcoordinatetransformation,
gdalrasterband, gdalcolortable;

type
TCharSet = set of char;

{ TfrmMain }

TfrmMain = class(TForm)
MainMenu1: TMainMenu;
Memo1: TMemo;
mnuCreate: TMenuItem;
mnuCCTest: TMenuItem;
mnuListDrivers: TMenuItem;
mnuCSInfo: TMenuItem;
mnuToolsTest1: TMenuItem;
mnuFileExit: TMenuItem;
mnuFileOpen: TMenuItem;
mnuTools: TMenuItem;
mnuFile: TMenuItem;
OpenDialog1: TOpenDialog;
SaveDialog1: TSaveDialog;
procedure FormCreate(Sender: TObject);
procedure mnuCCTestClick(Sender: TObject);
procedure mnuCreateClick(Sender: TObject);
procedure mnuCSInfoClick(Sender: TObject);
procedure mnuFileExitClick(Sender: TObject);
procedure mnuFileOpenClick(Sender: TObject);
procedure mnuListDriversClick(Sender: TObject);

procedure mnuToolsTest1Click(Sender: TObject);
private
{ private declarations }
function ReportCorner(CornerName : string; pixel, line : double;
gt : TDynArrayDouble; ct : TOGRCoordinateTransformation) : string;
procedure Tokenize(const tokens: TStrings; const sztarget: string; szdelims: TCharSet);

public
{ public declarations }
end;

var
frmMain: TfrmMain;

implementation

procedure TfrmMain.Tokenize(const 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 &lt;= 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.FormCreate(Sender: TObject);
begin
gdalcore.AllRegister;
gdalcore.SetConfigOption(‘GDAL_DATA’, ‘E:\SourceCodes\Lazarus\GDALTest1’);
end;

procedure TfrmMain.mnuCSInfoClick(Sender: TObject);
var
ds : TGDALDataset;
hds : TGDALDatasetH;
Geotransform : array[0..5] of double;
srs : TOGRSpatialReference;
latlong_srs : TOGRSpatialReference;
ct : TOGRCoordinateTransformation;
WKT : PCHAR;
szTopleft, szTopRight, szBottomLeft, szBottomRight, szCenter : string;
szpp : PPCHAR;
i : integer;
tokens : TStringList;
t1, t2 : PCHAR;
obj : TOGRCoordinateTransformationH;

begin
Memo1.Clear;
GDALCore.GDALAllRegister;
try
ds := TGDALDataset.Create;
hds := ds.Open(PCHAR(OpenDialog1.Filename), GA_ReadOnly);
ds.CInit(hds, true);
if (ds.IsValid) then
begin
try
srs := TOGRSpatialReference.Create;
latlong_srs := TOGRSpatialReference.Create;
ct := TOGRCoordinateTransformation.Create;
ds.GetGeoTransform(Geotransform);

Memo1.Lines.Add (‘Size: ‘ + inttostr(ds.XSize) + ‘P x ‘ + inttostr(ds.YSize) + ‘L’);

// report projection in pretty format.
WKT := ds.GetProjectionRef;
if Length(WKT) &gt; 0 then
begin
try
tokens := TStringList.Create;
Memo1.Lines.Add(‘Projection: ‘);
srs.SetFromUserInput(WKT);
szpp := srs.ExportToPrettyWkt(0);
t1 := szpp[0];
t2 := szpp[1];
tokenize(tokens, t1, [#10]);

for i := 0 to tokens.count -1 do
Memo1.Lines.Add(tokens.Strings[i]);

finally
tokens.Free;
end;

if (srs.GetAttrValue(‘PROJECTION’, 0) &lt;&gt; NIL) then
begin
latlong_srs := srs.CloneGeogCS;
obj := ogrcore.OCTNewCoordinateTransformation(srs.GetObjPtr,
latlong_srs.GetObjPtr);
if (obj &lt;&gt; NIL) then
ct.CInit(obj, 1);

end;
end;

Memo1.Lines.Add(‘Origin: ‘ + floattostr(Geotransform[0]) + ‘,’ + floattostr(Geotransform[3]));
szTopLeft := #13#10 + ReportCorner(‘Top Left ‘, 0, 0, Geotransform, ct) + #13#10;
szTopRight := ReportCorner(‘Top Right ‘, ds.XSize, 0, Geotransform, ct) + #13#10;
szBottomLeft := ReportCorner(‘Bottom Left ‘, 0, ds.YSize, Geotransform, ct) + #13#10;
szBottomRight := ReportCorner(‘Bottom Right ‘, ds.XSize, ds.YSize, Geotransform, ct) + #13#10;
szCenter := ReportCorner(‘Center ‘, ds.XSize / 2, ds.YSize / 2, Geotransform, ct) + #13#10;
Memo1.Lines.Add(‘Pixel Size: ‘ + floattostr(Geotransform[1]) + ‘x’ + floattostr(-1 * Geotransform[5])
+ szTopLeft + szTopRight + szBottomLeft + szBottomRight + szCenter);
finally
srs.Free;
latlong_srs.Free;
ct.Free;
end;
end
else
Memo1.Lines.Add(gdalcore.GetLastErrorMsg);
finally
ds.Free;
end;
end;

procedure TfrmMain.mnuFileExitClick(Sender: TObject);
begin
Close;
end;

procedure TfrmMain.mnuFileOpenClick(Sender: TObject);
begin
OpenDialog1.Execute;
Memo1.Clear;
Memo1.Lines.Add(‘Filename ‘ + OpenDialog1.Filename + ‘ selected.’);
end;

//Tested OK.
procedure TfrmMain.mnuListDriversClick(Sender: TObject);
var
drvm : TGDALDriverManager;
drvIndex, drvCount : integer;
szMsg : string;
begin
Memo1.Clear;
Memo1.Lines.Add(‘GDAL_DATA = ‘ + gdalcore.GetConfigOption(‘GDAL_DATA’, NIL));
try
drvm := TGDALDriverManager.Create;

if (drvm.GetDriverCount &lt; 1) then
gdalcore.AllRegister;

drvCount := drvm.GetDriverCount;
Memo1.Lines.Add(‘Count = ‘ + inttostr(drvCount));
//Wiered for loop &quot;i&quot; value run over drvCount-1
//for i := 0 to drvCount – 1 do
drvIndex := 0;
while (drvIndex &lt;= drvCount – 1) do
begin
drvm.Handle := drvm.GetDriver(drvIndex);
if (drvm.GetMetadataItem(DCAP_CREATE, NIL) = ‘YES’)
or (drvm.GetMetadataItem(DCAP_CREATECOPY, NIL) = ‘YES’) then
szMsg := ‘ (Read/Write)’
else
szMsg := ‘ (ReadOnly)’;
inc(drvIndex);
Memo1.Lines.Add(drvm.GetDescription + ‘: ‘ + drvm.GetMetadataItem(gdal.DMD_LONGNAME, NIL) + szMsg);
end;

finally
drvm.Free;
end;
end;

//Tested OK.
procedure TfrmMain.mnuToolsTest1Click(Sender: TObject);
var
ds : TGDALDataset;
hDS : TGDALDatasetH;
szFilename : string;
md : PPCHAR;
i, j, err, iColor : integer;
n1, n2 : double;
bandH : TGDALRasterBandH;
band : TGDALRasterBand;
//RawData : TDynArrayInt16;
ct : TGDALColorTable;
hct : TGDALColorTableH;
CEntry : array[0..3]of integer;
szBlckSX, szBlckSY, szOffset, szScale ,szMin, szMax, szDataType : string;

begin
Memo1.Clear;
szFilename := OpenDialog1.Filename;
GDALCore.GDALAllRegister;
try
ds := TGDALDataset.Create;
hDS := ds.Open(PCHAR(szFilename), GA_ReadOnly);
ds.CInit(hDS, true);
if (ds.IsValid) then
begin
Memo1.Lines.Add(‘Open succeeded.’);
Memo1.Lines.Add(‘Size: ‘ + inttostr(ds.XSize) + ‘x’ + inttostr(ds.YSize)
+ ‘x’ + inttostr(ds.BandCount));

md := ds.GetMetadata(NIL);
if (md &lt;&gt; NIL) then
begin
Memo1.Lines.Add(‘Metadata:’);
i := 0;
while md[i]&lt;&gt;NIL do
begin
Memo1.Lines.Add(‘ ‘ + md[i]);
Inc(i);
end;
end;

for i := 1 To ds.BandCount do
begin
bandH := ds.GetRasterBand(i);
try
band := TGDALRasterBand.create;
band.CInit(bandH);
szDataType := band.GetDataTypeName(band.DataType);
szBlckSX := inttostr(band.BlockXSize);
szBlckSY := inttostr(band.BlockYSize);
szOffset := floattostr(band.GetOffset);
szScale := floattostr(band.GetScale);
szMin := floattostr(band.GetMinimum);
szMax := floattostr(band.GetMaximum);

Memo1.Lines.Add(‘Band ‘ + inttostr(i)
+ ‘ BlockSize ‘ + szBlckSX + ‘x’ + szBlckSY+#13#10
+ ‘ DataType=’ + szDataType + #13#10
+ ‘ Offset=’ + szOffset + #13#10
+ ‘ Scale=’ + szScale+ #13#10
+ ‘ Min=’ + szMin + #13#10
+ ‘ Max=’ + szMax);

//setlength(RawData, ds.XSize);
//j := ds.XSize * sizeof(band.DataType);
//err := band.RasterIO(gdal.GF_Read, 0, 0, ds.XSize, 1, @RawData[0], band.DataType);
//n1 := RawData[0];
//n2 := RawData[10];
//Memo1.Lines.Add (‘ Random Data: ‘ + floattostr(n1) + ‘ ‘ + floattostr(n2));
try
ct := TGDALColorTable.Create;
ct := band.GetColorTable;
if (ct.IsValid) then
begin
Memo1.Lines.Add(‘ Has Color Table, ‘ + inttostr(ct.EntryCount) + ‘ entries’);
for iColor := 0 to ct.EntryCount – 1 do
begin
ct.GetColorEntryAsRGB(iColor, CEntry);
Memo1.Lines.Add(‘ ‘ + inttostr(iColor) + ‘: ‘
+ inttostr(CEntry[0]) + ‘,’ + inttostr(CEntry[1]) + ‘,’
+ inttostr(CEntry[2]) + ‘,’ + inttostr(CEntry[3]));
end;
end;
finally
ct.Free;
end;
finally
band.Free;
end;
end;
end
else
Memo1.Lines.Add(gdalcore.GetLastErrorMsg);

finally
ds.Free;
end;
end;

//Tested OK.
procedure TfrmMain.mnuCCTestClick(Sender: TObject);
var
SrcFilename, DstFilename : PCHAR;
Drv : TGDALDriver;
SrcDS : TGDALDataset;
DstDS : TGDALDataset;
hSrc : TGDALDatasetH;
hDRV : TGDALDriverH;
DrvM : TGDALDriverManager;

begin
Memo1.Clear;
SrcFilename := PCHAR(OpenDialog1.Filename);

GDALCore.GDALAllRegister;
try
SrcDS := TGDALDataset.Create;
DstDS := TGDALDataset.Create;
DrvM := TGDALDriverManager.Create;
Drv := TGDALDriver.Create;

hSrc := SrcDS.Open(SrcFilename, GA_ReadOnly);
SrcDS.CInit(hSrc, true);
//Erdas Imagine (*.img)
hDrv := DrvM.GetDriverByName(‘HFA’);

Drv.CInit(hDrv);

SaveDialog1.Title := ‘Save to…’;
SaveDialog1.Execute;

DstFilename := PCHAR(SaveDialog1.FileName);
DstDS := Drv.CreateCopy(DstFilename, SrcDS, True, NIL);
if (DstDS.IsValid = true) then
begin
ShowMessage(‘CreateCopy Succeeded, output is ‘ + DstFilename);
Memo1.Lines.Add (‘CreateCopy Succeeded, output is ‘ + DstFilename);
end
else
ShowMessage(‘Create Copy Failed: ‘ + gdalcore.GetLastErrorMsg);
finally
DstDS.Free;
DrvM.Free;
Drv.Free;
SrcDS.Free;
end;
end;

procedure TfrmMain.mnuCreateClick(Sender: TObject);
var
SrcFilename, DstFilename : PCHAR;
Drv : TGDALDriver;
dsSrc, dsDst : TGDALDataset;
hBand, hSrcBand, hDstBand : TGDALRasterBandH;
hDrv : TGDALDriverH;
hSrc : TGDALDatasetH;
DrvM : TGDALDriverManager;
err : TCPLErr;
gt : array[0..5] of double;
Scanline : array of double;
SrcBand, DstBand : TGDALRasterBand;
iBand, iLine : longint;
NoDataValue : double;
pSuccess : PINTEGER;
ct : TGDALColorTable;
ct2 : TGDALColorTable;
iColor : integer;
Tuple : array[0..3] of integer;
eDatatype : TGDALDatatype;
begin
Memo1.Clear;
try
dsSrc := TGDALDataset.Create;
dsDst := TGDALDataset.Create;
DrvM := TGDALDriverManager.Create;
SrcBand := TGDALRasterband.Create;
Drv := TGDALDriver.Create;
DstBand := TGDALRasterband.Create;

SrcFilename := PCHAR(OpenDialog1.Filename);

SaveDialog1.Title := ‘Save as to’;
SaveDialog1.DefaultExt := ‘tif’;
SaveDialog1.Execute;
DstFilename := PCHAR(SaveDialog1.FileName);

GDALCore.GDALAllRegister;
hSrc:= dsSrc.Open(SrcFilename, GA_ReadOnly);
dsSrc.CInit(hSrc, true);
//Use to find RasterBand Datatype
hBand := dsSrc.GetRasterBand(1);
SrcBand.CInit(hBand);
eDataType := SrcBand.GetRasterDataType(hBand);

if (hSRC = NIL) then
begin
Memo1.Lines.Add(gdalcore.GetLastErrorMsg);
exit;
end;
//ดูชื่อ Driver ได้จาก List Driver
hDrv := DrvM.GetDriverByName(‘GTiff’);
Drv.CInit(hDrv);
if (hDrv &lt;&gt; NIL) then
begin
dsDst := Drv.GDALCreate(DstFilename, dsSrc.XSize, dsSrc.YSize,
dsSrc.BandCount, eDatatype, NIL);

if (dsDst.Isvalid) then
Memo1.Lines.Add(‘Create Succeeded, file is ‘ + DstFilename)
else
begin
Memo1.Lines.Add(‘Create Failed: ‘ + gdalcore.GetLastErrorMsg);
exit;
end;
end;
//Copy geotransform

err := dsSrc.GetGeoTransform(gt);
if (err = CE_None) then
dsDst.SetGeoTransform(gt);

//Copy projection
dsDst.SetProjection(dsSrc.GetProjectionRef);

//Copy metadata.
dsDst.SetMetadata(dsSrc.GetMetadata(”), ”);

//Copy band info
setlength(Scanline, dsSrc.XSize);
for iBand := 1 to dsSrc.BandCount do
begin
hSrcBand := dsSrc.GetRasterBand(iBand);
hDstBand := dsDst.GetRasterBand(iBand);
SrcBand.CInit(hSrcBand);
DstBand.CInit(hDstBand);

DstBand.SetMetadata(SrcBand.GetMetadata(”), ”);
DstBand.SetOffset(SrcBand.GetOffset);
DstBand.SetScale(SrcBand.GetScale);

NoDataValue := SrcBand.GetNoDataValue(@pSuccess);
if (pSuccess &lt;&gt; NIL) then
DstBand.SetNoDataValue(NoDataValue);

//Copy Paletted if one present.
ct := SrcBand.GetColorTable;
if (ct.IsValid) then
// We manually copy the color table. This isn’t really
// necessary, but gives us a chance to try out all color
// table methods.
begin
ct2 := TGDALColorTable.Create;
for iColor := 0 to ct.EntryCount do
begin
ct.GetColorEntryAsRGB(iColor, Tuple);
ct2.SetColorEntry(iColor, Tuple);
end;
err := DstBand.SetColorTable(ct2);
end;
DstBand.DataType := SrcBand.DataType;
// Copy band raster data.
for iLine := 0 to dsSrc.YSize – 1 do
begin
SrcBand.RasterIO(gdal.GF_Read, 0, iLine, dsSrc.XSize, 1, @Scanline[0], SrcBand.DataType);
DstBand.RasterIO(gdal.GF_Write, 0, iLine, dsSrc.XSize, 1, @Scanline[0], DstBand.DataType);
end;
end; //for iBand
Memo1.Lines.Add(‘Copy seems to have completed.’);
finally
SrcBand.Free;
DstBand.Free;
Drv.Free;
dsDst.Free;
dsSrc.Free;
end;
end;

function TfrmMain.ReportCorner(CornerName : string; pixel, line : double;
gt : TDynArrayDouble;
ct : TOGRCoordinateTransformation) : string;
var
geox, geoy, longitude, latitude, Z : double;
latlong_valid : boolean;
begin
geox := gt[0] + pixel * gt[1] + line * gt[2];
geoy := gt[3] + pixel * gt[4] + line * gt[5];

latlong_valid := false;

if (ct.IsValid) then
begin
Z := 0;
longitude := geox;
latitude := geoy;
latlong_valid := ct.TransformEx(longitude, latitude, Z);
end;

if (latlong_valid) then
result := CornerName + floattostr(geox) + ‘,’ + floattostr(geoy)
+ ‘ ‘ + floattostr(longitude) + ‘,’ + floattostr(latitude)
else
result := CornerName + floattostr(geox) + ‘,’ + floattostr(geoy);
end;

initialization
{$I main.lrs}
end.
[/sourcecode]

  • จากโคีดจะเห็นว่าเรียกใช้คลาส TGDALDriverManager ผ่านตัวแปร drvm แล้วทำการ register ฟอร์แม็ตที่สนับสนุนก่อน ด้วยคำสั่ง ที่อยู่ในไฟล์ gdalcore.pas คือ gdalcore.AllRegister แล้วเรียกใช้ property ของ drvm เพื่อ iterate ดูฟอร์แม็ตแต่ละอย่าง แล้วตรวจสอบด้วยว่าแต่ละฟอร์แม็ตสนับสนุนทั้งอ่านและเขียน หรืออ่านอย่างเดียว ลองศึกษาโค๊ดดูครับไม่มีอะไรยาก

การแสดงระบบพิกัดของไฟล์ Raster

  • มาดูระบบพิกัดของไฟล์ georeference กันคลิกที่เมนู File > Open ลองเปิดไฟล์ Test_dem.tif ดู แล้วคลิกที่เมนู Tools > Coordinates System Information จะเห็นดังรูปด้านล่าง ระบบพิกัดที่แสดงจะอยู่ในรูปที่เรียกว่า Well Known Text

แสดงระบบพิกัดของไฟล์ Raster

  • นอกเหนือจากแสดงระบบพิกัดแล้วยังแสดงค่าพิกัดของรูป มุมบนซ้าย มุมบนขวา มุมล่างซ้ายและมุมล่างขวา เนื่องจากระบบพิกัดของไฟล์รูปเป็น UTM (Projection เป็น Transverse Mercator ค่า Central Meridain = 99) จึงแสดงเป็นค่าพิกัด UTM

การแสดงฟอร์แม็ตที่ GDAL สนับสนุน

  • ที่เมนูของโปรแกรมคลิกที่ Tools > List Drivers จะเห็นจำนวนฟอร์แม็ตที่ GDAL สนับสนุนทั้งหมด 75 ฟอร์แม็ต สังเกตุเห็นว่าคำย่อคำหน้าเช่น GTiff คือฟอร์แม็ต GeoTiff (*.tif), HFA คือ Erdas Imagine (*.img), JPEG คือ Jpeg (*.jpg) เป็นต้น

แสดงฟอร์แม็ตที่ GDAL สนับสนุน

อ่านทุดสอบไฟล์ Raster

  • ที่เมนูคลิกที่ File > Open เลือกไฟล์ Test_image.tif จากนั้นคลิกที่เมนู Tools > Read Test จะเห็นโปรแกรมแสดงผลดังรูปด้านล่าง

แสดง Metadata จากการอ่านทดสอบไฟล์รูป

การแปลงฟอร์แม็ต (CreateCopy)

  • คำสั่ง CreateCopy ผ่านฟังก์ชั่น GDALCreateCopy ของ TGDALDriver ใช้งานได้ง่ายผ่านคำสั่งไม่กี่คำสั่ง สำหรับข้อจำกัดบางอย่างสามารถอ่านได้ที่ www.gdal.org
  • ที่เมนูของโปรแกรมคลิก File > Open แล้วเลือกไฟล์ Test_Dem.tif คลิกที่เมนู Tools > CreateCopy Test เราจะลองแปลงเป็นฟอร์แม็ตของ Erdas Imagine(*.img) ถ้าดูในโค๊ดจะเห็นชื่อย่อของฟอร์แม็ตคือ HFA เมื่อโปรแกรมถามชื่อไฟล์ที่ต้องการเซฟหรือแปลง ป้อนชื่อเป็น Test.img

 

 

  • ลองใช้โปรแกรม ER Viewer เปิดไฟล์ Test.img ว่าเขียน DEM ได้ถูกต้องไหม

แปลงฟอร์แม็ตของ DEM จาก GeoTiff เป็น Erdas Imagine

การ Copy ไฟล์ใหม่แบบมีเงื่อนไข

  • การ copy ไฟล์ ผ่านฟังก์ชั่น GDALCreate ของ TGDALDriver มีทางเลือกมากกว่า GDALCreateCopy สามารถกำหนดขนาดของรูปได้ จำนวน Band และกำหนดขนาดของข้อมูลในแต่ละ pixel ได้ แต่เท่าที่ผมลองแล้วไม่ work ถ้าเป็นฟอร์แม็ตต่างกัน ผมไม่แน่ใจว่ารันใน Lazarus อาจจะไม่ work ลองไปรันโค๊ดต้นฉบับ VB6 ก็ไม่ work เหมือนกัน อย่างไรก็ตามจะพยายามอ่าน manual ดูว่าติดปัญหาอะไร ในตอนนี้ลอง copy ไปไฟล์ใหม่ในฟอร์แม็ตเดิมๆ
  • ที่เมนูคลิกที่ File > Open เลือกไฟล์ Test_Image.tif  แล้วคลิกที่เมนู Tools > Create Test ป้อนชื่อไฟล์เป็น out_create.tif ดูรูปด้านล่าง

แสดงการ copy ไฟล์

  • จากรูปด้านบนในโปรแกรมต้นฉบับ เมื่อ copy แล้วโปรแกรมจะแจ้งว่า “Copy seesm to have completed” หมายความว่าดูเหมือนจะสำเร็จ แสดงว่ามีอะไรบางอย่างที่หวังผลไมได้ 100% ผมลองเปิดดูไฟล์ out_create.tif ด้วย ER Viewer ปรากฎว่าใช้ได้
  • ต้องบอกกันนิดว่าโค๊ดบางส่วนยังไม่ clean นักอาจจะมี bug เพราะยังเป็นโค๊ดแปลงรุ่นแรกๆครับ