Update : โปรแกรม Surveyor Pocket Tools คำนวณความสูงจีออยด์ จากไฟล์ค่าพิกัดภูมิศาสตร์

ในกรณีที่ต้องการค่าความสูงจีออยด์จากจุดที่มีจำนวนมากตัวอย่างเช่นเป็นสิบจุดขึ้นไป การมานั่งคำนวณทีละจุดคงไม่ใช่เรื่องที่สะดวกนัก ผมปรับปรุงโปรแกรมให้สามารถอ่านไฟล์ค่าพิกัดภูมิศาสตร์ (ละติจูดและลองจิจูด) ในรูปแบบ CSV ที่ใช้ตัวแบ่งด้วยเครื่องหมายคอมมา “,” ค่าพิกัดละติจูดและลองจิจูด ต้องเป็นรูปแบบทศนิยม (degree) เท่านั้น การจัดเรียงค่าพิกัดของให้ขึ้นต้นด้วยค่าลองจิจูดตามด้วยเครื่องหมายคอมม่าและค่าละติจูด

ไฟล์ทดสอบ

ไฟล์ที่จะมาทดสอบโปรแกรม ผมสร้างจากโค้ดภาษาไพทอน ให้สุ่มจำนวนจุดค่าพิกัดขึ้นมา 10000 จุด โดยให้ค่าพิกัดที่สุ่มอยู่ในกรอบสี่เหลี่ยมนี้คือ ละติจูด เริ่มที่ 3 องศา ไปสุดที่ 22.9833333333333 องศา ค่าลองจิจูดเริ่มที่ 95 องศา ไปสิ้นสุดที่ 107.9833333333 องศา

ดาวน์โหลดและติดตั้ง

ผมปรับรุ่นโปรแกรมใหม่ในขณะนี้เป็นเวอร์ชั่น 1.03 build 651 ก็สามารถไปดาวน์โหลดได้ที่ลิ๊งค์นี้ เมื่อดาวน์โหลดแล้ว ก็ติดตั้ง ในกรณีมีรุ่นเดิมอยู่ให้ถอนออกเสียก่อน เมื่อเปิดโปรแกรมมาจะเห็นหน้าตาดังนี้

Surveyor Pocket Tools V1.03 build 651

เมื่อเปิดโปรแกรมมาแล้ว ลองไปดูว่าไฟล์ทดสอบอยู่ตรงไหน ที่เมนูหลักให้เปิด “Example folder”

เปิดโฟลเดอร์ตัวอย่างไฟล์

จะเป็นการเปิด Windows Explorer คลิกเข้าไปใน “example data” จะเห็นไฟล์ทดสอบชื่อ “tgm2017_test_points.csv” ดังรูปด้านล่าง

ไฟล์ทดสอบ tgm2017_test_points.csv

ผมลองเปิดไฟล์นี้ดูด้วย Notepad ++ เลื่อนไปที่บรรทัดสุดท้าย จะเห็นว่ามีค่าพิกัดสุ่มทั้งหมด 10000 จุด

ไฟล์ทดสอบมีจำนวนจุด 10000 จุด

ทดสอบคำนวณความสูงจีออยด์จากไฟล์แบบลากแล้ววาง

ที่เมนูหลักเปิดทูลส์ “Geoid Height” และเปิด Windows Explorer เพื่อจะลากไฟล์ทดสอบไปวางบนทูลส์ และไม่ลืมในการทดสอบนี้ผมต้องการลอง “TGM2017” ดังนั้นที่ทูลส์อย่าลืมเลือก EGM Model เป็น “TGM2017

ลากไฟล์ทดสอบไปวางบนทูลส์ “Geoid Height”
เปิด EGM Model เป็น “TGM2017” แล้วลากไฟล์ทดสอบมาวาง

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

ไฟล์ผลลัพธ์ชื่อ ” tgm2017_test_points_out.csv

ทดสอบการคำนวณจากการเปิดไฟล์

อีกวิธีถ้าไม่ใช้การลากแล้ววางสามารถคลิกที่ไอคอนของทูลส์ได้เพื่อเปิดไฟล์ทดสอบ

จะเห็นไดอะล็อกดังรูปด้านล่าง

คลิกที่ปุ่ม “Browse…” เพื่อเปิดไฟล์ CSV เปิดไฟล์ทดสอบชื่อ “tgm2017_test_points.csv”

คลิกชื่อไฟล์แล้วคลิกที่ปุ่ม “Open” จะเห็นเนื้อในไฟล์

คลิกที่ปุ่ม “Calculate” เพื่อทำการคำนวณหาความสูงจีออยด์ สุดท้ายทูลส์ “Geoid Height” จะสร้างไฟล์ชื่อ “tgm2017_test_points_out.csv” เมื่อเปิดไฟล์มาดูจะเห็นว่ามีคอลัมน์ ความสูงจีออยด์เพิ่มมาอีกคอลัมน์

ทดสอบกับไฟล์ขนาดใหญ่

สำหรับไลบรารี Proj4 ที่ผมนำมาใช้งานอยู่นี้ การคำนวณหาความสูงจีออยด์หรือจาก Vertical Datum ตัวหนึ่งไปอีก Vertical Datum แบบในอเมริกาที่มี Vertical Datum อยู่หลายตัวหลายรุ่น ไลบรารีนี้จึงถูกพัฒนาเพื่อให้คำนวณในเรื่องนี้ โดยเฉพาะการคำนวณจากไฟล์ไลดาร์ (Point cloud) ที่มีขนาดใหญ่ หลายสิบล้านจุด จากความสูงทรงรีมาเป็น Orthometric height (หรือระดับน้ำทะเลปานกลางในประเทศไทยเรา) ผมลองทดสอบกับไฟล์ขนาดที่มีจำนวนจุด หนึ่งแสนจุด, หนึ่งล้านจุดและสิบล้านจุดตามลำดับ และผมจับเวลาไว้ (เครื่องผม CPU Xeon E-2176M) ไฟล์ทดสอบนี้ผมสร้างจากไฟล์สุ่มเช่นเดียวกัน

จำนวนจุดเวลาประมวลผล (วินาที)
100,0001
1,000,0006
10,000,00084

จะเห็นว่าเวลาประมวลผล สิบล้านจุด ใช้เวลา 1 นาที 24 วินาที ผมดูว่าใช้เวลามากไปนิด ทั้งที่ก่อนหน้านี้ สิบล้านจุด ใช้เวลาสิบนาทีกว่าๆ ผมต้องกลับมาออปติไมซ์โค้ดอีกครั้งจนได้เวลาหนึ่งนาทีกว่าๆ อย่างที่เห็น ถ้าจะเอาเร็วกว่านี้คงต้องใช้ภาษา C/C++ แล้วครับ

เนื่องจากโปรแกรมที่ปรับปรุงในตอนนี้ยังไม่ได้ดักจับเรื่อง error มากนักโดยเฉพาะรูปแบบไฟล์ค่าพิกัดละติจูดและลองจิจูด ดังนั้นก่อนจะเปิดไฟล์มาทำการคำนวณให้พิจารณาดูให้ดีว่ารูปแบบได้ตามที่โปรแกรมต้องการไหม ครับตอนนี้ก็พอหอมปากหอมคอ ติดตามกันตอนหน้าต่อไปครับ

โปรแกรมรวมเครื่องมือฉบับกระเป๋าสำหรับช่างสำรวจ Surveyor Pocket Tools : Update คำนวณพื้นที่และ Scale Factor ด้วย TGM2017

สำหรับการคำนวณ Scale Factor ไม่ว่าจะเป็นจุดเดี่ยว (Point scale factor) หรือแบบเส้นตรงเฉลี่ย (Line scale factor) หรือไม่ว่าจะคำนวณพื้นที่จริงที่ทอนจากพื้นที่ตามระบบพิกัดฉากกริด ก็ตามผมเปิดโอกาสให้ผู้ใช้สามารถเลือกใช้ TGM2017 (Precise Geoid Model of Thailand 2017) เพื่อให้ค่าออกมาเนียนๆมากกว่าเดิม ที่เดิมทีใช้เฉพาะ EGM2008 เท่านั้น

การประยุกต์ใช้แบบจำลองความสูงจีออยด์ต่างๆ นอกจากใช้ในงานรังวัด GNSS ยังนำมาประยุกต์ใช้สำหรับคำนวณหาสเกลแฟคเตอร์ แบบ Elevation Scale Factor (ESF) ที่ภาษาไทยอาจจะเรียกว่า ตัวคูนมาตราส่วนความสูง เมื่อนำมาคูนกับ Grid Scale Factor (GSF) จะได้ตัวประกอบที่เรียกว่า Combined Scale Factor (CSF) = Elevation Scale Factor (ESF) x Grid Scale Factor (GSF) ตามไดอะแกรมด้านล่าง


การทอน scale factor จากพื้นโลก <–> ระบบพิกัดฉากยูทีเอ็ม

ตัวประกอบ CSF สามารถนำมาใช้ในการแปลงพื้นที่จากระบบพิกัดฉากเช่นยูทีเอ็มเป็นพื้นที่จริงๆบนพื้นโลก นอกจากนั้นในงานก่อสร้างระยะทางยาวๆเช่นงานก่อสร้างถนน สามารถนำค่านี้มาใช้การทอนระยะทางบนพื้นโลกที่วัดได้จากกล้องโททัลสเตชั่น ไปเป็นระยะทางบนระบบพิกัดฉาก เพื่อให้สอดคล้องตามแบบ drawing ที่ออกแบบไว้บนระบบพิกัดฉากอยู่แล้ว ค่า CSF อาจจะหาได้จากเส้นตรงนำมาเฉลี่ยหัวท้าย โดยใช้เครื่องมือ Line Scale Factor ที่อยู่ในเครื่องมือของโปรแกรมนี้ หรือใช้แบบจุดเดี่ยว Point Scale Factor ถ้าพื้นที่ก่อสร้างไม่กว้างใหญ่มากนัก

ระยะทางที่ต่างกันระหว่างระยะทางบนพื้นโลกและระยะทางบนระบบพิกัดฉากยูทีเอ็ม

สมัยปัจจุบันกล้อง total station มีความทันสมัยสามารถเอาค่า CSF มาป้อนลงที่กล้อง โดยที่ช่างสำรวจไม่ต้องกังวลเรื่องระยะทางที่ต่างกันระหว่างระยะทางบนพื้นโลกและระยะทางบนระบบพิกัดฉาก ตัวกล้องจะทำการคำนวณให้เอง สามารถนำกล้องไปเก็บงาน topographic หรือวางผังก่อสร้างได้

โปรแกรมรุ่น V1.02 Build 643

ก็ไปหน้าดาวน์โหลดกันได้ที่ ลิ๊งค์นี้ มองหารุ่นล่าสุด ปัจจุบันคือรุ่น 1.02 build 643 ผมได้ปรับปรุงเพิ่มเติมหลายๆอย่าง

ปรับรูปแบบการป้อนมุมให้ง่ายยิ่งขึ้น

รุ่นเก่าๆจะสังเกตุว่าการป้อนมุมจะแข็งขืน เมื่อผิดไม่สามารถลบถอยหลังแบบกดคีย์ Backspace ได้ แต่รุ่นใหม่มีความยืดหยุ่นมากกว่าเดิม ในด้านโปรแกรมมิ่งการป้อนค่าให้ตรงกับรูปแบบที่ตั้งไว้เรียกว่า Regular Expression (RegEx) โดยการสร้างกฎไว้ก่อนแล้วคอยดักว่าผู้ใช้กดคีย์ตรงกับที่กำหนดไว้หรือไม่ ในส่วนนี้ผมเองค่อนข้างจะมึนงงกับ RegEx อยู่พอสมควร ดังนั้นการป้อนค่าโดยเฉพาะมุมละติจูด ลองจิจูด อาจจะไม่ดีนักในรุ่นแรกๆ แต่จะดีและยืดหยุ่นขึ้นในรุ่นนี้ ยังไงอาจจะมีบั๊กบ้าง

เพิ่มค่าปริยายแบบจำลองความสูงจีออยด์

เพิ่ม Default Geoid Model คือตั้งแบบจำลองความสูงเป็นค่าปริยาย โดยสามารถไปตั้งค่าได้ที่ “settings

คลิก settings เพื่อตั้งค่าปริยาย “Default Geoid Model”

ไปที่แท็บ “Geoid Model” เป็น “TGM2017” เมื่อไปคำนวณด้วยทูลส์ตัวอื่นๆเช่น Point Scale Factor, Line Scale Factor หรือ Area โปรแกรมจะเลือกใช้แบบจำลองนี้เป็นค่าปริยาย

