ย้อนรอยวิธีสร้างไฟล์รูปแบบ PGM ของ TGM2017 สำหรับใช้ใน GeographicLib

ผมได้เขียนบล็อกเกี่ยวกับ TGM2017 (Thailand Geoid Model 2017) มาหลายตอนแล้ว ไม่นานมานี้ทางรุ่นพี่ที่เคารพอาจารย์ดร.ไพศาล สันติธรรมนนท์ ได้วานให้ตรวจสอบไฟล์ TGM2017-1.PGM ที่ทางอาจารย์ได้สร้างไว้ด้วยโค้ดไพทอนเพื่อนำมาใช้ในไลบรารี GeographicLib ผมทดสอบแล้วใช้งานได้ดี ในขณะเดียวกันผมเห็นว่าน่าสนใจเพราะสามารถเผยแพร่การใช้งาน TGM2017 ให้ใช้งานได้หลากหลายในวงกว้างยิ่งๆขึ้นไป ผมขอสรุปรูปแบบการใช้งานดังนี้

1. รูปแบบแอสกี้: TGM2017.ASC เป็นไฟล์รูปแบบ ASCII คำนวณโดย TGM2017.EXE ที่ทางผู้จัดทำแบบจำลอง TGM2017 ดร.พุทธิพลและทีมงาน ได้นำมาเผยแพร่ และผมขอขอบคุณทีมงานมา ณ ที่นี้
2. รูปแบบ NOAA (Datum Transformation Grid): TGM2017.GTX ที่ผมได้จัดทำไว้ สามารถใช้งานผ่านโปรแกรม Surveyor Pocket Tools ได้
3. รูปแบบ Trimble GGF (Geoid Grid File Format): TGM2017.GGF ที่ผมได้จัดทำไว้เช่นเดียวกัน สามารถใช้งานได้ในโปรแกรมตระกูลของ Trimble เช่น Trimble Business Center, Trimble Geo Office
4. รูปแบบ PGM (Portable Graymap): TGM2017.PGM ดร.ไพศาล สันติธรรมนนท์ ได้จัดทำไว้เพื่อใช้ในไลบรารี GeographicLib ขอขอบคุณอ.ดร.ไพศาล มา ณ ที่นี้ด้วยครับ

รูปแบบและข้อกำหนดของไฟล์ Geoid แบบ PGM

ผมว่ารูปแบบ PGM นั้นสร้างยากที่สุดเพราะต้องเอาข้อมูลแอสกี้ TGM2017 (1′ x 1′) ของดร.พุทธิพลและทีมงาน ขนาด 780 x 1200 มาปูลงทับ EGM2008-1.PGM (1′ x 1′) (ดาวน์โหลดไฟล์นี้ไว้ด้วยเพราะจะนำมาใช้ภายหลัง) ของไลบรารี GeographicLib ที่มีขนาดเต็ม 21600 x 10801 สองแบบจำลองนี้แต่ละช่องกริดห่างกัน 1 ลิปดาเหมือนกันจึงสามารถนำมาปูซ้อนทับกันได้ ตำแหน่งเริ่มต้นของ TGM2017 คือ Latitude, Longitude = 22.98333333. 95.0 ตรงกับ row, column = 5700, 4021 ไปสิ้นสุดที่ Latitude, Longitude = 3.0, 107.9833333 ตรงกับrow, column = 6479, 5220 หมายเหตุอย่าลืมครับว่าหนึ่งช่องนั้นกว้าง 1 ลิปดา

เอาละครับสิบปากว่าไม่เท่าตาเห็น ดูรูปด้านล่าง ให้จินตนาการให้นึกถึงตาราง excel ที่จำนวนคอลัมน์มีทั้งหมด 21600 ช่องและมีจำนวนแถวทั้งหมด 10801 แต่ละช่องกว้าง 1 ลิปดา ข้อมูลความสูงจีออยด์ (Geoid Height) จะอยู่ในช่องตารางนั้นแตกต่างกันไป จุดเริ่มต้นบนซ้ายสุดคือ latitude, longitude = 90N, 0E จุดสิ้นสุดด้านล่างซ้ายคือ latitude, longitude = 90S, 360E หรือ 90S, 0E (ย้อนกลับมาเส้นลองจิจูดเดิมเริ่มต้น)

กรอบใหญ่คือ EGM2008-1 ที่มีขนาดใหญ่เต็ม 21600 x 10801 แทนด้วยพื้นที่ทั้งโลก ถูกปูทับด้วย TGM2017 ขนาด 780 x 1200 แทนด้วยสี่เหลียมสีฟ้า

