การเขียนโปรแกรมคำนวณการแปลงค่าพิกัดระหว่าง UTM Grid และ Geographic (Lat/Long) ด้วย Lazarus และ Delphi (ตอนที่ 1)

ตัวอย่างโปรแกรมแปลงค่าพิกัด GeoCalc และ GeoTrans

  • ผู้ใช้งานด้าน GIS และคนที่ทำงานด้านสำรวจ (Surveying) คงหลีกหนีการแปลงพิกัดระหว่าง UTM Grid และ Geographic (Latitude, Longitude) ได้ยาก ทั้งทางตรงและทางอ้อม ทางตรงได้แก่ ใช้โปรแกรมแปลงพิกัดเช่น GeoCalc (ฟรี แต่ไม่ Opensource) สามารถดาวน์โหลดได้ที่ http://www.geocomp.com.au/geocalc/gcalc420.exe โปรแกรมนี้สามารถคำนวณค่าพิกัดข้าม Datum ไปอีก Datum ได้ (แต่ต้องระวังอีกนิดเพราะโปรแกรม Geocalc มี mistake อยู่ค่า Semi-major axis (a) ของทรงรี Everest 1830 ผิดต้องแก้ไขจาก 6377267.345 เป็น 6377276.345)
  • อีกโปรแกรมคือ GeoTrans ฟรีและ Opensource มี source code ทั้ง Java และ C สามารถ download ได้ที่ http://earth-info.nga.mil/GandG/geotrans/geotrans3.7/master.zip มีให้เลือก download ทั้ง Windows และ Unix (ผมยังไม่เคยนำไปใช้ใน Linux)  GeoTrans เป็นโปรแกรมที่สามารถใช้แปลงค่าพิกัดระหว่าง Datum สนับสนุน Datum & Map Projection ได้หมดทั่วโลก ตัว source code มีคนนำไปแปลงและใช้เช่น แปลงเป็น Delphi ใช้ในโปรแกรม TatukGIS (ไม่ฟรี) ซึ่งเป็น GIS Development Tool Kit มี ActiveX, .NET Winform และ Native VCL ให้เลือกใช้และใช้ได้บนวินโดส์เท่านั้น
  • ส่วนการใช้โปรแกรมแปลงค่าพิกัดที่เป็น Tool ติดมากับโปรแกรมอื่นๆ เช่นบน ERDAS Imagine, Terramodel,  Global Mapper หรือกระทั่งบน Web ก็มีให้เห็นมากพอสมควร

Library สำเร็จรูป “GeographicLib”

  • ถ้าเข้าไปดู source code ของ GeoTrans จะเห็นสูตรที่ใช้แปลงค่าพิกัดเหล่านี้ ยาวเหยียด แยกตามเส้นโครงแผนที่ (Map Projection) ถ้าเกิดจะมีใครพัฒนาโปรแกรมที่เกี่ยวข้อง GIS จาก Ground to top คงต้องศึกษา Library ของ GeoTrans ถ้าไม่คิดจะใช้ Library ของผู้อื่น พูดถึง Library ที่เป็น Opensource ก็ต้อง GeographicLib เขียนโดย Charles Karney ดาวน์โหลดได้ที่ Geographic.zip ตัว GeographicLib เขียนด้วย C++ ตัวนี้ปรับปรุงมาจาก GeoTrans  แต่มีฟังก์ชั่นมากกว่า GeoTrans เสียอีก เช่นการคำนวณเรื่อง Geoid หันมาดู source code ที่เขียนด้วย Delphi ผมหาไม่พบ ต้องยอมรับว่าเปอร์เซนต์ผู้ใช้ Delphi นั้นน้อยน่าใจหาย หลังจาก Borland ล้มหายตายจาก ผมไม่ทราบว่า Developer ของ Delphi จะมีเหลือมากน้อยเท่าใด ส่วน Lazarus ไม่ต้องพูดถึงว่าน้อยขนาดไหน