เลือก “TGM2017” เป็นค่าปริยาย

คำนวณ Point Scale Factor ใช้ TGM2017

ผมจะทดสอบเปรียบเทียบการใช้แบบจำลอง TGM2017 กับแบบจำลองจีออยด์ EGM2008 เพื่อให้เห็นความแตกต่าง ประเดิมด้วยการใช้ TGM2017 ที่จุดทดสอบดังนี้

ใช้แบบจำลองความสูงจีออยด์ TGM2017

เมื่อใช้แบบจำลองความสูงจีออยด์ EGM2008 ได้ค่าดังนี้


ใช้แบบจำลองความสูงจีออยด์ EGM2008

เมื่อใช้ TGM2017 ได้ค่า Combined Scale Factor (CSF) = 1.00012353591 ใช้ EGM2008 ได้ค่า CSF = 1.00012364867 จะเห็นว่าต่างกันที่ทศนิยมที่ 7 ซึ่งถือว่าต่างกันไม่มาก (ทศนิยมหกตำแหน่งไม่ต่างกันถือว่ายังใช้แทนกันได้) อย่างที่ผมเกริ่นไว้ตั้งแต่แรก ยังไงก็ตามค่า combined scale factor จาก TGM2017 ก็ยังถือว่าเนียนกว่า เพราะสังเกตุดีๆว่าค่าความสูงจีออยด์จุดนี้ต่างกันประมาณ เกือบ 80 ซม.


คำนวณ Line Scale Factor ใช้ TGM2017

มาลองทดสอบการคำนวณหา Line Scale Factor แบบเฉลี่ย เมื่อใช้แบบจำลองความสูง TGM2017 ได้ค่ามาดังนี้


ใช้แบบจำลองความสูงจีออยด์ TGM2017

เมื่อใช้แบบจำลองความสูงจีออยด์ EGM2008 ได้ค่าดังนี้


ใช้แบบจำลองความสูงจีออยด์ EGM2008

จุดทดสอบนี้ค่าความสูงจีออยด์ของ TGM2017 และ EGM2008 ต่างกันประมาณเกือบ 80 ซม. ค่า combined scale factor ต่างกันที่ทศนิยมที่เจ็ด ก็ยังถือว่าใช้แทนกันได้ แต่แนะนำให้ใช้ TGM2017 เพราะจะได้ค่าที่ถูกต้องกว่า เนียนกว่า


คำนวณพื้นที่

การคำนวณพื้นที่บนระบบพิกัดฉาก (Grid-based Area) เช่นยูทีเอ็ม จะต้องมีการทอนพื้นที่มาเป็นพื้นที่บนผิวโลก (Ground-based Area) โดยอาศัย combined scale factor (CSF) หลักการคือเอาพีื้นที่บนระบบพิกัดกริดนำมาตั้งและหารด้วยค่า CSF จะได้พื้นที่จริงๆ ที่เราใช้ซื้อขาย มาดูตัวอย่างในการคำนวณพื้นที่ เปรียบเทียบระหว่างแบบจำลองจีออยด์ TGM2017 และ EGM2008 ว่าจะมีความแตกต่างกันมากน้อยเพียงใด

การคำนวณผมเปิดไฟล์ตัวอย่างจากโฟลเดอร์ “example data” วิธีเข้าไปดูง่ายๆว่ามีไฟล์อะไรบ้าง คลิกที่ “Example folder” ที่รายการหลัก

จะเห็น Windows Explorer

วิธีคำนวณหาพื้นที่เปิดไฟล์ตัวอย่าง เครื่องคอมพิวเตอร์ผมอยู่ที่ “C:/Users/priabroy/AppData/Roaming/Surveyor Pocket Tools/example data/boundary2-utm47n-indian1975.csv” ไฟล์นี้ค่าพิกัดบันทึกในระบบพิกัด Indian 1975 ความสูงเฉลี่ยของพื้นที่ 203.500 เมตร (H) จากระดับน้ำทะเลปานกลาง หลักการคำนวณของโปรแกรม จะคำนวณหา centroid area ของพื้นที่นี้ให้ก่อน แล้วนำค่านี้ในระบบพิกัดฉากยูทีเอ็ม Indian 1975 ไปคำนวณหาค่าพิกัดในระบบพิกิดภูมิศาสตร์บนหลักฐาน WGS84 จากนั้นจะใช้ค่าพิกัดภูมิศาสตร์นี้เพื่อไปดึงค่าความสุงจีออยด์ (Geoid separation – N) จากแบบจำลองความสูง TGM2017

แปลงค่าระดับน้ำทะเลปานกลางที่กำหนดให้ (203.5 เมตร) เป็นระดับความสูงเมื่อเทียบกับทรงรี WGS84 (h) จากสูตร h = H + N = 203.5 – 30.8038 = 172.6962 และนำค่า h นี้ไปคำนวณหาค่า Elevation scale factor ได้ต่อไป ส่วนค่า Grid scale factor คำนวนได้จากสูตรเส้นโครงแผนที่อยู่แล้ว มาดูผลลัพธ์กัน


ใช้แบบจำลองความสูงจีออยด์ TGM2017

เมื่อใช้แบบจำลองความสูงจีออยด์ EGM2008 ได้ค่าดังนี้


ใช้แบบจำลองความสูงจีออยด์ EGM2008

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


ข้อจำกัดของ TGM2017

เนื่องจาก TGM2017 เป็นแบบจำลองจีออยด์ท้องถิ่นคือใช้เฉพาะบางพื้นที่ ในที่นี้คือครอบคลุมประเทศไทย ดังนั้นโปรแกรม Surveyor Pocket Tools พูดให้เจาะจงคือไลบรารี Proj4 ที่ผมใช้เป็นไลบรารีในการคำนวณอยู่ จะเป็นตัวตรวจว่าค่าพิกัดภูมิศาสตร์ที่ส่งมาคำนวณนั้นค่าอยู่ในโซนที่อนุญาตไหม ถ้าอยู่นอกจะ error และเตือนออกมา ตัวอย่างผมพยายามป้อนค่า latitude = 24 องศา และ longitude = 108 องศา เพื่อทำการคำนวณ

แสดงความผิดพลาดเมื่อป้อนจุดคำนวณอยู่นอกเหนือโซนที่แบบจำลองสนับสนุน

ขอบเขตที่ TGM2017 กำหนดไว้คือ ค่าละติจูดเริ่มจาก 3 องศา ไปจนสิ้นสุดที่ 22.98333333 องศา ส่วนค่าลองจิจูด เริ่มจาก 95 องศา ไปสิ้นสุดที่ 107.983333333 องศา ทดสอบได้ดังนี้

ถ้าเลยไปนิด จะ error

สรุปการใช้งาน

การใช้งานด้านคำนวณสเกลแฟคเตอร์อาจจะเห็นไม่ชัดนักเมื่อเปลี่ยนจาก EGM2008 มา TGM2017 แต่ก็แนะนำให้ใช้เพราะถูกต้องกว่า แต่การรังวัด GNSS ต่างๆต้องใช้เพราะค่าต่างระหว่างสองแบบจำลองจีออยด์นี้ในประเทศไทยมีค่าต่างประมาณ 70-80 ซม.

ยังติดค้างเรื่องคำนวณหาความสูงจีออยด์หรือ Geoid Separation (N) จากไฟล์ในกรณีมีหลายจุด ติดตามกันได้ต่อไปครับ

Update : โปรแกรม Surveyor Pocket Tools คำนวณความสูงจีออยด์ TGM2017

Update : โปรแกรม Surveyor Pocket Tools คำนวณความสูงจีออยด์ TGM2017

มาตามสัญญาที่ผมบอกว่าจะอัพเดท Surveyor Pocket Tools โปรแกรมช่างสำรวจฉบับกระเป๋า ให้สามารถใช้งานคำนวณความสูงจีออยด์ TGM2017 (Thailand Precise Geoid Model 2017) ดั้งเดิมสามารถคำนวณบนโมเดล EGM96 และ EGM2008 เพียงเท่านั้น

เปลี่ยนวีธีการคำนวณโดยใช้ไลบรารี Proj4

ดั้งเดิมตอนคำนวณหาความสูงจีออยด์บน EGM96 และ EGM2008 ผมพัฒนาโค้ดโปรแกรมมาจากการดัดแปลงของโค้ดดั้งเดิมภาษาซีของผู้พัฒนาท่านอื่นมาเป็นภาษาไพทอน ขณะนี้ได้ยกเลิกโค้ดชุดนี้หันมาใช้ไลบรารี Proj4 คำนวณให้ทั้งหมดเพราะสะดวกมาก และไลบรารี Proj4 สามารถจัดการเมมโมรีได้ดีกว่าของผมมาก สังเกตก่อนหน้านี้ถ้าใครเข้าไปดูใน Task Manager ของวินโดส์อาจจะตกใจที่โปรแกรม Surveyor Pocket Tools ของผมเขมือบเมมโมรีมหาศาล เพราะผมใช้วิธีอ่านไฟล์ของโมเดล EGM96 และ EGM2008 มาเก็บไว้ในเมมโมรีตอนเปิดโปรแกรม ถ้าใครใช้เครื่องคอมเก่าๆ อาจจะรำคาญตอนเปิดโปรแกรมเพราะรอนานมาก ขณะนี้ปัญหาเหล่านั้นได้หมดไป

ดาวน์โหลดและติดตั้ง

ตอนนี้ผมได้ปรับมาเป็นรุ่น 1.02 build 632 ก็ไปดาวน์โหลดกันได้ที่ลิ๊งค์นี้ มีทั้งสองเวอร์ชั่นให้เลือกคือ Windows 32 บิตและ Windows 64 บิต

Surveyor Pocket Tools V1.02 build 632

ตอนคำนวณความสูงจีออยด์ ไลบรารี Proj4 จะเปิดไฟล์โมเดลมาอ่านคำนวณหาความสูงจีออยด์ ซึ่งทำได้รวดเร็วมาก และโปรแกรม Surveyor Pocket Tools ในปัจจุบันก็กินเมมโมรีลงน้อยมาก ขณะเปิดโปรแกรมมายังไม่ได้คำนวณอะไรกินเมมโมรีประมาณ 95 MB ตอนคำนวณความสูงจีออยด์กินไปประมาณ 250 MB ก็ไม่มากมาย

โมเดล EGM96, EGM2008 และ TGM2017 ถ้าสังเกตุตอนติดตั้งโปรแกรมจะเห็นชื่อไฟล์คือ egm96_15.gtx, egm08_25.gtx และ tgm2017.gtx ไฟล์หลังสุดนี้ผมดัดแปลงจากไฟล์แอสกี้ของกรมแผนที่ทหารเป็นรูปแบบมาตรฐานของ NOAA ที่ไลบรารีต้องการ

ไฟล์โมเดลของจีออยด์หลังจากติดตั้งโปรแกรม

คำนวณความสูงจีออยด์

คลิกที่โปรแกรมย่อย “Geoid Height” จะเห็นไดอะล็อกหน้าตาดังรูปด้านล่าง เลือก “EGM Model” เป็น TGM2017

Geoid Model “TGM2017”

ป้อนค่า latitude, longitude ไปดังรูป คลิกลูกศรเพื่อทำการคำนวณจะได้ความสูงจีออยด์

ทดสอบคำนวณความสูงจีออยด์จุดที่ 1

ลองปักหมุด Google Maps ดู

ปักหมุดบน Google Maps

มาทดสอบจุดที่ 2 กันต่อ ป้อน latitude, longitude เข้าไปดังนี้

ทดสอบคำนวณจุดที่ 2

สิ่งที่จะพัฒนาต่อไป

ตอนนี้คำนวณแบบแมนวลจุดต่อจุด อนาคตผมจะให้สามารถเปิดไฟล์มาคำนวณได้เพื่อความสะดวก ส่วนโปรแกรมย่อยอื่นที่ใช้คำนวณหาความสูงจีออยด์เช่น Point Scale Factor, Line Scale Factor ยังใช้จีออยด์ EGM2008 เป็นค่าปริยาย ผมจะออกแบบพัฒนาโปรแกรมให้ผู้ใช้สามารถเลือกโมเดลจีออยด์ได้ ติดตามกันต่อไปครับ

ทดสอบคำนวณหาความสูงจีออยด์ TGM2017 ด้วยไลบรารี Proj.4

ทดสอบคำนวณหาความสูงจีออยด์ TGM2017 ด้วยไลบรารี Proj.4

ไม่นานมานี้มีผมดาวน์โหลดไฟล์โปรแกรมและข้อมูลของ 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 แถว
TGM2017 layout

สำหรับไลบรารี 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 เท่า