ก่อนจะไปต่อผมมาเกริ่นถึงรูปแบบ PGM นิดหนึ่ง ว่าเป็นรูปแบบการเก็บรูปภาพแบบ grayscale ที่มีมานานตั้งแต่ยุคคอมพิวเตอร์แรกๆ ลักษณะการเก็บแต่ละจุดภาพหรือ pixel ไล่ลำดับความหนักเบา แต่ละ pixel จะมีสีเทาแสดงด้วยตัวเลข 0 ถึงเลข 65535 ขึ้นอยู่กับความเข้ม แต่ละ pixel ใช้หน่วยความจำเท่ากับ 16 บิตหรือ 2 ไบต์ นั่นเอง การเก็บลักษณะแบบนี้ไม่ใช่เรื่องแปลกเพราะ DEM เราก็เก็บแบบนี้

การแทนที่ตัวเลขความสูง

ตัวเลขความสูงที่ค่าน้อยที่สุดของ EGM2008 (minimum value) = -106.910 m ส่วนความสูงที่มากที่สุด (maximum value) = 85.840 m ค่าต่ำสุดเลือกใช้ตัวเลขกลมๆคือ -108 m ตัวเลขนี้จะมาแทนด้วยสีเทาตัวเลข 0 ส่วนค่าสูงที่สุดคือ 108 m แทนด้วยสีเทาตัวเลข 65535 ตัวเลข -108 นี้จะเรียกว่า offset

ช่วงความสูงจีออยด์จากมากสุดลงมาหาต่ำสุดคือ 108+108 = 216 m แต่ช่วงของสีเทามีทัั้งหมด 65536 (นับเลข 0 ด้วย) ดังนั้น 216/65536 = 0.00329 เลือกใช้ 0.003 m ตัวเลขนี้เรียกว่า scale ทุกครั้งที่ตัวเลขเปลี่ยนตัวอย่าง 0 ไปหา 1 ความสูงจะเพิ่มทีละ 3 mm

ความสัมพันธ์คือ height = offset + scale. pixel (pixel คือ ตัวเลข 0-65535)

ลองคำนวณหา pixel จากความสูง -51.3520 m = (-51.352 + 108)/0.003 = 18832.667 ตัดทศนิยมทิ้งได้ 18832

เครื่องมือแก้ไขไฟล์ Hex Editor

ต่อไปจะมาย้อนรอยวิธีการสร้างไฟล์ PGM ด้วยโค้ดไพทอน ซึ่งผมใช้ทูลส์เครื่องมือแก้ไขไฟล์จำพวก Hex Editor มาช่วย ในกรณีที่สร้างมาแล้วสามารถมาเปิดดูค่าได้ตามตำแหน่งที่เราต้องการ ผมเลือกใช้ HxD เพราะฟรี ใช้งานง่ายและเก่งพอสมควร มาลองเปิดไฟล์ EGM2008-1.PGM ตั้ง byte order เป็น Big Endian ตั้ง byte per row = 32 ตั้ง byte group size = 2

ผมลองลาก header แล้ว copy มา paste ลง Notepad เพื่อให้อ่านง่ายจะเห็นตัวหัวไฟล์ดังนี้

สองไบต์แรกจะเป็น P5 อันเป็นการบอกว่านี่เป็นไฟล์ PGM แบบ 0-65535 gray scale และที่ขาดไม่ได้คือค่า Offset -108 และ Scale 0.003 จากนั้นตัวเลขความกว้างและความสูงของ EGM คือ 21600 และ 10801 ที่กล่าวไปแล้ว และบรรทัดสุดท้ายคือ 65535 คือตัวเลขสูงสุดที่เราใช้

เตรียมไฟล์ resource

ผมสร้างโฟลเดอร์ test โดยนำไฟล์ TGM2017.ASC และไฟล์ EGM2008-1.PGM มาไว้ดังรูปด้านล่าง

ผมเปลี่ยนชื่อไฟล์ EGM2008-1.PGM เป็นชื่อไฟล์ TGM2017-1.PGM จากนั้นผมเขียนโค้ดเพื่ออ่านไฟล์ TGM2017.ASC มาปูลงไฟล์ EGM2008-1.PGM ที่เราเปลี่ยนชื่อเป็น TGM2017-1.PGM เรียบร้อยแล้ว ตั้งชื่อไฟล์ของโค้ด rtsd2pgm.py

โค้ดไพทอน

กระบวนการทำงานคือ

