การเปลี่ยนแปลงครั้งใหญ่ของไลบรารี Proj.4

วันนี้มาพูดถึงไลบรารี Proj.4 แบบลึกๆกันหน่อย บทความตอนนี้จะเป็นเรื่องโปรแกรมมิ่งนะครับ ไลบรารีตัวนี้ผมใช้เป็นแกนหลักในโปรแกรมรวมเครื่องมือฉบับกระเป๋าสำหรับช่างสำรวจ (Surveyor Pocket Tools) เอามาแปลงพิกัดกับระบบพิกัดที่ใช้กันในโลกนี้ (อาจจะได้ไม่ทั้งหมด) และไม่นานนี้ผมได้นำมาคำนวณ Vertical Datum คือสามารถหาความสูงจีออยด์ได้ ในความเป็นจริงถ้ามี Vertical Grid Shift หลายๆอันสามารถแปลงค่าระดับข้ามไปมาได้แบบที่ใช้ในอเมริกา

ก็ใช้มาหลายปีแล้ว ในตอนนี้ไลบรารีนี้ ตั้งแต่รุ่น 5.0.0 มาถึงรุ่นปัจจุบันรุ่น 6.0.0 ทางผู้พัฒนา Proj.4 ได้เริ่มวางรากฐานใหม่ ผมในฐานะคนนำมาใช้งานก็ติดตามลุ้นด้วยความระทึกละครับ เพราะถ้าแก้ใหม่แบบถอนรากถอนโคน ผมก็ต้องออกแรงมากหน่อยเพื่อมาปรับให้โปรแกรมสามารถใช้ไลบรารีรุ่นใหม่ได้เต็มประสิทธิภาพ

กว่าสองทศวรรษของการพัฒนา

มากกว่าสองทศวรรษที่ไลบรารี Proj.4 ได้ถูกพัฒนาและนำมาใช้งานในโปรแกรมเริ่มในปี 1980 แต่สามารถใช้ได้ในวงจำกัด สามารถใช้แปลงพิกัดได้ในระบบสองมิติ (2D) เท่านั้น ต่อมาได้เริ่มพัฒนาและสามารถแปลงพิกัดได้ในระบบสามมิติ (3D) ด้วยระยะเวลาการพัฒนาอันยาวนาน เปลี่ยนผ่านผู้พัฒนามาแล้วหลายรุ่นตามความต้องการในแต่ละกาลเวลา ไลบรารีสามารถใช้งานได้ในระดับที่ไม่ลึกซึ้งนักด้าน geodetic ในสมัยแรกๆ

การคำนวณก่อนหน้านี้การแปลงค่าพิกัดจะคำนวณผ่าน WGS84 เป็นศูนย์กลาง (hub) ซึ่งนิยามของ WGS84 นั้นค่อนข้างจะอ่อนแอ (ill-defined) ทุกๆครั้งที่อ้างอิงถึง WGS84 ก็ต้องมาถามกันก่อนว่าเป็น realization ปีไหน ซึ่งความไม่แน่ชัดตรงนี้ทำให้เป็นจุดอ่อนของพื้นหลักฐาน WGS84 มาอย่างยาวนาน

ปัจจุบันมีความต้องการการแปลงพิกัด 3D ที่ต้องการความละเอียดสูงเพื่อนำมาใช้งานด้าน รังวัด GNSS ตลอดจนนำมาประยุกต์ด้าน geospatial สำหรับงานด้านการเกษตร งานก่อสร้างและงานขนส่งต่างๆ ทำให้ผู้พัฒนาเล็งเห็นว่าจะต้องมีการรื้อโครงสร้างไลบรารีขนานใหญ่เพื่อให้บรรลุวัตถุประสงค์ดังกล่าวข้างต้นนี้

Transformation Pipeline

นี่เป็นคำเรียกเฟรมเวิร์คใหม่สำหรับการแปลงค่าพิกัด แปลงค่าระดับ แปลงหน่วยทั้งหลาย จะอยู่ในส่วนนี้ทั้งหมดที่ทางผู้พัฒนาได้ออกแบบใหม่ทั้งหมด หลักการคือการแปลงทั้งหลายส่งไปคำนวณเหมือนน้ำที่ไหลไปตามท่อ การแปลงแต่ละขั้นจะอิสระแยกจากกัน แต่ผลลัพธ์ของการแปลงหนึ่งครั้งอาจจะเป็นข้อมูลนำเข้า input ของการคำนวณในขั้นต่อๆไป ไหลไปเป็นสาย