ไฟล์ tgm2017.gtx ขนาดประมาณ 3.57 MB

ตัวอย่างโปรแกรมไพทอนสำหรับคำนวณหาความสูงจีออยด์

ต่อไปโปรแกรมไพทอนที่จะนำมาทดสอบคำนวณหาความสูงจีออยด์ ผมตั้งชื่อง่ายๆว่า “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 ได้แก่การคำนวณความสูงจีออยด์ การคำนวณหาสเกลแฟคเตอร์

สรุป

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

แก้ไขจอดับ: Surveyor Pocket Tools แสดงแผนที่บน Google Maps (สัจธรรมโลกนี้ไม่มีอะไรฟรี)

ตั้งแต่วันที่ 16 มิถุนายน 2018 ที่ผ่านมา ถ้าผู้ใช้ใช้งาน Surveyor Pocket Tools ลองปักหมุดผ่านโปรแกรมนี้ จะเห็นว่าจอดับแสดงข้อความว่า “For development purpose only” เพราะว่ากูเกิ้ลเปลี่ยนมาเก็บเงินผู้ใช้โดยเฉพาะผู้พัฒนาโปรแกรม โดยที่ผู้พัฒนาจะต้องขอ API Key จากทางกูเกิ้ลก่อน แล้วมาลงนามในบัญชีพร้อมจะเป็นหนี้กูเกิ้ล ถ้าใช้เกินกำหนด ผมมานั่งพิจารณาหลายเดือนว่าจะเอาไงดี โปรแกรมแจกให้ใช้ฟรีไม่มีเงื่อนไข แต่การที่ต้องไปจ่ายเงินให้กูเกิ้ลเพื่อการนี้ก็ดูกระไรๆ (ถึงแม้จะมีแบบ premium ให้องค์กรไม่หวังกำไรใช้ฟรี แต่ก็ยุ่งยากเกินไปสำหรับผม)

จอดับกรณีที่กูเกิ้ลเปลี่ยนนโยบาย Google Maps API

เอาละครับเมื่อคิดได้ปลงตก ก็มาดูว่ากูเกิ้ลเขาให้ใช้ฟรีได้เท่าไหร่เมื่อใช้แผนที่แบบ static สำหรับโปรแกรมบน PC Desktop ตอนแรกกูเกิ้ลให้ใช้ฟรีประมาณ 2500 ครั้งต่อวัน (นับตัวเลขที่มีการ request ไปที่ server ของกูเกิ้ล) ไม่นาน Wikipedia กับ Foursquare สองเจ้านี้ย้ายไปใช้แผนที่ตัวอื่น ทำให้กูเกิ้ลต้องลดราคาลง และเพิ่มจำนวนให้ใช้ฟรีเป็น 25000 ครั้งต่อวัน ผมดูแล้วโปรแกรม Surveyor Pocket Tools ที่ผู้ใช้วันๆหนึ่งแล้วปักหมุดคงไม่มากเท่านี้ ก็เลยมาลองดูสักตั้ง

เปลี่ยนวิธีการแสดงผล

เดิมผมใช้วิธีการสร้างไฟล์ html ที่เครื่องโลคอลคือสร้างไว้ที่เครื่องคอมพิวเตอร์ของผู้ใช้โปรแกรม เวลาผู้ใช้ปักหมุดแค่เรียกใช้ web browser เช่น chrome หรือ firefox ที่ติดตั้งอยู่ก่อนแล้ว เมื่อเปิดใช้ไฟล์ html ก็จะสามารถออนไลน์ดูแผนที่ที่ปักหมุดได้บน Google Maps แต่เมื่อมาจอดับ ผมไม่สามารถใช้วิธีการนี้ได้อีกต่อไป เพราะในไฟล์ html ผมต้องใส่คีย์ (API Key) ของผมลงไปด้วย เผื่อมีคนเปิดไฟล์นี้เอาคีย์ของผมไปใช้เกินอัตราก็บรรลัยจ่ายเงินกันอาน

ในการพัฒนาโปรแกรมนี้ผมใช้ PySide2 ดูว่ามีอะไรแสดงผลเว็บได้ ก็พอดีที่ QtWebEngineView ที่ใช้เอนจิน Chrome ตัวใหม่ เอาละก็ต้องมาแก้ไขโค้ดโปรแกรมกันใหม่ออกแรงกันนิด สุดท้ายผมไล่แก้ในโปรแกรมย่อยๆทุกตัว ตรงที่ปักหมุดจะแสดงผลแบบรูปด้านล่างทั้งหมด สามารถซูมเข้าออกได้ตามความต้องการ

แสดงผล Google Maps จากเอนจินเว็บของ PySide2
ชยายละเอียดในโหมด Map
ดูภาพถ่ายดาวเทียมในโหมด Satellite

ผมปรับโปรแกรมเป็นรุ่น 1.01 ผมอัพโหลดโปรแกรม Survey Pocket Tools รุ่น 1.01 ตัวใหม่ให้ดาวน์โหลดไปใช้ ไปที่หน้าดาวน์โหลดได้ตามลิ๊งค์นี้ (มีเฉพาะรุ่น 64 บิต ส่วน 32 บิตมีปัญหา build ไม่ผ่าน) ติดตามกันต่อไปครับ

Surveyor Pocket Tools version 1.01

Update โปรแกรมคำนวณวงรอบ Traverse Pro รุ่น 2.73 (ฉลองครบรอบ 20 ปี)

ผมกลับมาอัพเดทโปรแกรมคำนวณวงรอบ Traverse Pro อีกครั้งหลังจากรุ่นล่าสุดคือรุ่น 2.63 ที่ทิ้งไว้หลายปี ก็ถือโอกาสมาปรับปรุงเพื่อฉลองครบรอบวันเกิดโปรแกรมนี้ 20 ปี ซึ่งเริ่มพัฒนาเมื่อปี 1999 ด้วย Delphi ในขณะนั้น ก่อนที่จะย้ายมาพัฒนาด้วย free pascal + Lazarus ในภายหลัง ถ้านับอายุของโปรแกรมแล้วเป็นหนุ่มเต็มตัว แต่ผมผู้พัฒนาโปรแกรมเองสังขารเริ่มจะโรยรา 🙂 ก็ว่ากันไปครับ โลกนี้ไม่มีอะไรจีรังยังยืน

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

โปรแกรมแสดงรายการคำนวณวงรอบรุ่นเก่า
ส่งรายการคำนวณมาที่ Excel

ปรับปรุงการแสดงผลลัพธ์ของรายการคำนวณวงรอบ

มาดูรายการที่ปรับปรุง อันดับแรกได้แก่การแสดงผลรายการคำนวณวงรอบบนโปรแกรม เดิมทีผมใช้ของฟรีที่มีคนแจกไว้ ในตอนนี้ผมเลือกใช้ไลบรารีหรือคอมโพเนนต์ของ Lazarus จากเครื่องมือพัฒนาโปรแกรมที่ผมใช้อยู่เอง เนื่องจากมีความสามารถมากกว่าที่ผมเห็นเมื่อห้าปีที่แล้ว มาลองทดสอบโปรแกรมดู เริ่มต้นด้วยเปิดไฟล์จากตัวอย่างที่ติดมากับตัวโปรแกรมคือไฟล์ “example2.trv” แล้วทำการคำนวณจะได้ผลลัพธ์ดังนี้

แสดงผลลัพธ์รายการคำนวณวงรอบจากโปรแกรมที่ปรับปรุงล่าสุด (version 2.73)

จะเห็นว่าดูแล้วสบายตา ด้านบนจะแสดงข้อมูลวงรอบ ตลอดจนความผิดพลาดทางมุม (Error angle) สุดท้ายเป็นความถูกต้องของวงรอบ (closing ratio) ส่วนในตารางจะแสดงเป็นแถบม้าลายตามสมัยนิยมเพื่อให้อ่านง่าย

ปรับปรุงการจัดเก็บไฟล์ผลลัพธ์

รุ่นเดิมสามารถส่งข้อมุลออก Microsoft Excel ได้อย่างเดียว แต่รุ่นใหม่นี้ผมจะไม่ใช้วิธีการยิงหรือส่งข้อมูลไปออก excel เพียงอย่างเดียว แต่สามารถจัดเก็บเป็นไฟล์ได้หลายๆรูปแบบ เช่น CSV, LibreOffice/Openoffice หรือ Microsoft Excel ซึ่งอาจจะมีผู้ใช้บางท่านที่ไมได้ใช้ Microsoft Excel แต่ใช้ชุดออฟฟิศของ LibreOffice หรือ OpenOffice ก็ยังเอาไปเปิดได้ ถ้าดูที่แถบเครื่องมือจะเห็นไอคอนเพิ่มมาจากรุ่นก่อนตามรูปด้านล่าง เมื่อคำนวณแล้วลองคลิกที่ไอคอนนี้

แถบเครื่องมือ แสดงไอคอนสำหรับจัดเก็บไฟล์ผลลัพธ์งานคำนวณที่วงไว้ด้วยวงสีแดง

จัดเก็บไฟล์ในรูปแบบ Microsoft Excel

เลือกรูปแบบไฟล์ที่ต้องการจัดเก็บไฟล์

ผมลองเลือกจัดเก็บไฟล์เป็น Microsoft Excel ก่อน ชื่อไฟล์คือ “Example2.xlsx” จากนั้นเปิดไฟล์นี้ด่วย Excel จะได้ผลลัพธ์ดังนี้

เปิดไฟล์ผลลัพธ์รายการตำนวณวงรอบด้วยโปรแกรม Microsoft Excel

ในตอนนี้ถ้าต้องการพิมพ์ออกเครื่องพิมพ์ ตอนเขียนไฟล์นี้ด้วยโปรแกรมวงรอบ จะตั้งขอบเขตพื้นที่พิมพ์ (print area) รวมทั้งระยะกั้นหน้ามาให้เรียบร้อยแล้ว รวมถึง footer ที่แสดงจำนวนหน้ามาด้วย ลองใช้เมนู File > Print ดูครับ


Print Preview ใน Excel

เลื่อนไปดูที่ Footer ด้านล่าง

Footer แสดงหมายเลขหน้า

จัดเก็บไฟล์ในรูปแบบ LibreOffice & OpenOffice Calc

เนื่องจากทั้ง LibreOffice และ OpenOffice เป็นโปรแกรมฟรีทั้งคู่ผ่านการพัฒนามาอย่างยาวนาน ในปัจจุบันความสามารถ ความใช้งานง่ายก็ใกล้เคียงกับโปรแกรมเชิงพาณิชย์อย่างไมโครซอฟท์ออฟฟิศ ต่อไปมาลองจัดเก็บไฟล์ในรูปแบบนี้ดู

เลือกรูปแบบเป็น LibreOffice/OpenOffice spreadsheet

ลองเปิดด้วย LibreOffice Calc

เปิดไฟล์ผลลัพธ์คำนวรวงรอบด้วย LibreOffice Calc

จากนั้นเมื่อต้องการพิมพ์ ก็ดูตัวอย่างการพิมพ์ก่อน ที่เมนูเลือก File > Print Preview…

สิ่งที่จะปรับปรุงในอนาคต

สำหรับ Traverse Pro ออกแบบมาคำนวณวงรอบแบบ single loop ก็คงจะไม่มีฟีเจอร์เพิ่มอะไรมากมายในอนาคต แต่เล็กๆน้อยๆที่มองๆไว้ คือเขียนไฟล์ผลลัพธ์ลง google earth KML แต่ก็จนใจยังหาไลบรารีไม่ได้ อนาคตถ้ามีคนเขียนไลบรารีพวกนี้และแจกฟรี ผมก็จะนำมาประยุกต์ใช้ในโปรแกรม Traverse Pro นี้ด้วย หรือถ้าพอมีเวลาว่างจากงานประจำ ผมอาจจะ port โค้ดจากไพทอนมาก็ได้

ดาวน์โหลดโปรแกรมวงรอบ Traverser Pro

ไปที่หน้าดาวน์โหลด (Download) ได้ มองหา Traverse Pro รุ่น 2.73 ถ้าท่านผู้อ่านอยากจะให้ปรับปรุงเพิ่มเติมอะไรก็มาโพสต์ลงได้ที่ท้ายบทความนี้ได้ครับ ถ้าไม่เหลือบ่ากว่าแรงก็จะปรับปรุงให้ในลำดับต่อไป พบกันตอนหน้าครับ

ทดสอบเขียนโปรแกรมไพทอน (Python)  บนเครื่องคิดเลข Casio fx-cg50 Prizm

ทดสอบเขียนโปรแกรมไพทอน (Python) บนเครื่องคิดเลข Casio fx-cg50 Prizm