1. เปิดอ่านไฟล์ TGM2017.ASC ด้วยฟังก์ชั่น readRTSDAsciiFile() จากนั้นกลับลำดับบรรทัดเพราะต้องการเรียงไฟล์ใหม่ให้ตรงกับรูปแบบ EGM2008
2. เปิดไฟล์ TGM2017-1.PGM (เดิมเป็น EGM2008-1.PGM) เขียน header ทับด้วย writePGMHeaderFile()
3. วนลูปสองลูป for ทำการแปลงค่าความสูงจีออยด์ที่อ่านจากไฟล์ในข้อ 1 ในแต่ละจุดเป็นเลข pixel 0-65535 เริ่มต้นที่คอลัมน์ 5700 แถว 4021สุดท้ายจะได้ไฟล์ TGM2017-1.PGM ที่พร้อมจะนำไปใช้งานได้

เปิด command prompt ของวินโดส์ ใช้คำสั่ง cd มาที่โฟลเดอร์ที่ติดตั้งโปรแกรมและไฟล์ resources ทำการรันด้วยคำสั่ง

python rtsd2pgm.py

ก็ลองดูโค้ดไพทอนของผมด้านล่าง โค้ดสั้นๆทำความเข้าใจได้ง่าย

import struct
HEADERSIZE = 404
W2008 = 21600  #derived from 180 degree x 2 x 60 second
ROW = 4021     #start row of TGM2017 based zero.
COL = 5700     #start column of TGM2017 based zero.

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 from lowest row (min latitude) to highest row (max latitude).
    #But EGM2008 needs organize from highest row (max latitude) to min row (min latitude).
    #Reordered line of TGM2017 to match EGM2008 requirement.	
    return reversed(geoid)
    
def writePGMHeaderFile(f):
    # write Header of file 404 bytes
    f.seek(0, 0)
    f.write("P5\n".encode("ascii")) #magic word is "P5" for indicating a Portable Gray Bitmap(PGM).
    f.write("# Geoid file in PGM format for the GeographicLib::Geoid class\n".encode("ascii"))
    f.write("# Description WGS84 TGM2017, 1-minute grid\n".encode("ascii"))
    f.write("# URL https://civil.eng.cmu.ac.th/people/puttipol-dumrongchai \n".encode("ascii"))
    f.write("# DateTime 2019-10-25 18:00:00\n".encode("ascii"))
    f.write("# MaxBilinearError 0.025\n".encode("ascii"))
    f.write("# RMSBilinearError 0.001\n".encode("ascii"))
    f.write("# MaxCubicError 0.003\n".encode("ascii"))
    f.write("# RMSCubicError 0.001\n".encode("ascii"))
    f.write("# Offset -108\n".encode("ascii"))
    f.write("# Scale 0.003\n".encode("ascii"))
    f.write("# Origin 90N 0E\n".encode("ascii"))
    f.write("# AREA_OR_POINT Point\n".encode("ascii"))
    f.write("# Vertical_Datum WGS84\n".encode("ascii"))
    f.write("21600  10801\n".encode("ascii"))
    f.write("65535\n".encode("ascii"))
    
def overWriteBlockofTGM2017toEGM2008BinaryFile(geoid):
    with open('tgm2017-1.pgm', 'r+b') as f:
        writePGMHeaderFile(f)
        r = ROW
        for row_tgm in geoid:
            offset = HEADERSIZE + 2*(W2008*r+(COL))
            r += 1
            c = 0
            for h in row_tgm:
                v = int((h + 108.0)/0.003) #convert to 0-65535 gray scale.
                f.seek(offset+c*2, 0)
                f.write(struct.pack(">H", v)) #H is unsigned short in C language.
                c += 1   
				
# Main program	
g = readRTSDAsciiFile()
overWriteBlockofTGM2017toEGM2008BinaryFile(g)

อาจจะมีคำถามว่าทำไมต้องอ่านไฟล์ tgm2017.asc แล้วไปปูทับ egm2008-1.pgm ให้เสียเวลาทำไม คำตอบคือเป็นข้อกำหนดของไลบรารี GeographicLib ( Models covering a limited area will need to be “inserted” into a world-wide grid before being accessible to the Geoid class )

ตรวจสอบไฟล์

ผมลองตรวจสอบไฟล์ TGM2017-1.PGM ที่สร้างขึ้นมาจากโค้ด จุดข้อมูลเริ่มต้นที่คอลัมน์ 5700 และคอลัมน์ 4021 สามารถคำนวณได้ดังนี้ 404 + 2 x (4021 x 21600 + 500) = 173719004 หมายเหตุ 404 คือขนาด header และคูนด้วยสอง เพราะหนึ่งจุดความสูงใช้พื้นที่ในการจัดเก็บ 2 ไบต์