ผมพยายามหาคำที่เป็นภาษาไทยมาแปลแต่สุดท้ายหาคำที่เหมาะสมไมไ่ด้ก็ขอเรียกทับศัพท์ เฟรมเวิร์คนี้ทางผู้พัฒนาได้ออกแบบให้การคำนวณได้ผลลัพธ์ที่มีความถูกต้องด้านตำแหน่งและค่าระดับที่ละเอียดมากขึ้น ในปัจจุบันทราบกันดีว่าเปลือกทวีป (Plate tectonic) ของเรามีการเลื่อนทุกๆปี แต่ละทวีปจะมีการเลื่อนในทิศทางต่างกัน ดังนั้นหมุดสำหรับงานสำรวจชั้นสูงก็จะมีการเปลี่ยนแปลงค่าพิกัดตลอดเวลา (Full time varying transformations) ดังนั้นเฟรมเวิร์คใหม่ชุดนี้จะสนับสนุนการคำนวณในลักษณะนี้เรียกว่า 4D spatiotemporal มิติที่ 4 ก็คือเรื่องเวลาที่เกี่ยวข้องนั่นเอง

ทดสอบการใช้งาน

ผมจะทดสอบการใช้งาน transformation pipeline โดยใช้ภาษาไพทอนผ่าน pyproj รุ่น 2.1.12 ที่เรียกใช้ Proj.4 รุ่น 6.0

ติดตั้งผ่าน pip ใช้คำสั่ง pip install pyproj หรือถ้ามี pyproj รุ่นเก่าอยู่สามารถ upgrade ได้ด้วยคำสั่ง pip install pyproj –upgrade

โจทย์ทดสอบ

มีหมุด GNSS มีค่าพิกัดละติจูด 16.8761696139 และลองจิจูด 100.3651157722 หน่วยเป็นดีกรี บนพื้นหลักฐาน WGS84 มีความสูงเทียบกับทรงรี (Ellipsoidal height) 43.451 เมตร ต้องการแปลงค่าพิกัดเป็นระบบพิกัดยูทีเอ็มบนพื้นหลักฐาน Indian 1975 และหาความสูงของจุดนี้เทียบกับระดับน้ำทะเลปานกลาง (รทก/MSL) โดยที่กำหนดให้ใช้แบบจำลองความสูงจีออยด์ความละเอียดสูง TGM2017 โปรแกรมทดสอบก็ไม่ได้ยาวเหยียดสั้นกระชับตามแบบฉบับของไพทอนที่เรียบง่ายทรงพลัง มาบวกกับไลบรารีเทพแบบ Proj4 สองแรงแข็งขัน

โค้ดทดสอบ

มาดูโค้ดโปรแกรมกัน มีกว่าสิบบรรทัดแค่นั้นเอง

from pyproj import Transformer

xin = 100.3651157722; yin = 16.8761696139; zhin = 43.451
pipeline_str = '''+proj=pipeline 
                  +step +proj=unitconvert +xy_in=deg +xy_out=rad 
                  +step +proj=longlat +ellps=WGS84 +step +proj=vgridshift +grids=d:/test/tgm2017.gtx +inv
                  +step +proj=push +v_3				   
				  +step +proj=utm +zone=47 +a=6377276.345 +b=6356075.41314024 +towgs84=204.4798,837.8940,294.7765,0,0,0,0 +units=m
				  +step +proj=pop +v_3'''
pipe_trans = Transformer.from_pipeline(pipeline_str)
x1975, y1975, zH2017 = pipe_trans.transform(xx=xin, yy=yin, zz=zhin)

print('Input => Latitude = {:.8f} Longitude = {:.8f} h = {:.4f} m (ellipsoidal height)'.format(yin, xin, zhin))
print('Output => Northing = {:.4f} Easting = {:.4f} H={:.4f} m(MSL)'.format(y1975, x1975, zH2017))

ผมสร้างโฟลเดอร์ชื่อ test เอาไฟล์แบบจำลองความสูงจีออยด์ความละเอียดสูงของประเทศไทยปี 2017 รวมท้ั้งไฟล์โค้ดโปรแกรม “test_proj4_vxi.py”

โฟลเดอร์จัดเก็บโปรแกรมและแบบจำลองความสูงจีออยด์

ถ้าติดตั้ง miniconda หรือ anaconda สามารถรันโค้ดได้ง่าย เครื่องผมติดตั้ง miniconda เมื่อเปิดมาแล้ว ผมใช้คำสั่ง cd เข้าไปในโฟลเดอร์ “test”