ไพทอนบนเครื่องคิดเลข

ช่วงนี้ผมมีโอกาสทำงานใกล้ชิดกับภาคสนาม ทำให้มีโอกาสได้จับและใช้เครื่องคิดเลขมากกว่าปกติ ในเวลาที่ผ่านมาไม่ถึงเดือนผมได้ซื้อเครื่องคิดเลข Casio fx-CG50 Prizm เคสสีขาว ที่ซื้อมาเพราะทราบว่าถ้า update OS เป็นรุ่น 3.20 จะสามารถใช้ ไพทอน (Python) ได้ ก็ขอหมายเหตุสักนิดว่าเป็นไมโครไพทอน (Micropython) ที่ทางทีมงาน Micropython ได้พอร์ตออกมาให้มีขนาดเล็กเพื่อเอาไปรันในบอร์ด iOT ได้ หรือบอร์ดที่เป็นไมโครคอนโทรลเลอร์ทั้งหลาย เน้นขนาดเล็ก หน่วยความจำต่ำ กินไฟน้อย ต่ออินเทอร์เน็ตได้ในตัว ผมจะไม่มุ่งไปทางนี้หรอกครับ ในบทความนี้ แต่จะพูดถึงเครื่องคิดเลขคาสิโอ ที่นำเอาไมโครไพทอนมาลงเครื่องคิดเลขรุ่นนี้ เพราะว่าไมโครไพทอนกินหน่วยความจำต่ำ ก็เลยเหมาะสมที่จะเอามารันในเครื่องคิดเลขที่มีทรัพยากรต่ำอยู่แล้ว ให้เกิดประสิทธิภาพมากยิ่งขึ้นไป

เป็นที่ทราบกันว่าไพทอนในปัจจุบันมีความนิยมมากยิ่งขึ้น เนื่องจากเรียนรู้ได้ง่าย เขียนง่าย อ่านง่าย ทรงพลัง แต่สำหรับไมโครไพทอน เนื่องจากออกแบบให้มีขนาดเล็ก ดังนั้นไลบรารีต่างๆที่ใช้ได้ในไพทอนรุ่นใหญ่ จะสามารถนำมาใช้กับไมโครไพทอนได้ต้องมีการพอร์ตใหม่ อาจจะลดฟังก์ชั่นหรือลดสิ่งไม่จำเป็นออก ทำให้มีขนาดเล็กลง ดังนั้นผมทำใจในจุดนี้อยู่แล้ว ว่าไม่สามารถเอาไลบรารีรุ่นใหญ่มารันได้แน่นอนเช่น pyproj, mathplotlib

สำหรับเครื่องคิดเลขที่ใช้ในแวดวงวิศวกรรม โปรแกรมที่เขียนด้วยไมโครไพทอนที่มีไลบรารี math หรืออาจจะเสริมด้วยไลบรารีพวกเมตริก (matrix) น่าจะพอนำมาเขียนใช้งานกันได้

 การพัฒนาโปรแกรมด้วยชุดพัฒนาโปรแกรมภาษาซี (Software Development Kit)

นี่เป็นความข้องใจของผมในฐานะแฟนเครื่องคิดเลขคาสิโอ รุ่น fx-9860G ทางคาสิโอจัดทำ SDK ให้สามารถเขียนโปรแกรมด้วยภาษาซี ที่ผมเขียนโปรแกรมมาแจกใช้ในแวดวงงานสำรวจแล้วหลายโปรแกรม แต่รุ่นนี้กลับไม่ทำมาให้  (ที่จริงไม่ทำมาให้ตั้งแต่ fx-CG10/fx-CG20)ไม่ทราบว่าเพราะเหตุใด เครื่องคิดเลขรุ่นนี้ไม่มี SD Card แต่ทดแทนด้วยการใส่ Flash memory มา 16 MB  ซึ่งก็พอจะใส่โปรแกรมใช้งานได้มากโขอยู่ หรือคิดว่ามีไมโครไพทอน มาให้แล้วน่าจะตอบโจทย์ได้หมด แต่ผมก็ไม่คิดอย่างนั้น ยังมีโปรแกรมเมอร์ภาษาซีอีกพอสมควร และในสภาพแวดล้อมของเครื่องคิดเลขจริงๆ โปรแกรมที่เขียนด้วยภาษาซีจะเร็วกว่าไพทอนอยู่แล้ว แต่ไพทอนได้เปรียบในด้านความง่าย

 เครื่องมือพัฒนาโปรแกรมของชุมชน

ยังมีชุมชนของนักพัฒนาที่สร้าง SDK  ขึ้นมาใช้งานเอง มีประมาณ 2-3 กลุ่มแต่สุดท้ายดูเหมือนไม่มีความเคลื่อนไหวกันมาหลายปีแล้ว เครื่องมือที่มีชื่อเสียงมากที่สุดคือ PrizmSDK และอีกอันคือ mini-SDK ผมเองใช้ไลบรารี MyLib แต่เผอิญผู้พัฒนาได้ทำไว้สำหรับเครื่อง fx-9860G เท่านั้น ไม่เป็นไรขอมุ่งลองไพทอนบนเครื่องคิดเลขรุ่น fx-CG50 นี้ก่อน ถ้าพัฒนาโปรแกรมด้วยภาษาซี ผมก็ยังมุ่งไปที่เครื่องคิดเลข fx-9860G เหมือนเดิม 

ผมลองเขียนโปรแกรมทดสอบเล็กๆลองดูด้วยเครื่องมือ PrizmSDK ก็ได้ดังรูปข้างล่าง (โปรแกรมไม่มีอะไรมีแต่เมนู) เทียบกับโปรแกรม System Manager ที่มากับเครื่อง

คุณสมบัติของเครื่องคิดเลข

โดยรวมรวมแล้วเครื่องคิดเลขนั้นเหมาะสำหรับนักศึกษามาก เพราะมีฟังก์ชั่นคณิตศาสตร์ การเงิน สถิติ มีกราฟมากมายให้ใช้ แต่สำหรับผมแล้วไม่มีอะไรต้องใช้เลย ยกเว้นเรื่องโปรแกรมบนเครื่องคิดเลขอย่างเดียว ถ้าไม่มีสิ่งนี้ก็ใช้เป็นที่ทับกระดาษได้เลย เครื่องรุ่นนี้ใช้โปรเซสเซอร์ตระกูล SH4 ขนาดหน้าจอ 384 x 216 จอ LCD จำนวนสี 65000 สี มีความสว่างพอสมควรและปรับได้ ความกว้างหน้าจอแบบทะแยง 3.17 นิ้ว หน่วยความจำของเครื่อง 60 KB มี Flash memory ที่สามารถเขียนอ่านได้ 16  MB ซึ่งจะเป็นที่เอาไว้เก็บโปรแกรมหรือข้อมูล ใช้ถ่าน AAA 4 ก้อน  เท่าที่ผมเปิดเครื่องใช้บ้างในแต่ละวันมาประมาณสองสัปดาห์ พบว่าแบตเตอรี่ลดลงมานิดหนึ่ง อนาคตอาจจะหาถ่านชาร์จมาใช้ ตอนนี้ใส่อัลคาไลน์ไปก่อน

ประเดิมโปรแกรมด้วยไพทอน

จะลองโปรแกรมทั้งทีผมพยายามให้โปรแกรมมีขนาดซับซ้อนมานิดหนึ่ง  และเรียกใช้โมดูลด้วย คิดไปคิดมาก็เลยจะลองโปรแกรมแปลงค่าพิกัดระหว่างค่าพิกัดภูมิศาสตร์กับค่าพิกัดยูทีเอ็ม เนื่องจากไม่ค่อยมีเวลาเขียน เลยลองหาไลบรารีที่ท่านอื่นได้ทำไว้ ผมเคยเกริ่นไปแล้วข้างต้นว่าไลบรารีรุ่นใหญ่เช่น pyproj ไม่สามารถเอามาใช้ได้ ลองค้นดูพบว่ามีไลบรารีไพทอนเล็กๆ ชื่อ utm มีสัญญาอนุญาตเป็น MIT-License ผมเอาโค้ดมาดัดแปลงนิดหน่อยให้เหมาะสมกับเครื่องคิดเลข แล้วเขียนไปอยู่ในไฟล์ utm.py เพื่อให้สะดวกเวลาเรียกใช้

การเขียนโปรแกรมไพทอน ถ้าไปเขียนบนเครื่องคิดเลข จะชักช้าเสียเวลาครับ เนื่องจากไมโครไพทอน พอร์ตไลบรารีเช่น math แล้ว ดังนั้นถ้าโปรแกรมของเราไม่ได้ใช้อะไรพิศดารมาก ก็สามารถมาเขียนโค้ดบน PyCharm หรือ Idle ได้ ผมเลือก PyCharm เมื่อทดสอบโปรแกรมเสร็จสามารถ โอนโปรแกรมเข้าไปไว้ในเครื่องด้วยการต่อเครื่องคิดเลขด้วยสาย USB เข้ากับคอมพิวเตอร์ที่รันวินโดส์ตามผังด้านล่าง

เมื่อรันโปรแกรมได้ตามความต้องการแล้วจากนั้นก็เอาเครื่องคิดเลขมาต่อกับคอมพิวเตอร์ วินโดส์จะมองเห็นเป็นไดรว์ สามารถใช้ File Explorer ก๊อปปี้โปรแกรมจากคอมพิวเตอร์ไปยังเครื่องคิดเลขได้ ตัวโค้ดโปรแกรมไลบรารีดูได้ด้านล่างครับ

import math
#from utm.error import OutOfRangeError

__all__ = ['to_latlon', 'from_latlon']

K0 = 0.9996

E = 0.00669438
E2 = E * E
E3 = E2 * E
E_P2 = E / (1.0 - E)

SQRT_E = math.sqrt(1 - E)
_E = (1 - SQRT_E) / (1 + SQRT_E)
_E2 = _E * _E
_E3 = _E2 * _E
_E4 = _E3 * _E
_E5 = _E4 * _E

M1 = (1 - E / 4 - 3 * E2 / 64 - 5 * E3 / 256)
M2 = (3 * E / 8 + 3 * E2 / 32 + 45 * E3 / 1024)
M3 = (15 * E2 / 256 + 45 * E3 / 1024)
M4 = (35 * E3 / 3072)

P2 = (3. / 2 * _E - 27. / 32 * _E3 + 269. / 512 * _E5)
P3 = (21. / 16 * _E2 - 55. / 32 * _E4)
P4 = (151. / 96 * _E3 - 417. / 128 * _E5)
P5 = (1097. / 512 * _E4)

R = 6378137

ZONE_LETTERS = "CDEFGHJKLMNPQRSTUVWXX"

class OutOfRangeError(ValueError):
    pass

def to_latlon(easting, northing, zone_number, hemi):

    northern = (hemi == 'N')

    x = easting - 500000
    y = northing

    if not northern:
        y -= 10000000

    m = y / K0
    mu = m / (R * M1)

    p_rad = (mu +
             P2 * math.sin(2 * mu) +
             P3 * math.sin(4 * mu) +
             P4 * math.sin(6 * mu) +
             P5 * math.sin(8 * mu))

    p_sin = math.sin(p_rad)
    p_sin2 = p_sin * p_sin

    p_cos = math.cos(p_rad)

    p_tan = p_sin / p_cos
    p_tan2 = p_tan * p_tan
    p_tan4 = p_tan2 * p_tan2

    ep_sin = 1 - E * p_sin2
    ep_sin_sqrt = math.sqrt(1 - E * p_sin2)

    n = R / ep_sin_sqrt
    r = (1 - E) / ep_sin

    c = _E * p_cos**2
    c2 = c * c

    d = x / (n * K0)
    d2 = d * d
    d3 = d2 * d
    d4 = d3 * d
    d5 = d4 * d
    d6 = d5 * d

    latitude = (p_rad - (p_tan / r) *
                (d2 / 2 -
                 d4 / 24 * (5 + 3 * p_tan2 + 10 * c - 4 * c2 - 9 * E_P2)) +
                 d6 / 720 * (61 + 90 * p_tan2 + 298 * c + 45 * p_tan4 - 252 * E_P2 - 3 * c2))

    longitude = (d -
                 d3 / 6 * (1 + 2 * p_tan2 + c) +
                 d5 / 120 * (5 - 2 * c + 28 * p_tan2 - 3 * c2 + 8 * E_P2 + 24 * p_tan4)) / p_cos

    return (180/math.pi*(latitude),
            180/math.pi*(longitude) + zone_number_to_central_longitude(zone_number))


