ไม่นานมานี้มีผมดาวน์โหลดไฟล์โปรแกรมและข้อมูลของ TGM2017 เรียกเต็มๆคือ Thailand Geoid Model 2017 ที่เป็นโครงการร่วมมือจากหลายๆฝ่ายของทางราชการ ผมยังไม่มีโอกาสได้นำไปใช้งาน โดยเฉพาะจะนำมาประยุกต์ใช้ในงานรังวัด GNSS ถือว่าเป็นสิ่งที่พวกเรารอคอยมานานที่จะได้มี local geoid model มาใช้งานกัน
โดยเฉพาะงานรังวัด GNSS เมื่อคำนวณแล้วจะได้ค่าพิกัดทางราบ และทางดิ่งจะได้ความสูง (h – Ellipsoid height) ที่เทียบกับทรงรี WGS84 ถ้าเรามีแบบจำลองความสูงจีออยด์สามารถคำนวณหาความสูง (H – Orthometric height) หรือระดับน้ำทะเลปานกลางได้ จากสูตร h = H + N โดยที่ไม่ต้องมีการเดินระดับไปหาจุดที่รังวัด GNSS แต่อย่างใด โดยที่โปรแกรมคำนวณหาค่าความสูงจีออยด์ (N – Geoid Height หรือ Geoid Separation) ที่ผมจะมาลองเขียนโปรแกรมเล็กๆเพื่อหาค่า N นี้ด้วยภาษาไพทอนในลำดับต่อไป
ในตอนนี้ผมจะใช้ไลบรารี Proj.4 เรียกไลบรารีในเวอร์ชั่นภาษาไพทอนนี้ว่า PyProj ตัวไลบรารี Proj.4 เองก็เพิ่งสนับสนุนการคำนวณ Vertical Datum ไปไม่นานนี้เอง สมมติว่าถ้าไม่ได้ใช้ไลบรารีตัวนี้เป็นตัวช่วยคงเป็นงานหนักหนาสาหัสเหมือนกัน
แปลงรูปแบบไฟล์ของแบบจำลองความสูงจีออยด์ TGM2017
ไฟล์แบบจำลองความสูงจีออยด์ TGM2017 ที่ผมดาวน์โหลดได้มานั้นเป็นแอสกี้ (Text file) ขนาดประมาณ 10.7 MB และมีโปรแกรมคำนวณมาด้วยชื่อไฟล์ execute คือ “TGM2017.EXE” สำหรับไฟล์แอสกี้ในคู่มือแสดงรายละเอียดดังนี้
- ค่าละติจูด (Latitude) เริ่มต้น : 3 องศา
- ค่าลองจิจูด (Longitude) เริ่มต้น : 95 องศา
- ความละเอียดเชิงพื้นที่ : 1 ลิปดา
- จานวนหลัก (Column) : 780 หลัก
- จานวนแถว (Row) : 1200 แถว
สำหรับไลบรารี PyProj นั้นถ้าต้องการคำนวณ Vertical Datum ต้องการไฟล์ในรูปแบบ gtx ที่กำหนดมาตรฐานโดย NOAA เพื่อให้การใช้งานมีประสิทธิภาพมาก ทางนี้แนะนำให้ใช้ไฟล์ไบนารีเพื่อประหยัดเมมโมรี สำหรับรูปแบบไฟล์ gtx นั้นไม่มีอะไรมาก กำหนดให้มีไฟล์ header เพื่อบอกจุดเริ่มต้นขอบเขตการใช้งาน ระยะห่างของละติจูดและลองจิจูดของแต่ละจุดในแนวตั้งและแนวนอน จำนวนแถวและจำนวนคอลัมน์ ที่เหลือจะเป็นความสูงจีออยด์ในแต่ละจุดมีขนาดเท่ากับ จำนวณแถวคูณจำนวนคอลัมน์ ตัวหัวไฟล์ผมสามารถเขียนได้ดังนี้
3.0 95.0 0.01666666666 0.01666666666 1200 780
อธิบาย : ละติจูดของมุมล่างซ้าย ลองจิจูดของมุมล่างซ้าย ระยะห่างของจุดในแนวละติจูด ระยะห่างของจุดในแนวลองจิจูด จำนวนแถว จำนวนคอลัมน์
ผมตรวจสอบไฟล์แอสกี้แล้วพบว่า ความสูงจีออยด์ถูกปูเรียงมาแล้วตามข้อกำหนดของ NOAA ดังนั้นผมจะอ่านไฟล์แอสกี้นี้แล้วแปลงเป็นไบนารี เพียงแต่เพิ่มส่วนหัวเข้าไป
ตัวอย่างโปรแกรมไพทอนสำหรับแปลงไฟล์แบบจำลองความสูงจีออยด์
ผมเขียนสคริปต์ภาษาไพทอน เผื่อมีใครสนใจ ผมตึ้งชื่อไฟล์ว่า “rtsd2gtx.py” โปรแกรมมีขนาดเล็กมากๆ ไม่กี่บรรทัด หลักการคืออ่านไฟล์แอสกี้ “tgm2017.asc” มาแล้วเขียนไฟล์ไบนารีในรูปแบบ gtx ข้อควรระวังคือไฟล์ไบนารี ต้องเป็นแบบ Big Endian ในโมดูล struct ผมจึงใส่เครื่องหมาย “>” ไปด้วย
ตัวอย่างเขียนขนิดข้อมูลเป็น double ใช้คำสั่ง f.write(struct.pack(“>d”, 3.0))
เขียนข้อมูลเป็น integer ใช้คำสั่ง f.write(struct.pack(“>i”, 1200))
<import struct
def readRTSDAsciiFile():
geoid = []
with open('tgm2017.asc', 'r') as f:
geoid = [[float(n) for n in line.split()] for line in f] #2D list, one list contain one row.
#TGM2017 was organized already from lowest row (min latitude) to highest row (max latitude).
#No need to reverse.
return geoid
def writeGTXBinaryFile(geoid):
with open('tgm2017.gtx', 'wb') as f:
# pack with ">" for big endian.
# write Header of file
f.write(struct.pack(">d", 3.0)) #lower left Latitude
f.write(struct.pack(">d", 95.0)) #lower left Longtitude
f.write(struct.pack(">d", 1.0/60.0)) #delta Latitude = 1 second
f.write(struct.pack(">d", 1.0/60.0)) #delta Longitude = 1 second
f.write(struct.pack(">i", 1200)) #number of rows
f.write(struct.pack(">i", 780)) #number of columns
# write height values.
for row in geoid:
for c in row:
f.write(struct.pack(">f", c))
# main program
g = readRTSDAsciiFile()
writeGTXBinaryFile(g)
ณ ตอนนี้ผมสร้างโฟลเดอร์ เอาไฟล์ tgm2017.asc มาวางไว้ และมีโปรแกรมสคริปต์ “rtsd2gtx.py” อยู่ในโฟลเดอร์เดียวกัน
จากนั้นทำการรันโปรแกรมด้วยคำสั่งที่คอมมานด์ไลน์ python rtsd2gtx.py รอแป๊ปหนึ่งไฟล์ gtx ที่เก็บแบบจำลองความสูงจีออยด์จะถูกสร้าง และมีขนาดเล็กลงจากไฟล์แอสกี้ดั้งเดิมประมาณ 3 เท่า
ตัวอย่างโปรแกรมไพทอนสำหรับคำนวณหาความสูงจีออยด์
ต่อไปโปรแกรมไพทอนที่จะนำมาทดสอบคำนวณหาความสูงจีออยด์ ผมตั้งชื่อง่ายๆว่า “test.py” ถ้าจะลองรันโปรแกรมนี้เพื่อทดสอบต้องติดตั้งไลบรารี pyproj ให้เรียบร้อยก่อน
โค้ดโปรแกรมเริ่มต้นจากสร้าง projection ต้นทาง (ตัวแปร wgs84 egm) โดยที่กำหนด +geoidgrids=tgm2017.gtx และสร้าง projection ปลายทาง (ตัวแปร wgs84) โดยที่เว้นไม่ใส่แบบจำลองความสูงจีออยด์
import pyproj # Import the pyproj module
# Define source coordinate system with TGM2017.
wgs84egm = pyproj.Proj("+init=EPSG:4326 +geoidgrids=tgm2017.gtx")
# Define target coordinate system without geoid model.
wgs84 = pyproj.Proj("+init=EPSG:4326")
lat, lon = 18.3353398540, 99.371215293000
x, y, z = pyproj.transform(wgs84egm, wgs84, lon, lat, 0.0)
print("1) Latitude = %.8f, Longitude = %.8f, Geoid Height = %.5f" % (y, x, z))
lat, lon = 7.75906518100, 98.303594340000
x, y, z = pyproj.transform(wgs84egm, wgs84, lon, lat, 0.0)
print("2) Latitude = %.8f, Longitude = %.8f, Geoid Height = %.5f" % (y, x, z))
lat, lon = 9.18558870000, 99.843706600000
x, y, z = pyproj.transform(wgs84egm, wgs84, lon, lat, 0.0)
print("3) Latitude = %.8f, Longitude = %.8f, Geoid Height = %.5f" % (y, x, z))
lat, lon = 6.96243169, 99.77398865
x, y, z = pyproj.transform(wgs84egm, wgs84, lon, lat, 0.0)
print("4) Latitude = %.8f, Longitude = %.8f, Geoid Height = %.5f" % (y, x, z))
lat, lon = 7.20852844, 99.72353542
x, y, z = pyproj.transform(wgs84egm, wgs84, lon, lat, 0.0)
print("5) Latitude = %.8f, Longitude = %.8f, Geoid Height = %.5f" % (y, x, z))
lat, lon = 7.35631365, 100.12797020
x, y, z = pyproj.transform(wgs84egm, wgs84, lon, lat, 0.0)
print("6) Latitude = %.8f, Longitude = %.8f, Geoid Height = %.5f" % (y, x, z))
lat, lon = 7.23806984, 100.55255020
x, y, z = pyproj.transform(wgs84egm, wgs84, lon, lat, 0.0)
print("7) Latitude = %.8f, Longitude = %.8f, Geoid Height = %.5f" % (y, x, z))
lat, lon = 7.09316728, 100.56368280
x, y, z = pyproj.transform(wgs84egm, wgs84, lon, lat, 0.0)
print("8) Latitude = %.8f, Longitude = %.8f, Geoid Height = %.5f" % (y, x, z))
lat, lon = 6.99643503, 100.43034140
x, y, z = pyproj.transform(wgs84egm, wgs84, lon, lat, 0.0)
print("9) Latitude = %.8f, Longitude = %.8f, Geoid Height = %.5f" % (y, x, z))
lat, lon = 13.0+42.0/60+38.64/3600, 100+29.0/60+13.02/3600
x, y, z = pyproj.transform(wgs84egm, wgs84, lon, lat, 0.0)
print("10) Latitude = %.8f, Longitude = %.8f, Geoid Height = %.5f" % (y, x, z))
เมื่อจัดเก็บไฟล์เรียบร้อย ในตอนนี้ ถ้าดูที่โฟล์เดอร์จะเห็นไฟล์ดังนี้
จุดทดสอบ
ถ้าดูในโค้ดจะเห็นจุดทดสอบทั้งหมด 10 จุด โดยที่ 9 จุดแรกผมนำมาจากไฟล์ INPUT.DAT ที่ติดมากับโปรแกรม TGM2017.EXE ถ้าเปิดดูด้วย Notepad จะเห็นดังนี้
ส่วนจุดที่ 10 ได้มาจากสไลด์ ชื่อ BM11SP Latitude = 13° 42′ 38.64″ N Longitude = 100° 29′ 13.02″E
คำนวณ
จากนั้นทำการรันโปรแกรมบนคอมมานด์ไลน์ ใช้้คำสั่ง python test.py จะได้ผลลัพธ์ดังนี้
ผมนำผลลัพธ์มาเปรียบเทียบกับที่คำนวณด้วย TGM2017.EXE พบว่ามีความต่างกันเล็กน้อยที่หลักมิลลิเมตร สาเหตุน่าจะเกิดตอน interpolation โดยการนำจุดรอบๆจากแบบจำลองความสูงจีออยด์ ซึ่งโปรแกรมอาจจะเลือกใช้แบบ BiCubic หรือ BiLinear ผมไม่แน่ใจว่า Proj.4 ใช้แบบไหน ตอนนี้ยังไม่มีโอกาสไปดูโค้ดต้นฉบับ
การพัฒนาโปรแกรมสำหรับการใช้งานในอนาคต
ผมวางแผนว่าจะนำแบบจำลองความสูงจีออยด์ TGM2017 ไปคำนวณในโปรแกรมรวมเครื่องมือฉบับกระเปาสำหรับช่างสำรวจ Surveyor Pocket Tools ได้แก่การคำนวณความสูงจีออยด์ การคำนวณหาสเกลแฟคเตอร์
สรุป
ขอขอบคุณคณะผู้จัดทำโครงการนี้ทุกๆท่าน สำหรับการจัดทำแบบจำลองความสูงจีออยด์ความละเอียดสูงสำหรับประเทศไทย อันจะเป็นประโยชน์เอนกอนันต์ต่อการประยุกต์ใช้ในแวดวงงานด้านสำรวจและด้านอื่นๆอีกมากในปัจจุบันและอนาคตอันใกล้นี้ ผมเองยังมีความหวังลึกๆว่าน่าจะมีการพัฒนาปรับปรุงให้แบบจำลองนี้ให้มึความถูกต้องมากๆยิ่งขึ้นในอนาคต
สามารถกดคีย์ F1 – Dn (Down) เพื่อเลื่อนภาพลง) F2 – Up (เลื่อนภาพขึ้น) F3 – Lt (Left เลื่อนภาพไปด้านซ้าย) F4 – Rt (เลื่อนภาพไปด้านขวา) หรือกด F5 – M-> (เลือกเมนูด้านขวา) เมนูด้านขว่าจะมี F1 -Z+ (เพื่อขยายภาพ) F2 – Z- (เพื่อย่อภาพ) เสร็จแล้วกดคีย์ F6 -Done เพื่อออก กดคีย์ F6 – Done อีกครั้งเพื่อออกมาเมนูหลัก
สังเกตว่าค่าพิกัดที่ได้
จะได้ค่าพิกัดผลลัพธ์ออกมาจากวิธีการของ Collins (1167.446E, 2016.916N) ซึ่งเท่ากับวิธีการของ Tienstra, Cassini
เช่นเดียวกันจะได้ค่าพิกัดผลลัพธ์ออกมาจากวิธีการของ Font-Llagunes (1167.446E, 2016.916N) ซึ่งเท่ากับวิธีการของ Tienstra, Cassini และ Collins แต่บางครั้งจะต่างกับวีธีการของ Kaestner-Burkhardt ในระดับหลักเมตร
ก็ได้รูปการวางตัวของจุดหลักฐาน 3 จุดและจุดที่ต้องการทราบค่าคือจุด P น่าเสียดายที่ผมไม่สามารถวาดมุมแสดง dimension เพราะข้อจำกัดของเครื่องคิดเลข และต้องเตือนกันอีกนิดว่าจุด 4 จุดนี้ต้องไม่อยู่บนวงกลมเดียวกัน (โอกาสเกิดน้อยมากๆครับ) ไม่ว่าสูตรไหนก็ไม่สามารถคำนวณได้ เรียกภาวะอย่างนี้ว่าเกิดเหตุการณ์สภาวะเอกฐาน Singularity
ค่า k คือ grid scale factor ค่า k นี้เราสามารถหาได้จาก สูตรด้านล่าง (เครดิตจาก 
แทนค่า
ก็ขอตั้งชื่อโปรแกรมเล็กๆสำหรับแปลงพิกัดบนพื้นหลักฐาน WGS84 นี้ว่า “GeoCalc Extra” ก่อนอื่นสูตรที่ใช้ในการคำนวณผมใช้ไลบรารีชื่อ 
เสร็จแล้วกดคีย์ F1 – Calc เพื่อคำนวณสเกลแฟคเตอร์ของจุดนี้