Datum Transformation Diagram
Datum Transformation Diagram
  • จากไดอะแกรมรูปข้างบนการแปลงค่าพิกัดกริดเช่นจาก UTM Grid (N,E,Z) บน WGS84 => UTM Grid (N,E,Z) บน Indian 1975 datum ได้สองเส้นทางหลัก

    • เส้นทางที่ 1 WGS84( N,E,Z)  ==> WGS84(Lat,Long,Height)  == Molodensky Transformation ==> Indian1975(Lat,Long,Height) ==> Indian1975(N,E,Z)
    • เส้นทางที่ 2 WGS84( N,E,Z)  ==> WGS84(Lat,Long,Height) ==>WGS84(X,Y,Z) ==7,5,3 Parameter Transformation==> Indian1975(X,Y,Z) ==> Indian1975(Lat,Long,Height) ==> Indian1975(N,E,Z)
  • เส้นทางที่ 1 จะสั้นกว่า จากค่าพิกัด Latitude,Long itude คำนวณโดยใช้ Molodensky Transformation หาค่าพิกัด Latitude,Longitude บนอีก datum ได้ ส่วนเส้นทางที่ 2 จะยาวกว่า จากค่าพิกัด Latitude,Longitude จะคำนวณไปหาค่าพิกัดในระบบ Cartesian แล้วใช้ค่าพารามิเตอร์ 3,5 หรือ 7 ค่าในการ Transformation ไปยังระบบ Cartesian บนอีก datum แล้วคำนวณขึ้นไปหา Latitude,Longitude จนถึงค่าพิกัดกริดยูทีเอ็ม ถ้าต้องการ
  • ซึ่งถ้าจะโปรแกรมกันจริงๆ ตามไดอะแกรมด้านบน ก็ใหญ่ขนาด GeographicLib ถ้าจะต้อง support ทุก datum ทุก Map Projection ที่กล่าวไปแล้ว แต่ในตอนนี้ผมขอบีบขอบเขตให้เล็กลง เราจะทำการคำนวณแปลงค่าพิกัดระหว่าง Lat/Long และ UTM Grid บน datum เดียวกันเท่านั้น เช่นคำนวณ Latitude,Longitude ของ WGS 84 ไปเป็นค่าพิกัดยูทีเอ็มกริด (Northing, Easting) บน WGS84 หรือในทางกลับกัน

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

  • ขอนำกลับมาที่ datum ก่อนจะไปต่อ datum เป็นพื้นหลักฐานที่แต่ละประเทศต่างก็ใช้แตกต่างกันไป เช่นประเทศในย่านเอเชียตะวันออกเฉียงใต้เช่นพม่า ไทย เขมร ลาว ใช้ Indian 1975 datum มี Spheroid (Ellipsoid) ทรงรีที่เราเรียกว่า Everest 1830 ทรงรีที่เราใช้แทนสัณฐานของโลกจะยุบที่ขั้วโลกเหนือใต้ โป่งออกตรงเส้นศูนย์สูตร สัณฐานของทรงรีจะกำหนดด้วยค่าSemi- major Axis (a) ความยาวของแกนหลัก ค่าการยุบตัว Flattening (f) และความยาวแกนรอง Semi-minor Axis (b) ค่า 3 ค่านี้จะสัมพันธ์กันหมดคือ b = a(1-f) สำหรับ Everest 1830 มีค่า a=6377276.345 , f = 1/300.8017
  • มาถึงยุคของ GPS พื้นหลักฐานของ GPS ในปัจจุบันคือ WGS84 ที่หลายประเทศรวมทั้งไทยเรา หันมาใช้ datum นี้กันมากขึ้น สัณฐานทรงรี WGS94 ค่า a=6378137 ค่า f=1/298.257 223 563

 รูปทรงของทรงรี (ภาพจาก Wikipedia)

  • ทรงรีที่เราใช้กันในโลกนี้มีไม่มาก ไม่น่าปวดหัว มีตัวหลักๆประมาณ 10 กว่าทรงรี
  • ถ้าจะเขียนโปรแกรมให้ใช้ได้ World wide ทั่วโลก สิ่งที่จะปวดหัวก็คือเส้นโครงแผนที่ Map Projection และ Transformation Parameters เพราะมีมากแต่ละประเทศจะใช้ค่าพารามิเตอร์จำนวน 3, 5, 7 เพื่อ Transformation ค่าพิกัดแตกต่างกันไป ตัวอย่างเช่นประเทศแถวตะวันออกกลาง จะใช้ Ellipsoid “International 1924” ใช้เส้นโครงแผนที่ UTM เหมือนกัน แต่จะต่างกันตรงที่ใช้ 3 Parameter Transformation แตกต่างกันไป แต่ถ้าใช้ในประเทศไทยเรา เส้นโครงแผนที่ที่ใช้ก็มีเพียงแต่ Transverse Mercator (TM) ซึ่งเส้นโครงแผนที่ TM ถ้าแบ่งโซนกันแน่นอนก็คือระบบพิกัด Universal Transverse Mercator (UTM) ที่เราคุ้นเคยกันอยู่

เริ่มสร้าง New Project ด้วย Lazarus

  • เปิดโปรแกรม Lazarus คลิกที่เมนู Project > New Project… จะเห็น Dialog เลือก Application คลิกที่ OK
สร้าง New Project
สร้าง New Project
  • จะเห็นฟอร์มเปล่าๆ ดังรูปข้างล่าง ในตอนนี้ผมรัน Lazarus บนวินโดส์
Lazarus รันบนวินโดส์
สร้างโปรเจคใหม่ (Lazarus บนวินโดส์)
  • ก่อนจะไปต่อตอนที่ 2 ผมขอแสดงโปรแกรมขนาดเล็ก ที่เราจะเขียน เพื่อคำนวณค่าพิกัดระหว่าง UTM และ Geographic (Lat/Long) ดังรูปด้านล่าง
โปรแกรมแปลงพิกัดระหว่าง UTM และ Geographic (Lat/Long)
โปรแกรมแปลงพิกัดระหว่าง UTM และ Geographic (Lat/Long)
  • ทิ้งท้ายไว้ตรงนี้ครับ ติดตามตอนต่อไป

Leave a Reply

Your email address will not be published.