def from_latlon(latitude, longitude, force_zone_number=None):
    if not -80.0 <= latitude <= 84.0:
        raise OutOfRangeError('latitude out of range (must be between 80 deg S and 84 deg N)')
    if not -180.0 <= longitude <= 180.0:
        raise OutOfRangeError('longitude out of range (must be between 180 deg W and 180 deg E)')

    lat_rad = math.pi/180*(latitude)
    lat_sin = math.sin(lat_rad)
    lat_cos = math.cos(lat_rad)

    lat_tan = lat_sin / lat_cos
    lat_tan2 = lat_tan * lat_tan
    lat_tan4 = lat_tan2 * lat_tan2

    if force_zone_number is None:
        zone_number = latlon_to_zone_number(latitude, longitude)
    else:
        zone_number = force_zone_number

    #zone_letter = latitude_to_zone_letter(latitude)
    if (latitude >= 0):
      hemi = 'N'
    else:
      hemi = 'S'

    lon_rad = math.pi/180*(longitude)
    central_lon = zone_number_to_central_longitude(zone_number)
    central_lon_rad = math.pi/180*(central_lon)

    n = R / math.sqrt(1 - E * lat_sin**2)
    c = E_P2 * lat_cos**2

    a = lat_cos * (lon_rad - central_lon_rad)
    a2 = a * a
    a3 = a2 * a
    a4 = a3 * a
    a5 = a4 * a
    a6 = a5 * a

    m = R * (M1 * lat_rad -
             M2 * math.sin(2 * lat_rad) +
             M3 * math.sin(4 * lat_rad) -
             M4 * math.sin(6 * lat_rad))

    easting = K0 * n * (a +
                        a3 / 6 * (1 - lat_tan2 + c) +
                        a5 / 120 * (5 - 18 * lat_tan2 + lat_tan4 + 72 * c - 58 * E_P2)) + 500000

    northing = K0 * (m + n * lat_tan * (a2 / 2 +
                                        a4 / 24 * (5 - lat_tan2 + 9 * c + 4 * c**2) +
                                        a6 / 720 * (61 - 58 * lat_tan2 + lat_tan4 + 600 * c - 330 * E_P2)))

    if latitude < 0:
        northing += 10000000

    return easting, northing, zone_number, hemi


def latitude_to_zone_letter(latitude):
    if -80 <= latitude <= 84:
        return ZONE_LETTERS[int(latitude + 80) >> 3]
    else:
        return None


def latlon_to_zone_number(latitude, longitude):
    if 56 <= latitude < 64 and 3 <= longitude < 12:
        return 32

    if 72 <= latitude <= 84 and longitude >= 0:
        if longitude <= 9:
            return 31
        elif longitude <= 21:
            return 33
        elif longitude <= 33:
            return 35
        elif longitude <= 42:
            return 37

    return int((longitude + 180) / 6) + 1


def zone_number_to_central_longitude(zone_number):
    return (zone_number - 1) * 6 - 180 + 3

โปรแกรมแปลงพิกัดภูมิศาสตร์ในภาคไพทอน

ผมเขียนไพทอนเป็นโมดูลอีกโมดูลเพื่อเรียกใช้ไลบรารี ตั้งชื่อว่า UTM2GEO.py โดยที่เขียนเมนูติดต่อการใช้งานง่ายๆ

from utm import *
   
def print_menu():
  print(5*'-',"MENU",5*'-')
  print('1: UTM to GEO')
  print('2: GEO to UTM')
  print('0: Exit')
  
def geo2utm(lon,lat):
  east,north,zn,hem=from_latlon(lat,lon)
  print("North={0:11.3f}".format(north))
  print("East={0:10.3f}".format(east))  
  print("UTM Zone No={0:0d}{1}".format(zn,hem))
    
def utm2geo(e,n,zoneno,hemi):
  lat,lon=to_latlon(e,n,zoneno,hemi)
  print("Latitude={0:11.7f}".format(lat))
  print("Longitude={0:10.7f}".format(lon))   
      
loop=True
while loop:      
  print_menu()
  choice=int(input('Selection[0-2]'))
  if (choice==0):
    loop=False
  elif (choice==1):
    loop=True
    y=float(input("Northing="))
    x=float(input("Easting="))
    zn=int(input("Zone No="))
    hem=input("Hemi (N/S)=")
    utm2geo(x,y,zn,hem)
  elif (choice==2):
    loop=True
    y=float(input("Latitude="))
    x=float(input("Longitude="))
    geo2utm(x,y)   
 

วิธีก๊อปปี้โปรแกรม

จากนั้นผมก๊อปปี้สองไฟล์คือ utm.py  และ UTM2GEO.py ลงบนไดรว์เครื่องคิดเลขดังนี้

ผมเก็บไว้ที่ไดเรคทอรี \SAVE-F\PROGRAM เวลาจะถอดสาย USB เพื่อเลิกการเชื่อมต่อต้อง Safely removal โดยการคลิกเมาส์ขวา “Eject” ที่ File Explorer จากนั้นมาที่เครื่องคิดเลขจากให้กดคีย์ “EXE” และ “EXIT” ตามลำดับ ถ้าไม่ทำไฟล์อาจจะไม่ได้ซิงค์กันอาจจะหายหรือไม่สมบูรณ์ได้ ที่เครื่องคิดเลขกดคีย์ “MENU” เลือก “Python

แปลงพิกัดจากค่าพิกัดภูมิศาสตร์เป็นค่าพิกัดยูทีเอ็ม

จากรูปด้านบนกดคีย์ F1-Run จะเห็นหน้าจอขึ้นเมนูติดต่อมาง่ายๆ

เราจะเลือกแปลงพิกัดจากค่าพิกัดภูมิศาสตร์ไปเป็นค่าพิกัดยูทีเอ็มเลือกกดคีย์ “2” ที่เครื่องคิดเลขแล้วกดคีย์ “EXE” ป้อนค่าพิกัด Latitude = 39.95259668 Longitude  = -75.15132081 (ป้อนเป็นหน่วยดีกรี ในตอนนี้ยังไม่รับค่าแบบอื่น) จะได้ผลการคำนวณออกมา เนื่องจากในตอนนี้ไม่มีคำสั่งเบรคการแสดงผลเมื่อเขียนด้วยไพทอน (เอาละจะมาบ่นทีหลัง ว่าใส่ไพทอนมาแล้วทางคาสิโอไม่ให้เครื่องมืออะไรมาเลย) การจะดูผลลัพธ์ ผู้ใช้ต้องกดคีย์ “0” เพื่อออกจากโปรแกรมและใช้ลูกศรกดขึ้นไปทางด้านบนเพื่อไปดูผลลัพธ์

จะได้ค่า Northing = 4422506.896 Easting = 487074.371 อยู่ในโซน 18N

แปลงพิกัดจากค่าพิกัดยูทีเอ็มเป็นค่าพิกัดภูมิศาสตร์

ทำการรันโปรแกรมใหม่อีกครั้ง ที่เมนูเลือกกดเลข “1” ป้อนค่าพิกัด Northing = 2642783.110 Easting =232030.949 UTM Zone No = 46 Hemi = N

ดูค่าผลลัพธ์ได้ (กดคีย์ “0” ออกจากเมนูก่อนแล้วเลื่อนขึ้นไปดู)

สรุปการใช้งาน

ตอนแรกผมคาดหวังจากที่ทางคาสิโอเอาไมโครไพทอนมาลงเครื่องคิดเลขรุ่นนี้ ยังไงการเขียนโปรแกรมใช้งานยังไงๆผลลัพธ์ที่ออกมาก็ต้องระดับเทพ เพราะไพทอนมันทรงพลังด้วยตัวของมัน แต่เมื่อลองแล้วผิดหวังมาก จนบัดนี้คาสิโอ้ยังไม่ได้ออกคู่มือแสดงฟังก์ชั่นที่ไพทอนสามารถเรียกมาใช้ได้ มีฟังก์ชั่น input กับ print สองฟังก์ชั่นนี้เท่านั้น เพียงแค่ผมค้นหาฟังก์ชั่น clear screen หน้าจอยังทำไม่ได้ ฟังก์ชั่นที่ต้องการสนับสนุนได้แก่การเขียนเมนูที่เรียกใช้ด้วยคีย์ F1 ถึง F6 การปริ๊นท์แสดงสีต่างๆ การเรียกใช้ฟังก์ชั่นกราฟต่างๆหรือพล็อทกราฟ หรือใช้งานเมตริก เป็นต้น

เอาละตอนนี้ไพทอนที่ปรากฎบน OS รุ่น 3.20 เพิ่งออกมาเตือนตุลาคม 2018 (ขณะที่เขียนบทความนี้เดือนพฤศจิกายน 2018) คงต้องให้เวลาสักพักว่าจะเป็นอย่างไร บอกตามตรงว่าคงต้องเอาเครื่องคิดเลขรุ่นนี้มาทับกระดาษอีกสักพักใหญ่ๆ

การเล็งสกัดย้อน (Resection) ด้วยการวัดมุมภายใน ระยะทางและและมุมแบริ่งด้วยวิธีการคำนวณแบบ Least Squares (ตอนที่ 2)

 ตั้งสมการ Observation Equation

ขอทบทวน ค่า aik, bik  เรียกว่า  direction coefficients และ  cik, dik เรียกว่า distance coefficients

ในกรณีวัดมุมเล็งสกัดย้อนจากสมการด้านบนและเอาแทนที่ในสมการด้านล่าง

เขียนให้ดูง่ายดังนี้  zi  คือค่าอะซิมัทเริ่มต้น

เราจะมาคำนวณหาค่า aik, bik กันก่อน มาคำนวณที่จุด P ค่าพิกัด N = 193939.897 E = 110879.464 ไปสถานีหลักฐานจุดที่ 1  ที่มีค่าพิกัด N = 192315.290 E = 120383.500 คำนวณระยะทางได้ 9641.890 เมตร คำนวณหาอะซิมัทได้ 99°42’1.1″ ดังนั้น aik = aP1 = -sin( 99°42’1.1″) / (9641.890 * 100) * 3600 * 180/3.141592654 = -0.2109 second (มุมแปลงเป็นหน่วยฟิลิปดา ระยะทางแปลงหน่วยเป็น ซม.)

bik = bP1 = cos(99°42’1.1″) / (9641.890 * 100) * 3600 * 180/3.141592654 = -0.0360 second

เราจะฟอร์มสมการในรูปแบบ v + Bx = f โดยที่

v คือเมตริกเวคเตอร์ของ residual

B คือเมตริกของค่า coefficient ของมุม ระยะทาง

f คือเมตริกความต่างของค่าที่คำนวณและค่าที่รังวัด

ผมคำนวณหาค่าaik, bik ทุกๆการรังวัดมุมมาดังนี้

มาทบทวนดูสมการระยะทางดังนี้้

คำนวณหา cik = cP1 = cos(99°42’1.1″) = -0.1685

dik= dP1 = sin(99°42’1.1″) = 0.9857

คำนวณทุกการรังวัดระยะทางมาได้ดังตาราง

สมการสุดท้ายคือการวัดเล็งสกัด (Bearing intersection)

จากสมการ v + Bx = f มาพิจารณาฝั่งซ้าย v + Bx ก่อน ผมฟอร์มเป็นเมตริกดังนี้

มาดูเมตริก เป็นเมตริกแสดงความต่างระหว่างค่าที่คำนวณกับค่ารังวัด f = Computed – Observation มาดูตามตาราง น่าจะเข้าใจได้ง่าย ค่า diff ในตารางก็คือค่า f นั่นเอง โปรดระวังหน่วยมุมจะเป็นฟิลิปดา (second) หน่วยระยะทางใช้เป็นซม.

ยกตัวอย่างค่าของเมตริก f เริ่มจากการวัดเล็งสกัดย้อน จะใช้มุมอะซิมัทที่ได้จากการคำนวณมาเป็นตัวเริ่มต้น ค่าจากจุด P ไปสถานีหลักฐานที่ 1 จะได้ค่าความต่างเท่ากับ 0.0 ต่อไปจากจุด P ไปสถานีหลักฐานที่ 2 คำนวณได้ 119.5116959 จากมุมที่วัดมา 99° 42′ 1.1″ + 19° 48′ 41″ = 119.5116957 ได้ค่าความต่าง = (119.5116959 – 119.5116957) * 3600 =  0.0008 second

ที่ง่ายที่สุดคือวัดระยะทางจากจุด P ไปสถานีหลักฐานจุดที่ 1 ได้ 9641.795 เมตร ส่วนการคำนวณจากค่าพิกัด P เริ่มต้นมายังค่าพิกัดสถานีหลักฐาน 1 ได้ค่าคำนวณ = 9641.890 เมตร ความต่าง = (9641.890 – 9641.795) * 100 = 9.5 ซม.