miniconda

รันโค้ดทดสอบด้วยไพทอน

ใช้คำสั่ง python test_proj4_vxi.py จะได้ผลลัพธ์ดังนี้

ผลลัพธ์
Input => Latitude = 16.87616961 Longitude = 100.36511577 h = 43.4510 m (ellipsoidal height) 
Output => Northing = 1866056.3407 Easting = 645746.8892 H=9.1400 m(MSL)

โค้ด Transformation pipeline

pipeline_str = '''+proj=pipeline 
                  +step +proj=unitconvert +xy_in=deg +xy_out=rad 
                  +step +proj=longlat +ellps=WGS84 +step +proj=vgridshift +grids=d:/test/tgm2017.gtx +inv
                  +step +proj=push +v_3				   
				  +step +proj=utm +zone=47 +a=6377276.345 +b=6356075.41314024 +towgs84=204.4798,837.8940,294.7765,0,0,0,0 +units=m
				  +step +proj=pop +v_3'''

สำคัญที่สุดคือบรรทัดนี้ เริ่มต้นจาก +proj=pipeline เป็นการกำหนดหัวว่าให้คำนวณแบบ pipeline ต่อมา +step คือการกำหนดให้คำนวณแต่ละขั้นตอน เริ่มจากให้คำนวณแปลงหน่วย +proj=unitconvert จากดีกรีเป็น เรเดียน เพื่อจะส่งต่อเป็น input ให้ +proj=longlat ในพื้นหลักฐาน WGS84

ขั้นต่อไปจะเป็นการคำนวณหาความสูงจีออยด์จากแบบจำลองความสูงจีออยด์ TGM2017 ด้วย +proj=vgridshift ระบุพาทที่อยู่ของแบบจำลองด้วย +grid=d:/test/tgm2017.gtx +inv สุดท้ายจะได้ความสูง orthometric หรือเทียบกับระดับน้ำทะเลปานกลาง

ขั้นต่อไป +proj=push +v_3 คือการเก็บค่าความสูงเอาไว้ก่อนใน stack ถ้าไม่สั่งขั้นต่อไปจะนำค่านี้แปลงกันเป็นความสูงบนวงรีของพื้นหลักฐาน indian 1975

ต่อไป +proj=utm +zone=47 +a=6377276.345 +b=6356075.41314024 +towgs84=204.4798,837.8940,294.7765,0,0,0,0 +units=m เป็นการแปลงจากพื้นหลักฐาน WGS84 มาเป็นระบบพิกัดฉากยูทีเอ็มบนพื้นหลักฐาน indian 1975 ตามพารามิเตอร์ที่กำหนดโดยกรมแผนที่ทหาร

ขั้นตอนสุดท้าย +proj=pop +v_3 คือการเอาค่าความสูงที่เก็บไว้ออกมาจาก stack

ผลลัพธ์การแปลงพิกัดและความสูง

ก็ได้ค่าตาม output ที่ปริ๊นท์ออกมาด้วยไพทอนคือเป็นค่าพิกัดยูทีเอ็มบนพื้นหลัก Indian 1975 N = 1866056.3407 E = 645746.8892 ความสูง = 9.140 เมตร (รทก.)

การ implement เพื่อใช้ใน Surveyor Pocket Tools

ก็คงสักพักใหญ่ละครับที่จะนำไปอิมพลีเมนต์เข้าใช้ในโปรแกรม Surveyor Pocket Tools เพราะผมเพิ่งศึกษาไลบรารี Proj.4 ในเชิงลึก เห็นว่ามีการเปลี่ยนแปลงไปมาก ประมาณว่าทุบอาคารตึกเก่าออก แล้วสร้างฐานรากใหม่สร้างอาคารใหม่ที่ไฉไลกว่าเดิมกันทีเดียว

ที่เขียนบทความนี้ก็เพื่อเป็นแนวทางให้ผู้ที่สนใจที่ใช้ไพทอนสามารถนำไลบรารี Proj.4 ที่ผมเห็นว่าน่าจะเป็นไลบรารีที่ดีที่สุดในด้านการแปลงพิกัดจากระบบพิกัดที่ต่างกันในโลกนี้ จะเห็นว่าการใช้งานไม่ได้ยากเย็นนัก ติดตามกันตอนต่อไปครับ

2 thoughts on “การเปลี่ยนแปลงครั้งใหญ่ของไลบรารี Proj.4”

  1. ขอบคุณสำหรับการแบ่งปันความรู้ครับ

Leave a Reply

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