เปิดไฟล์ด้วย HxD ใช้เมนู Search > Go to… ป้อนตัวเลขนี้ช่อง offset

จะเห็นค่าที่ตำแหน่งนี้เป็นเลขฐาน 16 คือ 49C2 เรารู้ว่าเป็นตัวเลข 16 บิตหรือ 2 ไบต์ คิดเป็นเลขฐานสิบคือ 18882 คำนวณกลับเป็นความสูงด้วยสูตร height = offset + scale x pixel = -108 + 0.003 x 18882 = -51.354 m

คำถามคือค่าความสูง -51.354 นี้อยู่ตรงไหนของ TGM2017.ASC คำตอบคือบรรทัดที่ 1200 บรรทัดสุดท้ายด้านซ้ายสุดคือตำแหน่งพิกัด 22.98333333N, 95E เริ่มต้นนั่นเอง แต่ค่าต่างกันนิดหน่อย 2 มม. นั่นเป็นเพราะระบบจัดเก็บเป็น 0-65535 pixel ที่ผมกล่าวไว้ข้างต้น อันจะทำให้ค่าความสูงขยับทีละ 3 มม.

การปูเรียงค่าความสูงของ TGM2017.ASC บรรทัดแรกจุดแรก จะเป็นค่าพิกัด 3N, 95E จากนั้นบรรทัดจะขยับไปทีละ 1 ลิปดา จากน้อยไปหามาก จนถึงบรรทัดที่ 1200 ที่จุดแรกจะเป็นค่าพิกัด 22.98333333N, 95E

ทดลองการใช้งาน GeographicLib

ความจริง GeographicLib นั้นมีหลายโมดูล ผมเองก็นำไปใช้ใน Surveyor Pocket Tools สำหรับคำนวณหา Geodesic distance บนทรงรี การคำนวณหา Geoid height นั้นเป็นส่วนหนึ่งเช่นเดียวกัน สามารถไปดาวน์โหลดไฟล์ไบนารีเพื่อมาทำการติดตั้งได้ที่หน้านี้ Using a binary installer ผมเลือกไฟล์ GeographicLib-1.50-win64.exe มาทำการติดตั้ง ตอนติดตั้งให้คลิกเลือกสร้าง path ไว้ให้ด้วยจะได้สะดวกในการเรียกใช้งาน ซึ่งโปรแกรมในชุดนี้เป็น command line ทั้งหมด

ดาวน์โหลดไฟล์ TGM2017-1.PGM

ไฟล์ที่ผมโฮสต์ไว้เป็นรุ่นที่จัดทำโดยดร.ไพศาล สันติธรรมนนท์ สามารถดาวน์โหลดได้ที่ TGM2017-1.PGM เมื่อดาวน์โหลดมาแล้วให้ copy ไปไว้ที่โฟลเดอร์ c:\programdata\geographiclib\geoids (ผมติดตั้ง GeographicLib ไว้ก่อนแต่เมื่อเปิดไปแล้วไม่เจอโฟลเดอร์ geographiclib\geoids ผมเลยสร้างใหม่) ไลบรารีเวลาประมวลผลจะดูที่โฟลเดอร์นี้อย่างเดียว

เปิด command prompt ของวินโดส์ทำการคำนวณโดยเรียกใช้ geoideval ดังตัวอย่างต่อไปนี้

geoideval -n tgm2017-1 --input-string "14.0 100.0"
-32.2897
geoideval -n tgm2017-1 --input-string "7 100.5"
-14.5592
geoideval -n tgm2017-1 --input-string "16.5 102.0"
-30.4109
geoideval -n tgm2017-1 --input-string "18.0 98.5"
-38.4128

สรุปสั้นๆถ้าใครต้องการคำนวณหาความสูงจีออยด์ TGM2017 ตอนนี้สามารถเลือกได้หลากหลาย ต้องการใช้ GeographicLip ก็ตามบทความนี้ ต้องการใช้แบบ GUI ก็ต้อง Surveyor Pocket Tools

ขอขอบคุณดร.ไพศาล สันติธรรมนนท์มาอีกครั้งที่ให้โอกาสทดสอบไฟล์ รวมทั้งสร้างโอกาสให้ศึกษารูปแบบการจัดเก็บไฟล์รูป PGM ในเชิงลึกรวมทั้งมีโอกาสได้เขียนโค้ดสั้นๆเพื่อย้อนรอยว่าไฟล์ PGM สร้างขึ้นมาได้อย่างไร โปรดติดตามกันตอนต่อไปครับ