สุดท้ายสามารถนำค่ามาเขียนเป็นเมตริก ดังนี้

v + Bx = f

จะเห็นว่าสมการนี้ติดค่า residual ของเมตริก v และเมตริก x ไม่สามารถคำนวณต่อไป แต่หัวใจของ least squares ตามชื่อเลยครับคือผลรวมค่ายกกำลังสองของ  residual ที่ได้ค่าน้อยที่สุด

ถ้าค่า weight หรือน้ำหนักของการรังวัดไม่เท่ากันจะต้องคูณน้ำหนักเข้าไปด้วย

จากสมการ v + Bx = f แทนค่า v = f – Bx ในสมการ

ค่า Φ จะมีค่าน้อยที่สุด ดังนั้นจะหาอนุพันธ์ (ดิฟ) โดยที่ให้ x มีค่าเท่ากับศูนย์

ที่นี้ก็จำง่าย Nx = t โดยที่ N = BTWB และ f = BTW

สุดท้ายสามารถหาค่าเมตริก x = N-1เมื่อได้ค่า x แล้วก็สามารถย้อนไปหาเมตริก residual (v) ได้ ตามสมการ v + Bx = f

ผลลัพธ์การคำนวณรอบที่ 1

ผมใช้ฟังก์ชั่นของ excel หาเมตริกได้ดังนี้


ได้เมตริก x คือค่าปรับแก้หน่วยเป็นซม. เอาพิกัดจุด P และมุมอะซิมัท เริ่มต้นมาปรับได้ดังนี้ N = 193939.897 – 0.090829= 193939.806 ค่า E = 110879.464 + 0.04695 =  110879.511 และค่ามุมอะซิมัท = 99°42’1.1″ – 2.1156″ = 99° 41′ 58.99″ (99.6997191)

ผลลัพธ์การคำนวณรอบที่ 2

เอาค่าพิกัดของจุด P ที่ได้จากรอบที่ 1 มาเป็นตัวเริ่มต้น พร้อมทั้งมุมอะซิมัทด้วย

ตั้งสมการเมตริก v + Bx = f ได้ดังนี้

แก้สมการ Nx = t ได้ดังนี้

จะได้เมตริก x ค่าใหม่ เอาพิกัดจุด P และค่าอะซิมัทมาปรับได้ดังนี้ N =193939.806  – 0.0016 = 193939.790  และ E = 110879.511 + 0.0085 = 110879.519  ค่าอะซิมัท = 99° 41′ 58.99″ + 0.137″ = 99°41′ 59.13″ (99.6997572)

ผลลัพธ์การคำนวณรอบที่ 3

ต่อไปฟอร์มรูปเมตริก v + Bx = f สังเกตว่าเมตริก B ค่าเปลี่ยนไปเล็กน้อยมาก

คำนวณหาเมตริก x จากสมการ Nx = t => x = N-1t

จะได้เมตริก x ค่าใหม่ เอาพิกัดจุด P และค่าอะซิมัทมาปรับได้ดังนี้ N = 193939.790- 0.0003 = 193939.787  และ E = 110879.519 + 0.0001 = 110879.521 ค่าอะซิมัท = 99° 41′ 58.99″ + 0.024″ = 99°41′ 59.15″ (99.6997639)

มาถึงตอนนี้จะเห็นว่าค่าปรับแก้ ในเมตริก x  น้อยมากอยู่ในระดับเศษส่วนของมิลลิเมตร ΔN = -0.285 cm ΔE = 0.152 cm ดังนั้นแสดงว่าค่าที่คำนวณมานั้น convergence แล้ว ดังนั้นผมสรุปว่าผลลัพธ์ดังนี้

ค่าพิกัดจุด P (Free Station)

N = 193939.787 E = 110879.521 ค่าอะซิมัทจากจุด P ไปสถานีหลักฐานที 1 = 99°41′ 59.15″

 ตรวจสอบผลลัพธ์การคำนวณด้วย Microsurvey StarNet

เพื่อให้มั่นใจว่าผลลัพธ์ที่ได้จะถูกต้อง ผมใช้ Microsurvey StarNet  มาเป็นตัวช่วย ข้อมูลสถานีรังวัดไม่เกิน ดังนั้นผมยังใช้เวอร์ชั่นทดลองใช้ได้อยู่ เมื่อเปิดโปรแกรมมาผมใช้เมนู Options -> Project ตั้งค่าดังนี้ เปลี่ยนหน่วยเป็นเมตรให้เรียบร้อยก่อนที่ Adjustment > Units > Linear > Meters

จากนั้นปรับ  standard error  ของกล้องวัดมุมและระยะทางให้สอดคล้อง ผมปรับระยะทาง Distance constant 4mm ปรับ Distance PPM = 20 ความหมายตัวนี้คือ 4mm ± 20mm/1,000,000 x L (m) ถ้าวัดระยะทาง 1000 เมตร error จะอยู่ประมาณ 24 mm (กล้อง Total Station สมัยปัจจุบันเรื่องระยะทางทำได้ดีกว่านี้มาก)

จากนั้นป้อนข้อมูลไปดังนี้ หมายเหตุว่า C = Control Point ไม่มีเครื่องหมาย ! !  ตามหลังแสดงว่าเป็นค่าเริ่มต้นหรือประมาณการ A = Angle, D = Distance และ B = Bearing ถ้าสนใจก็ไปดาวน์โหลดโปรแกรมของ Microsurvey StarNet  มาทดลองได้

แล้วคัดลอกข้อมูลด้านล่างแล้วไปวางลงในหน้า  input เพื่อลองคำนวณดูผลลัพธ์กันได้

# Resection 2D combined resection, distance and bearing intersection.

# Approximate coordinates for the unknown free station
#
C P 193939.897 110879.464

# Coordinates for the known stations
C 1 192315.290 120383.500 ! !
C 2 189545.730 118642.430 ! !
C 3 188084.770 112278.210 ! !
C 4 190640.940 109654.540 ! !
C 5 190044.860 108065.980 ! !
C 6 194455.370 110632.930 ! !
A P-1-2 19-48-41
A P-1-4 100-40-19
A P-1-5 116-8-36
A P-1-6 234-44-22

#Distance measurements
D P-1 9641.795
D P-3 6019.802
D P-5 4804.793

#Bearing measurements
B 1-P 279-41-59.5
B 3-P 346-33-52
B 5-P 35-50-34

จากนั้นใช้เมนู Run > Adjustment

สรุปสถิติการคำนวณดังนี้ จะเห็นว่ามีจำนวนข้อมูลรังวัดมา 10 (วัดมุมเล็งสกัดย้อน 4 มุม ระยะทาง 3 ระยะ วัดมุมเล็งสกัด 3 มุม จำนวนสิ่งที่ไม่รู้ค่า 2 ค่า คือค่าพิกัด  x,y ของจุด P

Adjustment Statistical Summary
==============================

Iterations = 2

Number of Stations = 7

Number of Observations = 10
Number of Unknowns = 2
Number of Redundant Obs = 8

 

Adjusted Coordinates (Meters)

Station N E Description
P 193939.788830 110879.521213
1 192315.290000 120383.500000
2 189545.730000 118642.430000
3 188084.770000 112278.210000
4 190640.940000 109654.540000
5 190044.860000 108065.980000
6 194455.370000 110632.930000

จะได้ค่าพิกัดจุด  P ดังนี้ N = 193939.789 E = 110879.521 ต่างจากที่ผมคำนวณมาเล็กน้อย อย่างแรกคือผมไม่ได้ใช้น้ำหนัก  weight  แต่ของ MicroSurvey Starnet บังคับใช้ โปรแกรมไม่ได้แสดงค่าสมการให้ดู จึงตรวจสอบไม่ได้ว่าเมตริกของ W เป็นอย่างไร อย่างที่สองคือการคำนวณวนลูปสองรอบเท่านั้น แสดงว่าอัลกอริทึ่มอาจจะใช้อนุกรมเทเลอร์ลำดับที่ 2 ด้วย ซึ่งผมใช้ลำดับเดียวจึงต้องวนลูปมากกว่า

เป็นอย่างไรบ้างครับ เรื่องคำนวณ  least squares adjustment ตอนสมัยเรียนมหาวิทยาลัยถือว่าเป็นวิชาที่ยากสำหรับนักศีกษา แต่ถ้ามีโอกาสได้ลองคำนวณดูไปทีละขั้นตอน สำคัญคือต้องลงมือเองด้วย จึงจะเข้าใจ บางทีรายละเอียดเช่นอนุกรมเทเลอร์ การหาค่าอนุพันธ์ อาจจะยาก ผมไม่ได้จำหรอกครับ เพียงแต่ติดก็กลับมาเปิดตำรา หลงๆลืมๆไปบ้างตามวัยของผม แต่ปัจจุบันความได้เปรียบคือตำราหาอ่านได้ง่ายในอินเทอร์เน็ต โปรแกรมช่วยคำนวณแบบ  Microsurvey Starnet หรือโปรแกรมค่ายอื่นๆก็มีมากมาย ที่จะมาช่วยทำให้การคำนวณง่าย แต่ถ้ามีพื้นฐานบ้างก็จะได้เปรียบ  ติดตามกันต่อไปครับ

การเล็งสกัดย้อน (Resection) ด้วยการวัดมุมภายใน ระยะทางและและมุมแบริ่งด้วยวิธีการคำนวณแบบ Least Squares (ตอนที่ 1)

จากที่เขียนโปรแกรมสำหรับเครื่องคิดเลขคำนวณเล็งสกัดย้อนสำหรับเครื่องคิดเลข Casio fx-9860G II SD ทำให้นึกถึงวิธีการคำนวณแบบ least squares ที่เป็นพื้นฐานเคยร่ำเรียนมา โดยเฉพาะการรังวัดในปัจจุบันที่การรังวัดระยะทางด้วยกล้องประมวลผลรวมทำได้ง่าย เมื่อรวมกับการรังวัดมุม จะทำให้มีค่าเกินหรือ redundant มาคำนวณในวิธีแบบ least square ได้ การรังวัดแบบเล็งสกัดย้อนบางตำราเรียกว่า free station

การวัดเล็งสกัดย้อนในยุคแรกจะวัดมุมจากหมุดสถานีหลักฐานกันเป็นหลัก และหมุดสถานีหลักฐานต้องมีค่าพิกัดและมีจำนวนอย่างน้อยสามหมุด ดังนั้นการวัดมุมอย่างน้อยสองมุมจะเพียงพอต่อการมาคำนวณ ในบทความนี้ผมจะพาไปทัวร์การคำนวณวิธี least squares ด้วยการใช้การรังวัดผสมประกอบด้วยการวัดมุมภายใน วัดระยะทาง และการวัดมุมแบริ่งหรือการวัดภาคทิศ จากนั้นจะยกตัวอย่างพร้อมทั้งวิธีการคำนวณ ปิดท้ายทดสอบผลลัพธ์การคำนวณด้วยโปรแกรม Microsurvey StarNet

การวัดแบบ 2D

การวัดแบบนี้จะวัดมุมราบและระยะราบก็พอนำมาคำนวณหาค่าพิกัดของ free station ได้ ตัวอย่างที่ผมจะแสดงการคำนวณในลำดับต่อไปจะเป็นการวัดแบบ 2D เพื่อลดความซับซ้อน

การวัดแบบ 3D

การวัดระยะทางจะวัดแบบ Slope distance วัดมุมดิ่ง ถ้าทราบค่าระดับของหมุดสถานีหลักฐานและวัดความสูงของเป้าที่ตั้งบนสถานีหมุดหลักฐาน จากนั้นวัดความสูงกล้อง เมื่อนำมาคำนวณร่วมกับมุมราบแล้วจะได้ค่าพิกัดของ free station รวมทั้งค่าระดับด้วย

แนะนำตำรา

เมื่อผมกลับไปดูวิธีการคำนวณด้วย least square แต่กลับพบกับความนะจังงัง ว่าความรู้ที่รับการประสิทธิประสาทมาบัดนี้ได้คืนท่านอาจารย์ที่มหาวิทยาลัยไปจนหมดแล้วสิ้นเชิง นั่งงงกับการตั้งสมการ Observation Equation อยู่นานพอสมควร เมื่อไม่เป็นผลก็ต้องกลับไปหาตำราค่อยๆพลิกความทรงจำกลับมาใหม่

ตำราที่จะแนะนำให้อ่านเป็นของ Mr.Rod Deakin (Rodney Edwin Deakin) ถ้าดูตามโปรไฟล์ก่อนจะเกษียณเป็นอาจารย์ให้มหาวิทยาลัยสอนเรื่อง Engineering Survey มาก่อน ตำราหรือบทความที่แต่งนั้นมีความหลากหลายมากในเว็บไซต์ส่วนตัวนี้ อ่านง่าย เพลินจนลืมไปว่าเรื่องที่อ่านนั้นยาก

ตามไปดูหน้า least squares ลองเลื่อนไปด้านล่างดูเรื่อง “Notes on Least Squares” คือตำราที่จะมาแนะนำกันมีทั้งหมด 10 บท เนื่องจากผมพอมีพื้นฐานมาบ้างเล็กน้อย จึงไม่ได้อ่านเรียงหน้าตั้งแต่บทที่ 1 แต่อาศัยข้ามไปอ่านบทที่ 7 ก่อน พอเจอเมตริกผมงงก็ข้ามไปอ่าน “Appendix A” เฉพาะเรื่องเมตริก บวก ลบ คูณ อินเวอร์ส แล้วค่อยกลับมาอ่านบทอื่นๆอีกที ดังนั้นคนที่ห่างเรื่องนี้นานๆแนะนำให้ไปดูเรื่องเมตริกอันดับแรก

บทที่ 7 อนุกรมเทเลอร์

สำหรับการคำนวณเล็งสกัดย้อนด้วยวิธี least square ปฐมบทจะอยู่ที่บทที่ 7 เริ่มต้นจากอนุกรมเทเลอร์ เหตุที่เราต้องใช้เพราะว่าสมการที่จะนำมาคิด  resection นั้นไม่ใช่สมการแบบเชิงเส้นหรือ linear equation ดังนั้นจะต้องมีการถอดสมการเชิงเส้นออกมาและอยู่ในรูปอนุกรม เพื่อให้สามารถนำมาคำนวณได้ จะมีตัวอย่างสมการ Observation Equation 2 ตัวอย่างคือ

  • วัดมุม (Measure direction)
  • วัดระยะทาง (Measure distance)

บทที่ 9 คำนวณเล็งสกัดย้อน (Least Squares Resection)

บทนี้จะเริ่มตั้งแต่ตั้งสมการ  Observation Equation (ไม่เชิงเส้น) ของการวัดมุม จาก free station ไปยังหมุดสถานีหลักฐาน จากนั้นจะถอดสมการไม่เชิงเส้นด้วยอนุกรมเทเลอร์ แล้วก็มีตัวอย่างแสดงวิธีการคำนวณ

บทที่ 10 คำนวณเล็งสกัด (Least Squares Bearing Intersection)

ผู้อ่านอาจจะสังเกตสองคำคือเล็งสกัดกับเล็งสกัดย้อน

  • เล็งสกัดย้อน (Resection) คือไปตั้งกล้องที่  free station แล้ววัดมุมภายในไปหาสถานีหมุดหลักฐาน

  • เล็งสกัด (Bearing Intersection) จะเป็นการวัดมุมแบริ่ง ตัวอย่างได้แก่การวัดดาวเหนือในสมัยก่อน คือเราจะมีหมุดคู่เป็น base line ต้องการทราบมุมแบริ่งของ base line นี้ก็ตั้งกล้องที่หมุดตัวแรกแล้วอาศัยวัดดาวเหนือกับหมุดอีกตัวบน base line ดังนั้นถ้าวัดแบริ่งสัก 2 base line ไปตัดกันก็จะได้ค่าพิกัด

การถอดสมการเล็งสกัดย้อนด้วยอนุกรมเทเลอร์

ต่อไปมาลองดูสมการที่เราใช้คำนวณหามุมและระยะทาง โดยที่ทราบค่าพิกัด ซึ่งเป็นสมการพื้นฐานของงานสำรวจ ในที่นี้เป็นสมการไม่เชิงเส้น เราจะมาถอดสมการออกมาเป็นเชิงเส้นด้วยอนุกรมเทเลอร์ แผนผังด้านล่างเป็นงานเล็งสกัดย้อน ตั้งกล้องที่จุด Pi วัดมุมไปหาสถานีหลักฐานที่ทราบค่าพิกัดคือจุด  P1, P2, P3, … Pk

การวัดมุมในทางปฏิบัติจะส่องไปหาสถานีแรก P1 ตั้งมุมราบเป็น 0 แล้วกวาดไปสถานีที่ 2 วัดมุมได αi1 กวาดไปหาสถานีหลักฐาน P2 วัดมุมได้ α12 ทำเช่นนี้เรื่อยๆ จนถึงสถานีหลักฐาน Pk จะได้มุม αik แต่ในงานสำรวจด้านตรีโกณมิติเราจะไม่คำนวณโดยใช้มุมภายในแต่จะใช้มุมอะซิมัท ดังนั้นมุมภายในเหล่านี้ต้องเปลี่ยนเป็นเทอมของมุมอะซิมัททั้งหมด ซึ่งไม่ได้ยุ่งยากอะไร ถ้าทราบมุมอะซิมัทที่ไปหาสถานีแรก P1 ก็จะทราบอะซิมัทสถานีที่เหลือทั้งหมด แต่ในความเป็นจริงเราไม่ทราบ เพราะว่าค่าพิกัด P1 ยังไม่ทราบในตอนนี้

สำหรับสมการคำนวณหาอะซิมัทและระยะทางระหว่างจุด Pi และ Pk เมื่อทราบค่าพิกัดคำนวณได้ดังสมการ

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

โดยที่ E = E0 + ΔE และ N = N0 + ΔN ส่วน E0 และ N0 คือค่าเริ่มต้นเป็นค่าประมาณการ ส่วนค่า ΔE, ΔN เป็นค่าปรับแก้ (correction) ซึ่งสองตัวนี้จะคำนวณได้จากสมการ least squares ที่กำลังจะว่ากันต่อไป และถ้าดูสมการจะเห็นว่าติดค่าอนุพันธ์ (ดิฟ) ซึ่งค่าดิฟหาได้ดังนี้

จากสมการด้านบน aik, bik เรามีค่าประมาณการเริ่มต้น Ei, Ni ดังนั้นสามารถหาค่าพวกนี้ได้ มาดูสมการระยะทางค่าดิฟ หาได้ดังนี้

จากสมการด้านบน cik, dik เรามีค่าประมาณการเริ่มต้น Ei, Ni ดังนั้นสามารถหาค่าพวกนี้ได้เช่นเดียวกัน

ตั้งสมการ Observation Equation สำหรับการวัดมุมภายใน (เล็งสกัดย้อน)

สมการ Observation Equation เป็นสมการที่เป็นหัวใจของการคำนวณ least squares ถ้าตั้งสมการนี้ได้ก็สามารถคำนวณต่อได้ แต่งานรังวัดถ้าวัดเกินมาจะมี residual เป็นส่วนเกิน สมการ Observation Equation ของการวัดมุมอะซิมัทแสดงได้ดังนี้

vik คือ residual คงไม่ลืมกันนะว่า มุมอะซิมัท Øik จากสถานีตั้งกล้องไปยังแต่ละสถานีหมุดหลักฐานเท่ากับ อะซิมัทไปหมุดสถานีหลักฐานตัวแรก Zi บวกด้วยมุมราบที่กวาดไป αik เขียนในรูปอนุกรมเทเลอร์ดังนี้

ในกรณีงานรังวัดเล็งสกัดย้อน ΔNk และ ΔEk ทราบค่าพิกัดดังนั้นเทอมตัวนี้จะเป็นศูนย์ กลายเป็นสมการนี้

ตั้งสมการ Observation Equation สำหรับการวัดระยะทาง

คล้ายๆกับสมการสำหรับการวัดทิศทาง แสดงได้ดังนี้

เขียนในรูปอนุกรมเทเลอร์

ตั้งสมการ Observation Equation สำหรับการวัดแบริ่ง (เล็งสกัด)

ถือว่าเป็นอาหารเสริมก็แล้วกัน สมัยนี้ยุค GNSS  คงหายากแล้วสำหรับการวัดภาคทิศอะซิมัทด้วยการวัดดาว

สมการเดียวกับที่เราใช้เล็งสกัดย้อน แต่ในที่นี้เราทราบค่าพิกัด Nk, Ek  ทำให้ ΔNk และ ΔEk เทอมตัวนี้จะเป็นศูนย์ และสมการจะคงเหลือดังนี้

ตัวอย่างงานรังวัดเล็งสกัดย้อนแบบรวมวัดมุม ระยะทาง และรังวัดภาคทิศ

มาถึงตอนนี้ผู้อ่านบางท่านเห็นสมการแล้วอาจจะงงงวยเหมือนผมในตอนแรกๆ มาลองแก้สมการดูจากตัวอย่างโจทย์ ที่นี้จะทำให้เข้าใจมากขึ้น ดูรูปด้านล่าง ตั้งกล้องที่จุด P ทำการรังวัดโดยที่วัดมุมไปยังสถานีหลักฐาน โจทย์ตัวนี้เป็นโจทย์ที่ผมเจอในอินเทอร์เน็ต เมื่อค้นหาการคำนวณเล็งสกัดย้อนแบบวิธีคำนวณ least square ดูได้ที่ลิ๊งค์นี้ โจทย์ที่จะคำนวณเปิดไปดูได้หน้าท้ายๆ

กำหนดค่าพิกัดหมุดสถานีหลักฐานดังตารางด้านล่าง วัดมุมโดยที่เล็งไปที่สถานี 1 ตั้งมุมราบเป็นศูนย์แล้ววัดมุมไปที่สถานี 2, 4, 5 และ 6 ทำการวัดระยะทางไปที่สถานี 1, 3  และ 5  รังวัดภาคทิศ (วัดดาว) ที่สถานี 1, 3 และ 6 ดังตารางด้านล่าง

คำนวณหาค่าพิกัดจุดตั้งกล้องโดยประมาณ

ผมใช้เครื่องคิดเลข fx 9860 GII SD คำนวณหาค่าพิกัดของ P  ด้วยโปรแกรม Resection ป้อนค่าพิกัดและมุมได้ดังนี้

จุด A, B และ C  ของเครื่องคิดเลขก็คือจุดที่ 1, 2 และ 6 ของโจทย์ด้านบนตามลำดับ

ผมจะใช้ค่าจุด P ที่ค่าพิกัดโดยประมาณ N = 193939.897 E = 110879.464 เพื่อเป็นค่าเริ่มต้นในการคำนวณด้วย least squares

เนื่องจากบทความยาวมาก ติดตามกันตอนที่ 2 ต่อนะครับ

Update: โปรแกรมคำนวณเล็งสกัดย้อนฉบับปรับปรุง (Resection) สำหรับเครื่องคิดเลขเทพ Casio fx 9860G II SD

Update: โปรแกรมคำนวณเล็งสกัดย้อนฉบับปรับปรุง (Resection) สำหรับเครื่องคิดเลขเทพ Casio fx 9860G II SD

ผมเขียนเรื่องการคำนวณเล็งสกัดย้อน (Resection) จากตอนก่อนหน้านี้ด้วยอัลกอริทีมใหม่ของ Josep M. Font-Llagunes อ่านได้ที่

ติดปีกเครื่องคิดเลขเทพ Casio fx 9860G II SD ด้วยโปรแกรมภาษาซีบน AddIn ตอนที่ 5 โปรแกรมคำนวณ Resection ด้วยอัลกอริทึ่มสมัยใหม่


ในตอนนี้ผมจะมาเรียบเรียงโปรแกรมด้วยไลบรารี MyLib ที่ใช้ก่อนหน้านี้มาหลายโปรแกรมแล้ว โดยในครั้งนี้จะใช้อัลกอริทึมหลายๆตัวมาคำนวณเพื่อเปรียบเทียบกันดู เมื่อคำนวณแล้วมีรูปกราฟฟิคหยาบๆให้ดูพอกล้อมแกล้ม

การคำนวณเล็งสกัดย้อน หรือ Resection Problem มีมานานแล้วหลักๆแล้ววิธีการคำนวณแบ่งเป็น 2 วิธีหลักๆคือ แบบการตรีโกณมิติ (Trigonometric Solution) และวิธีการแบบเรขาคณิตวงกลมตัดกัน (Geometric Circle Intersection) วิธีการนี้นอกจากนำมาใช้ในงานสำรวจแล้ว ในแวดวงหุ่นยนต์ที่ใช้ในสถานที่ในร่ม ผู้อ่านอาจจะคุ้นกับการแข่งขันหุ่นยนต์ของน้องๆนักศึกษา การหาตำแหน่งในร่มของหุ่นยนต์จะอาศัยวิธีการนี้ จะอาศัยตั้งสถานีฐานที่รู้ค่าพิกัด 3 สถานี (Beacon) แล้วหุ่นยนต์จะมีอุปกรณ์จำพวก Goniometer ที่คอยหมุนกวาดหามุม 3 มุม จากตำแหน่งที่หุ่นยนต์ตั้งอยู่กับสถานีฐานนั้น จากนั้นจะคำนวณหาพิกัดตัวเองแบบ realtime

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

รายชื่อโปรแกรมเล็งสกัดย้อน

1.Kaestner-Burkhardt Method ใช้ตรีโกณมิติมาช่วยในด้านการคำนวณ วิธีการคำนวณนี้คิดค้นกันมานมนานแล้วและวิธีการคำนวณจะตั้งชื่อตามผู้ที่คิดค้น วิธีการนี้ได้มีคนนำมาปรับปรุงกันหลายครั้ง

2.TienStra Method ใช้วิธีการด้านตรีโกณมิติมาช่วยในการคำนวณ

 

3.Collins Method ใช้วิธีการด้านตรีโกณมิติมาช่วยในการคำนวณ

4.Cassini Method ใช้วิธีการด้านตรีโกณมิติมาช่วยในการคำนวณ

5.Font-Llagunes ใช้วิธีเรขาคณิตโดยใช้วงกลมมาตัดกันมาช่วยในการคำนวณ

ผมประยุกต์โดยการนำสูตรมาแปลงเป็นโค้ดภาษาซี  โดยได้นำวิธีการคำนวณทั้งหมด 5 วิธี แต่คงไม่สามารถอธิบายวิธีการคำนวณได้หมด

ดาวน์โหลดโปรแกรม (Download)

ไปที่หน้าดาวน์โหลด (Download) หรือไปที่โปรแกรม Resectionจะได้ไฟล์ชื่อ “RESCTNEX.G1A” จากนั้นทำการก็อปปี้ไฟล์เข้าเครื่องคิดเลขผ่านโปรแกรม Casio FA-124 หรือ copy ผ่าน SD Card ก็สะดวกเหมือนกัน

เริ่มใช้งานโปรแกรม

ที่เครื่องคิดเลขกดปุ่ม “MENU” จะเห็นลิสต์ของไอคอนของโปรแกรมต่างๆ ตัวโปรแกรม Resection อยู่ด้านล่างๆใช้ปุ่มลูกศร เลื่อนลงมาแล้วกดคีย์ “EXE

เมนูหลัก  (Main Menu)

ก่อนจะไปต่อวิธีการใช้งานมาดูเมนูหลักของโปรแกรม เนื่องจากเป็นคำย่อต้องทำความเข้าใจกันเสียก่อน

Met – Method (F1) เลือกวิธีการคำนวณที่ผมกล่าวมาข้างต้นจำนวน 5 วิธี

Coor – Coordinates (F2) ป้อนค่าพิกัดของเป้าหลัง (control point) จุดที่ทราบค่าพิกัด 3 จุด (A, B และ C ตามลำดับ)

Ang – Angle (F3) ป้อนค่ามุมสองมุม (α และ β ตามลำดับ)

Calc – Calculate (F4) คำนวณเล็งสกัดย้อนหาค่าพิกัดจุดตั้งกล้อง

Info – Information (F5) แสดงเครดิตไลบรารีที่ใช้งาน

Exit – Exit (F6) ออกจากโปรแกรมนี้

เลือกวิธีการคำนวณ (Method)

ที่เมนูหลักของโปรแกรมกดคีย์ F1 – Met เลือกวิธีการคำนวณ จะเห็นหน้าจอแสดง Resection Method

กดคีย์ “EXE” จะเห็น drop down list ลงมาให้เลือกดังรูปด้านล่าง ในที่นี้ลองเลือก “TienStra

กด F6 – OK เพื่อออก

ตัวอย่างที่ 1 กำหนดค่าพิกัดจุด A (1000E, 2000N) ค่าพิกัดจุด  B (1078.331E, 2077.869N) และค่าพิกัดจุด C (1172.191E, 2154.753N) วัดมุม α = 40°8’24” และมุม β = 57°36’0″ หาค่าพิกัดจุด P

ป้อนค่าพิกัดหมุดหลักฐาน (Input Control Points Coordinates)

ที่เมนูหลักกดคีย์ F2-Coor เพื่อป้อนค่าพิกัดหมุดหลักฐานที่ทราบค่าพิกัด 3 จุด ดังนี้

 

กดคีย์ F6 – OK เพื่อออก

ป้อนค่ามุม (Input Angles)

ที่เมนูหลักกดคีย์ F3 – Ang เพื่อป้อนมุม 2 มุมดังนี้ สำหรับการป้อนมุมให้ใช้เครื่องหมายลบ (-) คั่น

กดคีย์ F6 – OK เพื่อออก

คำนวณหาค่าพิกัด (Calculate Resection)

ที่เมนูหลักกดคีย์ F4 – Calc เพื่อคำนวณเล็งสกัดย้อน จะได้ผลลัพธ์ดังรูปด้านล่าง กดคีย์ F1 – PgUp หรือ F2 – PgDn เพื่อเลื่อนดูผลลัพธ์ โปรแกรมจะแสดงค่าพิกัดหมุดหลักฐาน 3 จุด พร้อมทั้งมุมที่ป้อนไปเพื่อทบทวนว่าผู้ใช้ป้อนเข้าไปถูกต้องตามที่ต้องการหรือไม่ และสุดท้ายก็แสดงค่าพิกัดจุดเล็งสกัดย้อน (Resection) ที่เราต้องการคือ (1167.446E, 2016.916N)

แสดงแผนผัง (Plot)

ต้องการดูแผนผังการจัดวางตัวของจุดหลักฐาน 3 จุดและผลลัพธ์ก็สามารถดูได้คร่าวๆ กดคีย์ F5 – Plot 

สามารถกดคีย์ F1 – Dn (Down) เพื่อเลื่อนภาพลง) F2 – Up (เลื่อนภาพขึ้น) F3 – Lt (Left เลื่อนภาพไปด้านซ้าย) F4 – Rt (เลื่อนภาพไปด้านขวา) หรือกด F5 – M-> (เลือกเมนูด้านขวา) เมนูด้านขว่าจะมี F1 -Z+ (เพื่อขยายภาพ) F2 – Z- (เพื่อย่อภาพ) เสร็จแล้วกดคีย์ F6 -Done เพื่อออก กดคีย์ F6 – Done อีกครั้งเพื่อออกมาเมนูหลัก

คำนวณด้วยวิธีการอื่น

การคำนวณด้วยวิธี Tienstra ได้นำเสนอไปแล้วลองมาใช้วิธีการอื่นมาลองคำนวณบ้าง ถ้าผู้ใช้ออกจากโปรแกรมด้วยการกดคีย์ F6 – Exit โปรแกรมจะทำการเก็บค่าต่างๆที่ป้อนไปทังหมดไว้ เมื่อเข้ามาใช้โปรแกรมอีกรอบ ถ้าต้องการใช้ค่าเดิมจะไม่ต้องป้อนใหม่อีกครั้ง

Kaestner-Burkhardt Method

ที่เมนูหลักกดคีย์ F1 -Met เพื่อเลือกวิธีการคำนวณ ในที่นี้เลือก Kaestner กดคีย์ F6 – OK เพื่อออก

ที่เมนูหลักกดคีย์ F3 – Calc เพื่อทำการคำนวณ (ไม่ต้องป้อนค่าพิกัดและค่ามุมใหม่) จะได้ผลลัพธ์มาดังนี้ กดคีย์ F1 – PgUp หรือ F2 – PgDn เพื่อเลื่อนดูผลลัพธ์

สังเกตว่าค่าพิกัดที่ได้ (1168.228E, 2014.113N) มีความต่างจากวิธี Tienstra (1167.446E, 2016.916N)ไปพอสมควร ผมไม่แน่ใจเมื่อตรวจสอบโค้ดที่แปลงจากสูตรก็ยังคิดว่าไม่มีที่ผิด ลองไปเทียบดูจากวิธีอื่นอีกที

 

Cassini Method

ที่เมนูหลักกดคีย์ F1 – Met เพื่อเลือกวิธีการคำนวณ เลือก “Cassini” กดคีย์ F6 – OK  เพื่อออก

ที่เมนูหลักกดคีย์ F3 – Calc เพื่อทำการคำนวณ (ไม่ต้องป้อนค่าพิกัดและค่ามุมใหม่) จะได้ผลลัพธ์มาดังนี้ กดคีย์ F1 – PgUp หรือ F2 – PgDn เพื่อเลื่อนดูผลลัพธ์

 

จะได้ค่าพิกัดผลลัพธ์ออกมาจากวิธีการของ Cassini (1167.446E, 2016.916N) ซึ่งเท่ากับวิธีการของ Tienstra

Collins Method

ที่เมนูหลักกดคีย์ F1 – Met เพื่อเลือกวิธีการคำนวณ เลือก “Collins” กดคีย์ F6 – OK  เพื่อออก

ที่เมนูหลักกดคีย์ F3 – Calc เพื่อทำการคำนวณ จะได้ผลลัพธ์มาดังนี้

จะได้ค่าพิกัดผลลัพธ์ออกมาจากวิธีการของ Collins (1167.446E, 2016.916N) ซึ่งเท่ากับวิธีการของ Tienstra, Cassini

 

 

Font-Llagunes

เป็นวิธีการที่ผมชอบมากที่สุดและเคยเขียนลงแบบละเอียดในตอนก่อนหน้านี้ การใช้งานทำได้ตามขั้นตอนดังต่อไปนี้ ที่เมนูหลักกดคีย์ F1 – Met เพื่อเลือกวิธีการคำนวณ เลือก “Font-Llagunes” กดคีย์ F6 –  OK เพื่อออก

กดคีย์ F3 – Calc เพื่อทำการคำนวณ จะได้ผลลัพธ์มาดังนี้ กดคีย์ F1 – PgUp หรือ F2 – PgDn เพื่อเลื่อนดูผลลัพธ์

เช่นเดียวกันจะได้ค่าพิกัดผลลัพธ์ออกมาจากวิธีการของ Font-Llagunes (1167.446E, 2016.916N) ซึ่งเท่ากับวิธีการของ Tienstra, Cassini และ  Collins แต่บางครั้งจะต่างกับวีธีการของ Kaestner-Burkhardt ในระดับหลักเมตร

 

 

ตัวอย่างที่ 2 กำหนดค่าพิกัดจุด A (1000E, 5300N) ค่าพิกัดจุด  B (2200E, 6300N) และค่าพิกัดจุด C (3100E, 5000N) วัดมุม α = 109°18’16.2″ และมุม β = 115°3’7.2″ หาค่าพิกัดจุด P

การคำนวณผมจะใช้การรวบรัดขอเสนอแต่วิธีการของ Font-Llagunes เท่านั้น ที่เมนูหลักกดคีย์ F1 – Met เพื่อเลือกวิธีการคำนวณ เลือก Font-Llangunes จาก drop down list กดคีย์ F6 – OK เพื่อออก

ที่เมนูหลักกดคีย์ F2- Coor เพื่อป้อนค่าพิกัด เสร็จแล้วกดคีย์ F6 – OK

ที่เมนูหลักกดคีย์ F3 – Ang เพื่อป้อนมุม เสร็จแล้วกดคีย์ F6 – OK

ที่เมนูหลักกดคีย์ F4 – Calc เพื่อทำคำนวณ กดคีย์ F1 – PgUp หรือ F2 – PgDn เพื่อเลื่อนดูผลลัพธ์

จะได้ค่าพิกัดจุดเล็งสกัดย้อนมาเท่ากับ (2128.989E, 5575.375N) กดคีย์ F5 – Plot เพื่อดูแผนผัง

ก็ได้รูปการวางตัวของจุดหลักฐาน 3 จุดและจุดที่ต้องการทราบค่าคือจุด P น่าเสียดายที่ผมไม่สามารถวาดมุมแสดง dimension เพราะข้อจำกัดของเครื่องคิดเลข และต้องเตือนกันอีกนิดว่าจุด 4 จุดนี้ต้องไม่อยู่บนวงกลมเดียวกัน (โอกาสเกิดน้อยมากๆครับ) ไม่ว่าสูตรไหนก็ไม่สามารถคำนวณได้ เรียกภาวะอย่างนี้ว่าเกิดเหตุการณ์สภาวะเอกฐาน Singularity

 

สำหรับโปรแกรมในชุดเครื่องคิดเลข Casio fx-9860G II SD คงจะมีอีกหลายตอนครับ สำหรับเครดิตของโปรแกรมในชุดนี้ก็กดคีย์ F5 – Info จากเมนูหลัก ก็ติดตามกันตอนต่อไปครับ