Category: Windows

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

การเปลี่ยนแปลงครั้งใหญ่ของไลบรารี 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 ที่ผมเห็นว่าน่าจะเป็นไลบรารีที่ดีที่สุดในด้านการแปลงพิกัดจากระบบพิกัดที่ต่างกันในโลกนี้ จะเห็นว่าการใช้งานไม่ได้ยากเย็นนัก ติดตามกันตอนต่อไปครับ

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 เป็นค่าปริยาย ผมจะออกแบบพัฒนาโปรแกรมให้ผู้ใช้สามารถเลือกโมเดลจีออยด์ได้ ติดตามกันต่อไปครับ

การเล็งสกัดย้อน (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 ต่อนะครับ

การออกแบบเส้นโครงแผนที่ความเพี้ยนต่ำ (Low Distortion Projection) ตอนที่ 2 (กรณีศึกษาออกแบบเส้นโครงแผนที่ความเพี้ยนต่ำสำหรับกรุงเทพมหานครและปริมณฑล)

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

แนะนำการใช้เส้นโครงแผนที่ความเพี้ยนต่ำ (Low Distortion Projection)

และ

การออกแบบเส้นโครงแผนที่ความเพี้ยนต่ำ (Low Distortion Projection) ตอนที่ 1

เรื่องเส้นโครงแผนที่ความเพี้ยนต่ำเป็นเรื่องใหม่สำหรับประเทศไทย แต่ในต่างประเทศบางประเทศได้ประยุกต์ใช้งานมานานแล้ว ประโยชน์ของเส้นโครงแผนที่ความเพี้ยนต่ำเมื่อประยุกต์ใช้แล้วคือ ความต่างระหว่าง Ground Distance และ Grid Distance จะน้อยมากจนสามารถละเลยไปได้ ไม่เหมือนกับการใช้แผนที่ระบบพิกัดยูทีเอ็ม (UTM) ที่ค่าระยะทางบนพื้นโลกกับระยะทางบนแผนที่ต่างกันมาก (ตัวอย่างระยะทางประมาณ 1 กม. สองระยะทางนี้อาจจะต่างกันประมาณ 40-80 ซม.แล้วแต่พื้นที่) แต่ข้อเสียคือจะต้องมีการกำหนดใช้เส้นโครงแผนที่ความเพี้ยนต่ำแบ่งเป็นพื้นที่หรือเป็นโซน ที่ค่าพิกัดศูนย์กำเนิดจะต่างกันไป อาจจะทำให้ช่างสำรวจหรือผู้ใช้งานสับสนได้ แต่ข้อเสียนี้สามารถลดลงได้ถ้ารัฐหรือหน่วยงานของรัฐได้กำหนดและประกาศใช้เป็นทางการ โดยที่มีเอกสารและไฟล์ projection สำหรับแปลงพิกัดจากระบบพิกัด UTM ไปยังระบบพิกัดที่ใช้ LDP ในแต่ละโซน ผู้ใช้งานสามารถนำค่าพารามิเตอร์นี้หรือนำไฟล์ projection (ตัวอย่างเช่นไฟล์ prj ของ Shape file) ไปแปลงพิกัดได้บนโปรแกรมด้าน GIS หรือนำไปตั้งค่าบนเครื่องมืออุปกรณ์เช่น GNSS RTK ที่สามารถแปลงพิกัดได้แบบ real time

ตัวอย่างการประยุกต์ใช้งาน

ผมขอยกตัวอย่างอีกครั้งเช่นรัฐโอเรกอนของอเมริกาที่มีการออกแบบ LDP และประกาศใช้กันมานานแล้วดังรูปด้านล่าง 

พื้นที่รัฐโอเรกอน ประมาณครึ่งหนึ่งของประเทศไทย (ประมาณ 255,000 ตร.กม.)

ออกแบบเส้นโครงแผนที่ความเพี้ยนต่ำสำหรับกรุงเทพมหานครและปริมณฑล

ก็เป็นกรณีศึกษาก็แล้วกันนะครับ ผมจะออกแบบคร่าวๆให้พอมองเห็นภาพในภาพรวม ผมจะไล่ไปตามขั้นตอนที่ได้กล่าวไว้ในตอนที่ 1 และผมจะตั้งเป้าว่า ความเพี้ยน (Distortion) ไม่เกิน 20 ppm ก็มาดูกันว่าในพื้นที่ศึกษานี้ ค่าความเพี้ยนจะอยู่ในเกณฑ์นี้ไหม 20 ppm ก็คือระยะทางจริงๆบนพื้นโลก (Ground Distance)  1 กม. ระยะทางบนระนาบเส้นโครงแผนที่ LDP (Grid Distance) จะต่างกันไม่เกิน 20 มม. (20 มม. ต่อ 1 ล้านมิลมิเมตร หรือ 1 กม. นั่นเอง)

1.กำหนดพื้นที่ขอบเขตและหาค่าตัวแทนความสูงเฉลี่ยเหนือทรงรี (h0)

สำหรับขอบเขตก็ตามหัวข้อคือประกอบไปด้วยจังหวัดกรุงเทพมหานคร สมุทรปราการ นนทบุรี และปทุมธานี ขนาดพื้นที่ประมาณ 85 กม.ในแนวเหนือใต้ และกว้างประมาณ 75 กม. ในแนวตะวันออกตะวันตก หรือกล่าวโดยย่อพื้นที่ 85 กม. x 75 กม.

ต่อไปจะหาค่าระดับที่เป็นตัวแทนความสูงเฉลี่ยเหนือทรงรี (h0) ข้อมูลที่จะนำมาในการหาค่าเฉลี่ยจะใช้แผนที่ของกรมแผนที่ทหาร ปี 2553 ชื่อ “แผนที่แสดงค่าหมุดระดับในเขตกรุงเทพมหานครและปริมณฑล” เนื่องจากแผนที่ไม่สามารถหาแหล่งดาวน์โหลดทางการได้ จึงได้ดาวน์โหลดจากกระดานสนทนาจากเว็บไซต์ ที่ความคมชัดน้อย บางครั้งตัวเลขค่าระดับอาจจะแตกต่างค่าจริงไปบ้าง แต่ผมคิดว่าคงไม่ได้ทำให้การออกแบบ LDP กรณีศึกษานี้มีความด้อยลง  ผมนำแผนที่ชุดนี้มา ทำ rubber sheet เพื่อขึงพิกัดให้เข้ากับเส้นโครงแผนที่ UTM จากนั้นทำการ digitize จุดแต่ละจุดระดับลง ไม่ได้เอาทุกจุด แต่เลือกจุดประมาณ 10 กม.ต่อหนึ่งจุด ค่าระดับนี้เป็นค่าระดับน้ำทะเลปานกลาง (Orthometric Height) ซึ่งเราจะแปลงค่าระดับนี้ไปเป็นค่าระดับเหนือทรงรี (Ellipsoid Height) ในขั้นตอนต่อไป

จากนั้นทำการจัดเก็บจุดค่าระดับเป็นไฟล์ shape file กำหนดระบบพิกัดเป็นภูมิศาสตร์ (Geographic) เพื่อสะดวกต่อการใช้ค่าพิกัดนี้ในภายหลัง

นำไฟล์รูปที่ขึงแล้วและไฟล์จุดค่าระดับเข้าโปรแกรม QGIS ใช้ฟังก์ชั่น vector ทำการหา Basic Statistics for fields จำนวนจุดทั้งหมด 365 จุด ค่าระดับต่ำสุด 0.000 เมตร ค่าระดับสูงที่สุด  9.956 เมตร ค่าเฉลี่ย Mean 2.988 เมตร ผมจะนำค่าเฉลี่ยนี้ไปใช้งาน ค่านี้ขอใช้ตัวย่อเป็น H0 = 2.988 เมตร

ค่าระดับ H0 = 2.988 เมตร นี้จะนำมาแปลงเป็นความสูงเทียบกับทรงรี (h0) การประยุกต์ใช้ LDP ก็คือนำระนาบมาวางแตะค่าระดับนี้ โดยที่กำหนดโซนความกว้างทางราบ และช่วงค่าระดับความสูงที่ยังสามารถใช้ได้

ไดอะแกรมแสดงเส้นโครงแผนที่ความเพี้ยนต่ำที่ระนาบพิกัดฉากสัมผัสที่ความสูงเฉลี่ย h0

2.เลือกเส้นโครงแผนที่และกำหนด Central Meridian ที่จุดใกล้จุดศูนย์กลางพี้นที่

เส้นโครงแผนที่ที่นิยมนำมาทำ LDP มี 3 ประเภทคือ Transverse Mercator (TM), Lambert Conformal Conic (LCC 1SP) และ Oblique Mercator (OM) โดยที่แนวทางการเลือกถ้าพื้นที่ยาวไปในทิศทางตะวันออกตะวันตกเลือก LCC ถ้าพื้นที่ยาวไปในทิศทางเหนือใต้เลือก TM ถ้าพื้นที่เอียงไปในแนวทะแยงมุมกันทิศเหนือใต้และตะวันออกตะวันตกให้เลือก OM ในเคสนี้พื้นที่ยาวในทิศทางเหนือใต้ก็เลือกเป็น TM ที่เราคุ้นเคยกันดี

ต่อไปกำหนด Central Meridian (CM) ที่จุดกึ่งกลางของพื้นที่ (Centroid) ผมไปดาวน์โหลดไฟล์ shape file ที่รวมเอาเส้นขอบเขตของจังหวัดในประเทศไทย เมื่อโหลดมาแล้วเปิดใน QGIS จากนั้นทำการรวมพื้นที่ 4 จังหวัดนี้แบบ Dissolve เพื่อให้เหลือเส้น polygon เส้นเดียว แล้วจะนำไปหาจุดศูนย์กลางพี้นที่ โดยใช้ฟังก์ชั่นด้าน vector ของ Geometry tools เพื่อหา centroid ได้จุดมาดังรูปด้านล่าง

ได้ค่าพิกัดภูมิศาสตร์ของจุด Centroid ดังนี้ latitude: 13.852166 longitude: 100.629706 หน่วยเป็นดีกรี แปลงเป็นหน่วย DMS ได้ latitude: 31°51’7.8″N longitude: 100°37’46.94″E การวาง central meridian จะนิยมเลือกลิปดา (second) ที่เป็นจำนวนเต็ม ผมเลือก ค่าเต็มๆคือ Latitude = 31°51′ และ  CM =  100°38′

3.คำนวณหาค่าสเกลแฟคเตอร์ k0 ที่แกน Central Meridian

ก่อนหน้านี้ในข้อ 1. ผมได้ค่าระดับน้ำทะเลปานกลางของพื้นที่เฉลี่ย (H0) 2.988 เมตร จะนำค่านี้ไปแปลงเป็นค่าระดับเทียบทรงรี (h0) เตือนกันนิดว่าทรงรีที่เราใช้เป็น WGS84 การหา k0 ไม่ได้ยากตามสูตรนี้

Formula 1: Calculate Axis Scale Factor

การหาค่า h0 ก็ไม่ได้ยุ่งยากอะไรในทูลส์ Surveyor Pocket Tools ก็มีโปรแกรม Geoid Height ให้ใช้งาน h0 = H0 + N  โดยที่ N คือ Geoid Separation แต่จะก่อนคำนวณทีละขั้นตอนแบบนี้แบบแมนวลผมจะขอเสนอวิธีที่สะดวกกว่านั้น ผมจะใช้ทูลส์ชือ Init Design LDP ที่อยู่ใน Surveyor Pocket Tools มาช่วย

 คำนวณหาค่า k0 ด้วย Init Design LDP

เปิดโปรแกรม Init Design LDP จะเห็นหน้าตาโปรแกรมดังรูปด้านล่าง

ป้อนข้อมูลค่าระดับเฉลี่ย H0 = 2.988 เมตร ป้อนค่า Latitude  of project center =  13°51′ และ  Longitude of project center (CM )=  100°38′ ที่ได้จากข้อ 2.

จากนั้นคลิกที่ไอคอนลูกศรชี้ลง (เลข 3) เพื่อทำการคำนวณจะได้ผลลัพธ์ดังนี้

โปรแกรมจะคำนวณหาค่า h0 ให้และนำค่านี้ไปแทนในสมการด้านบน สุดท้ายจะคำนวณหา k0 = 0.999996 (แนะนำให้ใช้ทศนิยม 6 ตำแหน่ง) ผมพยายามขยับ CM ไปทางตะวันออกและตะวันตกครั้ง 15″ แต่ค่า k0 ยังเกาะที่ค่า 0.999996 นี้ ผมเลยเลือก CM = 100°38′

4.ตรวจสอบความเพี้ยน (Distortion) ตลอดทั้งพื้นที่

เป็นขั้นตอนที่สำคัญมาก คือถ้าเราเลือก Central Meridian มาหลายๆอันจะต้องเอาค่า k0 มาคำนวณหา Distortion

ค่า k คือ grid scale factor ค่า k นี้เราสามารถหาได้จาก สูตรด้านล่าง (เครดิตจาก Map Projection – A Working Manual ของ John P. Snyder หน้า 61)

Formula 2: Calculate Grid Scale Factor

สูตรก็เป็นสูตรเดียวที่เราใช้หา grid scale factor สำหรับ UTM เพียงแต่ค่า k0 ที่เราใช้จะเป็นค่า k0 ที่ได้จากการคำนวณด้วยโปรแกรม “Init Design LDP” คือ k0 = 0.99996, แทนค่า λ0 ด้วย 1.756383004 เรเดียน (λ0 คือ Central Meridian = 100°38′)

เลือกจุดทดสอบ

ผมเลือกจุดมาทั้งหมด 10 จุด ตำแหน่งให้อยู่ขอบๆ เป็นที่ทราบกันดี ว่าถ้าเป็นเส้นโครงแผนที่ TM ตัว grid scale factor จะเปลี่ยนจากด้านตะวันออก-ตะวันตกเท่านั้น (ไม่มีผลกับเลื่อนไปทางเหนือ-ใต้)

เรากำหนด CM ค่อนข้างจะกลางของพื้นที่ ดังนั้นจะมาดูกันว่าด้านขอบนั้นมี distortion จะยังอยู่ในเกณฑ์ไหม

การตรวจสอบจุดหาความเพี้ยนจะอาศัยการคำนวณจากสูตรที่ผมลงมาให้ก่อนหน้านี้ ลองมา workshop กันดูสักนิด  เริ่มจากจุดที่ 1  แต่จุดต่อๆไปผมจะใช้ทูลส์อีกตัวในชุด Surveyor Pocket Tools มาช่วย

Point No 1 Lat (ɸ)= 13.5079570 Long (λ) = 100.8674784 H = 6.394

คำนวณหาความสูงเทียบกับทรงรี

ได้ค่าh = H + N แทนค่าในสูตร h = 6.394 -29.9632 = -23.5692 เมตร

คำนวณ RG

RGเป็นค่าที่ขึ้นอยู่กับ latitude และพารามิเตอร์ของทรงรี (Ellipsoid) WGS84 ดังนี้

a = 6378137, e = 0.08181919084262149, e’ = 0.08209443794969568

แทนค่า

ใช้สูตรที่ 1 (formula 1) ได้ค่า  RG= 6359074.928

คำนวณค่า k

ใช้สูตรที่ 2 (formula 2) คำนวณได้ค่า T = 0.057708361, C = 0.006371791, A = 0.003973557 สุดท้ายคำนวณหาค่า k = 1.000003945

คำนวณหาค่าความเพี้ยน (Distortion)

แทนค่า k, RG และค่า h ลงไป

จะได้ค่า 7.65 x 10-6 เขียนให้ง่ายคือ 7.65 ppm (7.65 มม. ต่อหนึ่งล้านมม. ซึ่งก็เท่ากับ 7.65 มม.ต่อ 1 กม.) สรุปได้ว่าที่จุดที่ 1

Point No 1 Lat (ɸ)= 13.5079570 Long (λ) = 100.8674784 H = 6.394

มีค่าความเพี้ยน (Distortion) = 7.65 ppm จะเห็นว่ายังไม่เกินค่า 20 ppm ที่ผมตั้ง tolerance ไว้

สร้างเส้นโครงแผนที่ความเพี้ยนต่ำ

มาลองสร้างเส้นโครงแผนที่ความเพี้ยนต่ำโดยอาศัยทูลส์ Create LDP มาช่วย ทูลส์ตัวนี้นอกจากจะสร้าง LDP ได้แล้วยังสามารถตรวจเช็คค่าความเพี้ยนได้ทันที พร้อมทั้งแปลงค่าพิกัดจาก latitude/longitude ไปยังค่าพิกัดใน LDP ได้ เปิดทูลส์ Surveyor Pocket Tools คลิก Create LDP

ตามรูปบนผมสร้างเส้นโครงแผนที่ความเพี้ยนต่ำ โดยอันดับแรกเลือก Projection ก็คือ Transverse Mercator กำหนดใช้ Latitude of origin = 13°51′ และ Central Meridian = 100°38′ ผมกำหนดค่า False Northing (FN) = 500000 และ False Easting (FE) = 200000 หน่วยเป็นเมตร รับรองว่าค่าพิกัดขอบของพื้นที่จะไม่มีค่าติดลบ (พื้นที่ 85 กม. x 75 กม. หรือ 85000 ม. x 75000 ม.) และที่สำคัญมากคือค่า Scale factor at grid origin (k0) = 0.999996 ที่คำนวณไว้ตั้งแต่ตอนแรกๆ ค่า FN และ FE ผมหลีกเลี่ยงเลือกค่าที่ใกล้เคียงกัน และพยายามจะไม่ให้ค่ามากจนไปใกล้เคียงกับ UTM อันจะก่อให้เกิดความสับสน

คำนวณหาค่าความเพี้ยนด้วยทูลส์ Create LDP

ลองป้อนพิกัดจุดที่ 1 เข้าไปในกรอบที่ 2 ดังรูป

ทำการคำนวณด้วยการคลิกไอคอนรูปลูกศรจะได้ผลลัพธ์

ค่าความเพี้ยน 7.65 ppm ตรงกับที่เราคำนวณด้วยมือ ผมใช้โปรแกรมช่วยคำนวณมาทั้ง 10 จุดได้ผลลัพธ์ดังนี้

จะเห็นว่าจุด ที่ 3 มีความเพี้ยน (distortion) อยู่ที่ 13.18 ppm เพราะอยู่ชายขอบด้านตะวันออกสุด และจุดที่ 9 ค่าความเพี้ยนมากถึง 17.27 ppm อยู่เกือบขอบด้านตะวันตก แต่ยังไงก็ไม่เกิน 20 ppm ตัวเลขที่ตั้งไว้ ส่วนที่ลองจิจูดใกล้เคียงกับ CM ได้แก่จุดที่ 4 จะเห็นค่าความเพี้ยนเล็กมาก ขนาด -0.6 ppm แค่นั้นเอง

5.กำหนดพารามิเตอร์เส้นโครงแผนที่ความเพี้ยนต่ำให้เรียบง่าย

จากที่คำนวณและออกแบบมาตั้งแต่ต้น สามารถกำหนดพารามิเตอร์เส้นโครงแผนที่นี้ให้เรียบง่ายและอ่านง่ายได้ดังนี้
Projection: Transverse Mercator
Latitude of grid origin: 13° 51’ 00” N
Longitude of central meridian: 100° 48’ 00” E
Northing at grid origin: 500,000 m
Easting at central meridian: 200,000 m
Scale factor on central meridian: 0.999996 (exact)

ค่าพารามิเตอร์นี้ต้องติดไว้ข้างแผนที่ที่ใช้เส้นโครงแผนที่ความเพี้ยนต่ำนี้เสมอ

6.กำหนดหน่วยระยะทางและพื้นหลักฐานให้ชัดเจน

Linear unit:  Meter

Ellipsoidal datum :  World Geodetic System 1984 (WGS84)

Vertical datum:  Mean Sea Level (MSL)

System:  Bangkok Metropolis Low Distortion Projection Coordinate System

Zone:  Bangkok Metropolis Area

หน่วยระยะทางและพื้นหลักฐานต้องติดไว้ข้างแผนที่ที่ใช้เส้นโครงแผนที่ความเพี้ยนต่ำนี้เสมอ

จัดเก็บเส้นโครงแผนที่ความเพี้ยนต่ำเข้าฐานข้อมูล

ทูลส์ Create LDP นอกจากสามารถคำนวณหาค่าความเพี้ยนและแปลงพิกัดได้แล้ว ยังสามารถจัดเก็บเส้นโครงแผนที่ที่เราออกแบบ เข้าไปเก็บในฐานข้อมูล เพื่อความสะดวกสามารถนำมาใช้ในภายหลังได้ เมื่อเปิดโปรแกรม Surveyor Pocket Tools

คลิกที่ LDP Database จะเห็นฐานข้อมูลของ LDP สำหรับเครื่องผมแล้วมีฐานข้อมูลเก็บเส้นโครงแผนที่ความเพี้ยนต่ำที่ผมนำมาศึกษาดังนี้

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

จากนั้นกลับมาดูที่ LDP Database คลิกที่ไอคอนลูกศรวนเพื่อ “Refresh” จะเห็นฐานข้อมูลอัพเดทดังนี้ ผมวงสีแดงเน้นให้ดูพารามิเตอร์ที่ป้อนไป

แปลงพิกัดด้วยทูลส์ Transform Coordinates

เมื่อจัดเก็บฐานข้อมูลเส้นโครงแผนที่ความเพี้ยนต่ำแล้ว ต้องการแปลงพิกัดระหว่าง UTM/Geographic ไปยังเส้นโครงแผนที่ความเพี้ยนต่ำก็สามารถทำได้ดังตัวอย่างต่อไปนี้ เริ่มจากเปิดโปรแกรม Transform Coordinates มาก่อน

จากตัวอย่างคำนวณค่าพิกัดยูทีเอ็ม N = 1,561,926.0937, E = 639,807.9239 Zone: 47N บนพื้นหลักฐาน WGS84 ไปยังพื้นหลักฐาน Bangkok Metropolis (LDP) จะได้ค่าพิกัด N = 530,441.5826, E = 163,493.4658

เครดิตสูตรที่นำมาใช้

เส้นโครงแผนที่ความเพี้ยนต่ำจริงๆแล้วก็คือเส้นโครงแผนที่ตัวหนึ่ง ไม่มีอะไรพิเศษพิศดาร เพียงแต่ยกระนาบขึ้นมาแตะทึ่ค่าระดับเฉลี่ย h0 จุดนี้เองที่พิเศษเพราะว่าจะได้ค่า k0 ตัวใหม่ที่ไม่ใช่ 0.9996 แบบยูทีเอ็ม ดังนั้นสูตรที่นำมาใช้เพื่อคำนวณเส้นโครงแผนที่ความเพี้ยนต่ำ ทั้งการหา Grid Scale Factor หรือแปลงพิกัด ก็ยังเป็นสูตรของ Transverse Mercator

ในโปรแกรม Surveyor Pocket Tools ผม implement สูตรการคำนวณจาก Map Projection – A Working Manual by John P. Snyder เป็นโค้ดภาษาไพทอนเฉพาะการคำนวณที่เกี่ยวข้องกับ LDP ส่วนการแปลงพิกัดข้ามพื้นหลักฐานอื่นๆยังใช้ไลบรารีจาก Proj4 ถ้าสนใจสูตรของ Transverse Mercator ให้ไปดูได้ที่หน้า 60-63 ของตำราเล่มนี้

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

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

 

เมื่อผมลองดี โมดิฟายด์แอพ DJI GO 4 ด้วย Deejayeye-Modder สำหรับโดรน DJI Spark

จากลองของเป็นลองดี

จากตอนที่แล้วเมื่อต้องลองของเอาโดรนเซลฟี่ DJI Spark มาบินทำแผนที่ภาพถ่ายทางอากาศ ตอนนี้เลยคำว่า”ลองของ”ไปไกลหลายช่วงตัวแล้ว เพราะแพ็ตช์หรือโมดิฟายด์แอพ DJI GO 4 เรียกง่ายว่าเป็นการ “ลองดี” ไปแล้ว เพื่อนำโดรนมาบินทำแผนที่ภาพถ่ายทางอากาศโดยที่เปิดหรือใช้ฟีเจอร์ที่ซ่อนไว้ ตอนนี้ฟีเจอร์เทียบได้กับโดรนรุ่นใหญ่เช่น Phantom 4 Pro หรือ Mavic

คำเตือน

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

ต้นทุนความเสี่ยง

มาแจกแจงราคาของสปาร์คกันนิด ผมซื้อสปาร์คแบบ fly combo ราคาในตอนนั้น 23,900 บาท ซื้อแบตเตอรีมาเพิ่มหนึ่งลูกราคา 2000 บาทถ้วน (รวมแบตเตอรีทั้งหมด 3 ลูก) ซื้อกระเป๋าสำหรับใส่โดรนมาครั้งแรก 650 บาทแต่ไม่ถูกใจเพราะใส่ของได้ไม่เยอะ ซื้อมาอีกใบราคา 1850 บาท สาย OTG 200 บาท รวมทั้งหมด 28,600 บาท ใช้เงินขนาดนี้สามารถซื้อสมาร์ทโฟนรุ่นดีๆ ได้สักหนึ่งเครื่องเลยทีเดียว

ทำไมต้องโมดีฟายด์แอพ

ความจริงโดรนแต่ละรุ่นของ DJI แอพจะเป็นส่วนหนึ่งที่เป็นตัวกำหนดว่าตัวไหนมีความสามารถอะไรบ้าง โดรนแต่ละรุ่นจะถูกจัดให้ Product อยู่ในเกรดไหนระดับไหน ตามกำลังซื้อของลูกค้า ลูกค้าที่ต้องการ Premium พันธุ์แท้มีกำลังซื้อสูงก็ควรจะมาซื้อรุ่นเทพ Inspire ที่ราคาเหยียบสองแสนบาทไปใช้ ส่วนใครมีกำลังซื้อลดหลั่นกันมาก็ไล่ตั้งแต่ Phantom 4 Pro, Mavic Pro, Mavic Air จนกระทั่ง DJI Spark โดรนเซลฟี่ราคาถูกที่สุดในไลน์การผลิดของ DJI

ฟีเจอร์ของ DJI Spark ที่ถูกซ่อนคือบินตามเส้นทางที่กำหนด (Waypoint) กับบินเป็นวงโคจร (Orbit) ก่อนหน้านี้ผมพยายามเข้าไปดูข่าวสารตามที่ DJI เคยบอกไว้ว่าอาจจะเอาการบินตามเส้นทางที่กำหนดมาใส่ให้สปาร์ค แต่ก็ตามข่าวมาสี่ ห้าเดือนแล้วก็ไม่มีความคืบหน้าอะไร และในตอนนี้ที่เขียนบทความนี้ ทางดีเจไอก็ออกโดรนรุ่น Mavic Air ตัวใหญ่กว่าสปาร์คนิดหนึ่ง ราคาแพงกว่าประมาณหมื่นกว่าบาท บินได้นานขึ้น gimbal ได้สามแกน (แต่สำหรับคนทำแผนที่ภาพถ่ายทางอากาศสองแกนก็พอครับ) แต่ก่อนจะขายบอกว่าสนับสนุนการบินด้วย waypoint แต่พอออกขายจริงๆ กลับตัดฟีเจอร์นี้ไปดื้อๆ ซะอย่างนั้น ก็ได้ก้อนอิฐไปตามระเบียบไปหลายก้อน

นักโมกลุ่ม Deejayeye-Modder

ผมเข้าไปติดตามข่าวสารบ่อยๆตามฟอรั่มของผู้ใช้โดรน ก็เลยค้นไปเจอว่ามีกลุ่มโมดิฟายด์ (modder) ชื่อรวมๆคือ Deejayeye โมดิฟาย์แอพ DJI Go 4 เพื่อเปิดฟีเจอร์ที่ซ่อนนี้สำหรับโดรนสปาร์ค เป็นโค้ดเปิดใน Github โดยการเอาไฟล์ที่ใช้สำหรับติดตั้งในแอนดรอยด์คือไฟล์ APK (Android Package) มาทำการ patch เพื่อเปิดฟีเจอร์ที่ถูกซ่อนนี้ ในตอนแรกก็ไม่มั่นใจแต่ลองเอามาทำด้วยเครื่องมือที่ทีมงานแนะนำมาพร้อมโค้ด ก็ลองดูกลับทำได้สำเร็จ ลองเอาไฟล์นี้มาติดตั้งในโทรศัพท์แอนดรอยด์ของผมก็สามารถใช้งานได้

ย้อนรอยตอนลองของ

ขอเล่าย้อนรอยสั้นๆก่อนหน้านี้ ผมใช้แอพ Litchi ผมพอใจการทำงานของ Litchi ปัญหาเมื่อโดรนต่อกับรีโมทคอนโทรล (Remote controller) สัญญานหลุดจากกันแต่ก็ reconnect ให้ใหม่ได้เร็ว ไม่ต้องลุ้นมาก แอพมีความเสถียรพอสมควร นานๆจะแครชที แต่ข้อเสียคือสุ่มเสี่ยงกับการไม่ได้รับประกันจาก DJI และข้อเสียอีกอย่างคือไม่สนับสนุนการบินตามเส้นทางที่กำหนด (waypoint) สำหรับโดรน DJI Spark ที่ไม่สนับสนุนเพราะว่าเครื่องมือ SDK (Software Development Kit) ของ DJI ไม่มีไว้ให้นั่นเอง คำว่าไม่สนับสนุนคือสั่งให้บินไปตามเส้นทางไม่ได้ แต่สามารถสร้างเส้น waypoint นี้ได้

ผมกำหนดเส้นทางการบิน (waypoint) จากที่มีในแอพ เรียกว่าสามารถใส่จุดล่วงหน้าจากแอพโดยการจิ้มแต่ละจุดไปบน map แต่ข้อเสียคือไม่สามารถนำเข้าไฟล์จาก kml ได้ เครื่องมือสร้าง auto waypoint จากรูปปิดที่กำหนดให้ก็ไม่มี ตอนบินก็เรียกเส้นทางการบินนี้ให้แสดงบนแผนที่ จากนั้นก็บังคับโดรนบินไปตามเส้นทางเหล่านี้ คุณภาพของการบินให้ตรงกับเส้นทางที่กำหนดนี่ต้องอาศัยทักษะพอสมควร แต่ส่วนใหญ่ผมเป๋ไปด้านซ้ายด้านขวาบ้าง แต่ผมออกแบบให้เส้นแนวบิน Overlap กันมากขึ้นคือให้เส้นแนวบินเข้ามาใกล้กันมากกว่าทฤษฎีจะได้ชดเชยเมื่อบินไม่ตรงแนว ตอนบินก็สั่งให้ gimbal ก้มด้วยมุม 85 องศา (เต็มที่แล้ว ไม่ได้ 90 องศา) จากนั้นก็สั่งถ่ายรูปด้วย interval ทุกๆ 3 วินาที ลองดูรูปด้านล่างจะเห็นว่าแนวบินไม่เป็นเส้นตรงเฉไปเฉมา

ตั้งแต่ผมซื้อโดรนมา เคยเอามาเซลฟี่ประมาณ 4-5 ครั้งแค่นั้นเอง จากนั้นเอามาบินทำแผนที่ภาพถ่ายทางอากาศทั้งหมด ดังนั้นฟีเจอร์ที่ผมต้องการที่สุดก็คือบินตามเส้นทางที่กำหนดเท่านั้นเพราะจะได้เบาแรงเมื่อเอาโดรนขึ้นไปบนฟ้า เพิ่มเติมอีกนิดผมเคยใช้แอพทางการของ DJI GO 4 อยู่ไม่กี่ครั้ง จากนั้นก็ใช้ Litchi มาตลอด

เมื่อ DJI กลับลำในกรณี Mavic Air

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

พบแสงสว่างบนทางแห่งความเสี่ยง

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

แพ็คเกจ APK สามารถหาดาวน์โหลดได้เฉพาะ DJI GO 4 รุ่น 4.1.4 ถึง 4.1.9 แต่รุ่นหลังจากนี้ทาง DJI ได้ทำการ encrypted ไว้ ไม่สามารถหาได้บนดิน ต้องมุดใต้ดินไปหาครับ

ปัญหา OTG (On The Go)

ปัญหาของ OTG บนแอนดรอยด์สำหรับใช้บนแอพ DJI GO 4 ถือว่าเป็นปัญหาอลเวง เดี๋ยว DJI สนับสนุนเดี๋ยวถอดออก ผมซื้อสายมา 200 บาทแต่ไม่ได้ใช้เลย แต่ที่ได้ยินมาคือฝั่ง IOS ใช้ได้ดีไม่ปัญหาเหมือนฝั่งแอนดรอยด์

OTG คือสายเคเบิลที่เอามาต่อระหว่างรีโมทคอนโทรลกับโทรศัพท์มือถือที่ลงแอพ DJI GO 4 คือที่ได้ยินมาตลอด ถ้ามันต่อได้ดีเสถียร จะเป็นประโยชน์มากกว่าการต่อด้วยไวไฟเพราะไวไฟนี้อาจจะไปกวนไวไฟที่ต่อระหว่างรีโมทคอนโทรลกับโดรนที่บินกลางอากาศ

เลือกเวอร์ชั่นของแอพ DJI GO 4

ผมหาแพ็คเกจ apk ของแอพตอนแรกเลือก 4.1.22 รุ่นล่าสุดในขณะที่เขียนบทความ ทำการแพ็ตช์ ติดตั้งลงโทรศัพท์มือถือปรากฎว่าสาย OTG ใช้ไม่ได้ครับ เพราะตอนแรกอ่านผ่านๆว่ารุ่นหลังๆสนับสนุนแล้ว พอไปอ่านดีๆกลับผมว่า DJI ไม่สนับสนุนในแอพรุ่นนี้ สุดท้ายพบว่ารุ่น 4.1.15 ที่สนับสนุน ก็เลยต้องไปหาแพ็คเกจ apk รุ่นนี้มาทำการแพ็ตช์ แล้วทำการติดตั้งลงมือถือ เลือก Flight mode คือปิดการสื่อสารทุกอย่างบนโทรศัพท์มือถือ เปิดแอพ DJI GO 4 รุ่น 4.1.15 ในที่สุด OTG ใช้ได้

เปิดฟีเจอร์ Waypoint

เมื่อเปิดแอพ DJI GO 4 ทำการต่อสายรีโมทคอนโทรล จากนั้นก็ต่อกับโดรน ช่่วงนี้ลองในบ้านได้ ทำการปรับค่า config ต่างๆ แล้วลองแท็บไปในโหมด Fly intelligent จะเห็นว่า การบินด้วยเส้นทางที่กำหนด (Waypoints) เปิดมาเป็นที่เรียบร้อยพร้อมที่จะไปลองทดสอบ

ปัญหาเมื่อเส้นทางกำหนดการบินไม่สามารถกำหนดล่วงหน้าได้

ผมได้บอกไปแล้วว่าด้าน hardware โดรนในค่ายของ DJI ทำได้ยอดเยี่ยม แต่ด้านซอฟท์แวร์ กลับได้ยินเสียงอึงคะนึงหลายๆอย่าง ที่หนักๆแรงๆคือการบินด้วย waypoint ไม่สามารถป้อนหรือใส่ได้ก่อนแบบที่ Litchi หรือแอพตัวอื่นทำได้ จะวางแผนการบินล่วงหน้าทาง DJI ก็มีให้คือเสียเงินไปซื้อ DJI Ground Station สำหรับ IOS แต่ไม่มีสำหรับฝั่งแอนดรอยด์ที่ผมใช้อยู่

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

ปัญหาอีกอย่างคือในการบินตามเส้นทางที่กำหนด ต้องทำการเป็นลำดับ 1, 2, 3, 4, …. (ไม่ต้องถึงไลน์สุดท้ายก็สามารถเรียกโดรนกลับมาได้) แต่ตัวอย่างเช่นต้องการบินจาก 3, 4 ,5, 6 ทำไม่ได้ครับ ต้องไล่ตั้งแต่ 1, 2, 3 ซึ่งเป็นข้อเสียอย่างแรง

ตามล่าหาวิธีการสร้างเส้นทางการบินที่กำหนดล่วงหน้า

จากปัญหาที่พบผมไม่ต้องไปค้นที่ไหนไกล ในกลุ่มนักโมดิฟายด์ ก็มีคนเห็นปัญหานี้ เขาใช้เครื่องที่ root ตามเข้าไปในดูในโฟลเดอร์ของแอพบนมือถือก็พบว่ามีไฟล์ที่เป็นฐานข้อมูล sqlite ข้างในมีตารางที่แสดง waypoint ด้วยสตริงแบบ json จากนั้นก็เพิ่มฟีเจอร์เข้าไปใน DJI GO 4 เพื่อทำให้ไฟล์นี้เป็น public คือเครื่องไม่ root สามารถไปแก้ไขได้ ในขณะนี้รู้ว่าแอพทำการเก็บข้อมูลที่ไหนและรู้รูปแบบคือเป็นฐานข้อมูล sqlite

Mission Planner สุดยอดโปรแกรมออกแบบเส้นทางการบิน

โปรแกรมนี้บน Desktop PC เป็นโปรแกรมฟรี ใช้งานง่ายมาก ไปดาวน์โหลดได้ที่ ลิ๊งค์นี้  โปรแกรมนี้สามารถสร้าง Polygon ด้วยการจุดบน map ได้ง่ายๆ หรือจะนำเข้าจาก shape file ก็ได้ เมื่อได้พื้นที่แล้วสามารถสร้าง waypoint ได้จากการป้อนค่ามุมอะซิมุทแนวการบิน ความสูงของโดรน และค่าอีกหลายๆค่า จากนั้นโปรแกรมจะสร้าง waypoint ให้ ซึ่งเราสามารถจัดเก็บเป็นไฟล์ในรูปแบบรหัสแอสกี้ (text file) ได้ แต่ตอนนี้ได้ไฟล์มาแต่ยังไม่โดนใจแอพ DJI GO 4 ต้องทำการแปลงเป็นฐานข้อมูล sqlite ก่อน

มาลองดูไฟล์ waypoints ที่ได้จากโปรแกรม Mission Planner

QGC WPL 110
0 1 0 16 0 0 0 0 23.864033 90.366674 8.000000 1
1 0 0 16 0.00000000 0.00000000 0.00000000 0.00000000 23.86296860 90.36620300 50.000000 1
2 0 0 16 0.00000000 0.00000000 0.00000000 0.00000000 23.86312430 90.36512590 50.000000 1
3 0 0 16 0.00000000 0.00000000 0.00000000 0.00000000 23.86336720 90.36516850 50.000000 1
4 0 0 16 0.00000000 0.00000000 0.00000000 0.00000000 23.86321230 90.36623980 50.000000 1
5 0 0 16 0.00000000 0.00000000 0.00000000 0.00000000 23.86345610 90.36627670 50.000000 1
6 0 0 16 0.00000000 0.00000000 0.00000000 0.00000000 23.86361010 90.36521110 50.000000 1
7 0 0 16 0.00000000 0.00000000 0.00000000 0.00000000 23.86385300 90.36525370 50.000000 1
8 0 0 16 0.00000000 0.00000000 0.00000000 0.00000000 23.86369980 90.36631360 50.000000 1
9 0 0 16 0.00000000 0.00000000 0.00000000 0.00000000 23.86394360 90.36635050 50.000000 1
10 0 0 16 0.00000000 0.00000000 0.00000000 0.00000000 23.86409600 90.36529640 50.000000 1
11 0 0 16 0.00000000 0.00000000 0.00000000 0.00000000 23.86433890 90.36533900 50.000000 1
12 0 0 16 0.00000000 0.00000000 0.00000000 0.00000000 23.86418730 90.36638740 50.000000 1
13 0 0 16 0.00000000 0.00000000 0.00000000 0.00000000 23.86443110 90.36642420 50.000000 1
14 0 0 16 0.00000000 0.00000000 0.00000000 0.00000000 23.86458180 90.36538160 50.000000 1
15 0 0 16 0.00000000 0.00000000 0.00000000 0.00000000 23.86482470 90.36542420 50.000000 1
16 0 0 16 0.00000000 0.00000000 0.00000000 0.00000000 23.86467480 90.36646110 50.000000 1
17 0 0 16 0.00000000 0.00000000 0.00000000 0.00000000 23.86491860 90.36649800 50.000000 1
18 0 0 16 0.00000000 0.00000000 0.00000000 0.00000000 23.86506770 90.36546680 50.000000 1
19 0 0 16 0.00000000 0.00000000 0.00000000 0.00000000 23.86531060 90.36550940 50.000000 1
20 0 0 16 0.00000000 0.00000000 0.00000000 0.00000000 23.86516230 90.36653490 50.000000 1

ArduMissionToDJISQL ทูลส์มาช่วยแปลงข้อมูล

และนักโมดิฟายด์ท่านนี้ก็ได้สร้างทูลส์แปลงจาก textfile เป็นไฟล์ sql ซึ่งเป็น text file เช่นเดียวกันแต่เป็นไฟล์โค๊ดภาษา sql เพื่อเตรียมปั๊มป์ข้อมูลนี้เข้าไปในฐานข้อมูล ไปหาดาวน์โหลดได้ที่ฟอรั่มตามลิ๊งค์นี้

ไฟล์ sql ที่ได้จากการแปลง waypoints ก็แบบนี้ครับ

INSERT INTO dji_pilot_dji_groundstation_controller_DataMgr_DJIWPCollectionItem ( distance, pointsJsonStr, autoAddFlag, createdDate )
VALUES (1328.54, ‘{“points”:[{“craftYaw”:-81,”gimbalPitch”:-85,”gimbalYaw”:0,”height”:50.00000000000000,”lat”:23.86296860000000,”lng”:90.36620300000000},{“craftYaw”:9,”gimbalPitch”:-85,”gimbalYaw”:0,”height”:50.00000000000000,”lat”:23.86312430000000,”lng”:90.36512590000000},{“craftYaw”:99,”gimbalPitch”:-85,”gimbalYaw”:0,”height”:50.00000000000000,”lat”:23.86336720000000,”lng”:90.36516850000000},{“craftYaw”:8,”gimbalPitch”:-85,”gimbalYaw”:0,”height”:50.00000000000000,”lat”:23.86321230000000,”lng”:90.36623980000000},{“craftYaw”:-81,”gimbalPitch”:-85,”gimbalYaw”:0,”height”:50.00000000000000,”lat”:23.86345610000000,”lng”:90.36627670000000},{“craftYaw”:9,”gimbalPitch”:-85,”gimbalYaw”:0,”height”:50.00000000000000,”lat”:23.86361010000000,”lng”:90.36521110000000},{“craftYaw”:99,”gimbalPitch”:-85,”gimbalYaw”:0,”height”:50.00000000000000,”lat”:23.86385300000000,”lng”:90.36525370000000},{“craftYaw”:8,”gimbalPitch”:-85,”gimbalYaw”:0,”height”:50.00000000000000,”lat”:23.86369980000000,”lng”:90.36631360000000},{“craftYaw”:-81,”gimbalPitch”:-85,”gimbalYaw”:0,”height”:50.00000000000000,”lat”:23.86394360000000,”lng”:90.36635050000000},{“craftYaw”:9,”gimbalPitch”:-85,”gimbalYaw”:0,”height”:50.00000000000000,”lat”:23.86409600000000,”lng”:90.36529640000001},{“craftYaw”:99,”gimbalPitch”:-85,”gimbalYaw”:0,”height”:50.00000000000000,”lat”:23.86433890000000,”lng”:90.36533900000001},{“craftYaw”:8,”gimbalPitch”:-85,”gimbalYaw”:0,”height”:50.00000000000000,”lat”:23.86418730000000,”lng”:90.36638739999999},{“craftYaw”:-81,”gimbalPitch”:-85,”gimbalYaw”:0,”height”:50.00000000000000,”lat”:23.86443110000000,”lng”:90.36642420000000},{“craftYaw”:9,”gimbalPitch”:-85,”gimbalYaw”:0,”height”:50.00000000000000,”lat”:23.86458180000000,”lng”:90.36538160000001},{“craftYaw”:99,”gimbalPitch”:-85,”gimbalYaw”:0,”height”:50.00000000000000,”lat”:23.86482470000000,”lng”:90.36542420000001},{“craftYaw”:8,”gimbalPitch”:-85,”gimbalYaw”:0,”height”:50.00000000000000,”lat”:23.86467480000000,”lng”:90.36646110000000},{“craftYaw”:-81,”gimbalPitch”:-85,”gimbalYaw”:0,”height”:50.00000000000000,”lat”:23.86491860000000,”lng”:90.36649800000001},{“craftYaw”:9,”gimbalPitch”:-85,”gimbalYaw”:0,”height”:50.00000000000000,”lat”:23.86506770000000,”lng”:90.36546679999999},{“craftYaw”:99,”gimbalPitch”:-85,”gimbalYaw”:0,”height”:50.00000000000000,”lat”:23.86531060000000,”lng”:90.36550939999999},{“craftYaw”:99,”gimbalPitch”:-85,”gimbalYaw”:0,”height”:50.00000000000000,”lat”:23.86516230000000,”lng”:90.36653490000000}]}’,1,1524200609000)

SQLite Studio โปรแกรมช่วยบริหารฐานข้อมูล SQLite

SQLite Studio เป็นโปรแกรมฟรีและเปิดโค้ด ผมใช้มานานแล้ว เมื่อติดตั้งแล้ว ผมจะใช้โอนไฟล์จากแอพ DJI GO 4 ชื่อไฟล์คือ dji_mod_4_1_15.db ในเครื่องผมอยู่ที่ /mnt/storage/DJI/ จากนั้น copy มาฝั่ง desktop pc ทำการแก้ไขด้วย SQLite Studio ด้วยการ Add Database แล้วเลือกไฟล์นี้ จากนั้นมองหาตารางข้อมูลชื่อ “dji_pilot_dji_groundstation_controller_DataMgr_DJIWPCollectionItem” เปิด Sqlite Editor จากนั้นเปิดไฟล์ sql  ด้วย Notepad จากขั้นตอนที่แล้วทำการ copy เนื้อหาไปยังคลิปบอร์ด แล้วทำการ paste ที่ Sqlite Editor นี้แล้วทำการ Execute Query ข้อมูลจะถูกปั๊มป์เข้าไปในฐานข้อมูล จากนั้นจะ copy ไฟล์ ฐานข้อมูลนี้กลับไปยังโทรศัพท์มือถือเพื่อไปทับกับไฟล์เดิม

ทดสอบบินตามเส้นทางที่กำหนด

วันแรกที่ผมต้องการลองบินแบบบินตามเส้นทางที่กำหนด (waypoint) รู้สึกใจไม่ค่อยดี หนึ่งนั้นที่เคยบอกไปว่าไม่ค่อยคุ้นกับแอพ DJI GO 4 เพราะไม่ค่อยได้ใช้ และอีกอย่างคือไฟล์เส้นทางการบินที่เตรียมมาไม่สามารถเปิดดูได้จากแอพ DJI GO 4 ก่อน ว่ามีอะไรบ้าง นี่ก็เป็นข้อเสียที่ร้ายแรงอีกอย่างครับ ทำไมให้เปิดดูไม่ได้ก่อน ต้องเอาโดรนบินขึ้นไปบนฟ้าเท่านั้นจึงจะดูได้ แหมถ้าเป็นยุคซามูไร ทีมงานผู้พัฒนาชุดนี้ควรจะต้องคว้านท้องตัวเอง (เอาฮานะครับ) เมื่อขึ้นไปบนท้องฟ้าแล้วผมก็ดูไฟล์แนวบินที่ผมสร้างเอาไว้ ในฐานข้อมูลนั้นจะมี location ด้วยผมใส่คำอธิบายไปด้วยสั้นๆ เพื่อให้หาไฟล์ง่าย เพราะ DJI GO 4 เลือกที่จะเน้นแสดงผลเส้นทางการบินนี้ด้วยวันเวลาที่สร้างซึ่งดูยากมาก

ครั้งแรกที่ลองบินตาม waypoint ผมเลือกเอาเส้นทางการบินตามแนวรถไฟฟ้าที่กำลังก่อสร้างนอกเมือง บริเวณนั้นมีหนองน้ำเยอะพอสมควร สลับด้วยทุ่งหญ้า เมื่อบินขึ้นท้องฟ้าปรากฎว่าผมไปยืนอยู่ผิดตำแหน่งคือห่างจาก waypoint ไปประมาณ 400 เมตร (ถือว่าไกลสำหรับสปาร์ค) ด้วยความที่ไม่อยากเสียเวลาก็เลย Apply เพื่ออัพโหลด waypoint ชุดนี้เข้าโดรน จากนั้นโดรนก็บินไปตาม waypoint จุดที่ 1 พอบินไปได้สัก 6 เส้น ผมลอง stop เพื่อจะเรียกเครื่องกลับ เป็นเรื่องครับ สัญญาณกลับขาดหายไปดื้อๆ ทั้งๆตอนที่บินอยู่ในตอนนั้นโดรนส่งรูปที่ถ่ายมาให้โทรศัพท์มือถือดูได้ตลอด ก่อนสัญญาณจะหลุดผมดูแล้วมีแบตเตอรีเหลือ 38% เมื่อติดต่อกันไม่ได้ผมก็เดินจ้ำอ้าวไปยังจุดที่โดรนบินค้างอยู่กลางอากาศแต่ไม่ทราบว่าอยู่ตรงไหน ด้วยอารามรีบร้อนกลับเลยไปอีกทางหนึ่ง เมื่อตั้งสติได้ก็มาดูที่แอพอีกทีดูจุดสุดท้ายบนแผนที่ ก็เลยวิ่งกันมาที่จุดนั้น เดชะบุญปรากฎว่าโดรนแลนดิ้งลงมาเอง อยู่กลางถนนลาดยางเล็กๆที่ไม่มีใครใช้ แบตเตอรีกระพริบอยู่ พอผมไปถึงแบตเตอรีก็หมดพอดี ผมลองเปลี่ยนแบตเตอรีลูกใหม่ใส่เข้าไปทดลองบินขึ้นอีกครั้งก็ยังดีอยู่ ไม่มีอะไรเสียหาย เมื่อกลับไปผมเอารูปที่ถ่ายมาดู จุดสุดท้ายที่ถ่ายเป็นทุ่งหญ้า แต่ก็สงสัยว่าเครื่องมาแลนดิ้งลงที่ถนนลาดยางห่างออกไป 20 เมตรได้อย่างไร และปลอดภัยด้วย ถือเป็นความโชคดีมากๆเนื่องจากในบริเวณนั้นส่วนใหญ่เป็นหนองน้ำ

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

ไม่พบปัญหาในการทดสอบในภายหลัง

วันต่อมาผมก็เอาโดรนพร้อมทั้งทำ waypoint ประมาณ 3 ชุด แต่ละชุดผมคำนวณระยะทางและความเร็วของการบินโดรนซึ่งผมไม่ให้เกิน 12 กม.ต่อชม. ใช้เวลาประมาณแค่ 9-10 นาทีเท่านั้น (แบตเตอรีตามสเป็คแล้วใช้ได้ 15 นาที) และจุด home ที่ไปยืนปล่อยโดรนต้องเป็นจุดที่อยู่ในพื้นที่นั้นๆ ป้องกันโดรนบินออกไปไกลทั้งขาไปและขากลับ วันที่สองนี้ไม่มีปัญหาครับ สัญญาณไม่หลุด นิ่ง และที่สำคัญคือในขณะนั้นลมแรงไปนิด เพราะมีเมฆฝนตั้งเค้าไกลๆ แต่โดรนสามารถบินไปตามแนวได้อย่างตรงแนว ผิดกับที่ผมบินด้วย manual แบบหน้ามือกับหลังมือ และแอพก็เตือนตลอดว่าลมแรงให้บินด้วยความระมัดระวัง เมื่อบินจบแนวเส้นทางการบินที่กำหนดไว้แล้ว เครื่องก็บินกลับมาจุดปล่อย (home) ได้อย่างปลอดภัย ลองดูรูปด้านล่างจะเห็นว่าแนวบินเป็นเส้นตรง

สรุป

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

แนะนำโปรแกรมมิ่งภาษาซีบนเครื่องคิดเลข Casio fx-9860G II SD ด้วยเครื่องมือพัฒนา SDK ของ Casio

แนะนำโปรแกรมมิ่งภาษาซีบนเครื่องคิดเลข Casio fx-9860G II SD ด้วยเครื่องมือพัฒนา SDK ของ Casio

เคยเกริ่นมาก่อนว่าต้องการเขียนบทความนี้ขึ้นมาเพื่อวงการศึกษาบ้านเราที่สนใจเรื่องโปรแกรมมิ่งบนเครื่องคิดเลขสามารถจะพัฒนาโปรแกรมภาษาซีบน Casio fx-9860G II SD หรือรุ่นที่ใกล้เคียงนี้ได้ โดยที่มีไม่มีข้อจำกัดด้านภาษาโปรแกรมมิ่ง เหมือนกับภาษา casio basic อาจจะส่งผลให้ในอนาคต มีโปรแกรมที่พัฒนาโดยบุคคลากรท่านอื่นๆ เข้ามาสู่วงการนี้มากขึ้น และได้ตัวโปรแกรมงานสำรวจที่มีความหลากหลายและความสามารถมากขึ้นทั้งนี้เพื่อขยายขีดความสามารถโปรแกรมบนเครื่องคิดเลขให้สามารถคิดงานที่ยาก ซับซ้อนได้ บางครั้งเกือบจะเทียบเท่าโปรแกรมที่ใช้งานบนคอมพิวเตอร์

เครื่องมือพัฒนา Software Development Kit (SDK)

เครื่องมือตัวนี้เดิมทีสามารถดาวน์โหลดที่เว็บไซต์ตามลิ๊งค์นี้ได้ http://edu.casio.com/support/en/agreement.html#2 ขั้นตอนแรกยอมรับเงื่อนไขแล้วเลื่อนหน้าไปด้านล่างๆจะเห็นเครื่องคิดเลขรุ่น fx-9860 เมื่อคลิกลิ๊งค์ SDK เข้าไปจะเห็นว่าลิ๊งค์เครื่องมือพัฒนาโปรแกรมขาด แต่คู่มือยังสามารถดาวน์โหลดมาอ่านศึกษาได้ ผมอาศัยลงใต้ดินที่มีคนปล่อยให้ดาวน์โหลด (ถ้าใครอยากได้เครื่องมือตัวนี้ก็ขอมาหลังไมค์กันได้ครับ) และดาวน์โหลดโปรแกรม FA-124 มาด้วยอยู่ในหมวด Support Software/PC Link software

ติดตั้งเครื่องมือพัฒนาโปรแกรม

เปิดไฟล์ zip ของเครื่องมือพัฒนาจะเห็นไฟล์ข้างในดังนี้

จากนั้นทำการติดตั้งเครื่องมือลงคอมพิวเตอร์  ข้อสำคัญคือพาทของโฟลเดอร์หรือไดเรคทอรีที่ติดตั้งจะต้องไม่มีช่องว่าง ดังนั้นให้ติดตั้งไปที่รากของไดรว์ C: ตัวอย่างผมใช้ชื่อว่า fx-9860-sdk 

ถ้าพาทของโฟลเดอร์ที่ติดตั้งมีช่องว่างการ compile & build จะไม่ผ่านเลย เมื่อติดตั้งแล้วจะมีไอคอนที่หน้า desktop

เริ่มต้นใช้งาน

เมื่อเปิดโปรแกรมจากไอคอนที่ desktop จะเห็นหน้าตาเครื่องมือพัฒนาโปรแกรม ดังรูป อาจจะดูทื่อๆเพราะเครื่องมือตัวนี้ออกมานานแล้วตั้งแต่วินโดส์รุ่นก่อนหน้านี้ ที่เมนู “Project” คลิกเลือก “New” ผมสร้างโฟลเดอร์ชื่อ “FX9860GIISD” ไว้ที่ไดรว์ D: และตั้งชื่อโฟลเดอร์สำหรับทดสอบการเขียนโปรแกรมนี้ว่า “Test” และตั้งชื่อโปรแกรมว่า “Hello”

จะเห็นหน้าตาเครื่องมือพัฒนาประมาณรูปด้านล่าง และตัวอีมูเลเตอร์ “Display” และ “Keyboard” ของ fx-9860G

Project แรกเริ่ม

เมื่อเราสร้าง Project ใหม่จะเห็นโครงร่างที่เครื่องมือเขียนมาให้ดังนี้

/*****************************************************************/
/*                                                               */
/*   CASIO fx-9860G SDK Library                                  */
/*                                                               */
/*   File name : Hello.c                                 */
/*                                                               */
/*   Copyright (c) 2006 CASIO COMPUTER CO., LTD.                 */
/*                                                               */
/*****************************************************************/
#include "fxlib.h"

//****************************************************************************
//  AddIn_main (Sample program main function)
//
//  param   :   isAppli   : 1 = This application is launched by MAIN MENU.
//                        : 0 = This application is launched by a strip in eACT application.
//
//              OptionNum : Strip number (0~3)
//                         (This parameter is only used when isAppli parameter is 0.)
//
//  retval  :   1 = No error / 0 = Error
//
//****************************************************************************
int AddIn_main(int isAppli, unsigned short OptionNum)
{
    unsigned int key;

    Bdisp_AllClr_DDVRAM();

    locate(1,4);
    Print((unsigned char*)"This application is");
    locate(1,5);
    Print((unsigned char*)" sample Add-In.");

    while(1){
        GetKey(&key);
    }

    return 1;
}

//****************************************************************************
//**************                                              ****************
//**************                 Notice!                      ****************
//**************                                              ****************
//**************  Please do not change the following source.  ****************
//**************                                              ****************
//****************************************************************************

#pragma section _BR_Size
unsigned long BR_Size;
#pragma section

#pragma section _TOP

//****************************************************************************
//  InitializeSystem
//
//  param   :   isAppli   : 1 = Application / 0 = eActivity
//              OptionNum : Option Number (only eActivity)
//
//  retval  :   1 = No error / 0 = Error
//
//****************************************************************************
int InitializeSystem(int isAppli, unsigned short OptionNum)
{
    return INIT_ADDIN_APPLICATION(isAppli, OptionNum);
}
#pragma section

จากโค้ดด้านบนจะเห็นว่าไม่เห็นจุดเริ่มต้นเข้าโปรแกรม (entry point) ฟังก์ชัน main() เช่นภาษาซีทั่วๆไป แต่จะมี AddIn_main() มาแทนให้รู้ว่าเป็นจุดเริ่มต้นของโปรแกรม AddIn สำหรับเครื่องคิดเลขนี้

มาดูโค้ดกันสักนิด ที่ #inlucde “fxlib.h” จะเป็น header ของ Casio สำหรับการแสดงผลเช่นฟังก์ชัน Print เพื่อแสดงบนจอภาพเครื่องคิดเลข Bdisp_AllClr_DDVRAM(); จะเป็นฟังก์ชันเคลียร์หน้าจอภาพให้ว่างเปล่า  locate(1,4); เลื่อนเคอเซอร์มาที่คอลัมน์ 1 และบรรทัดที่ 4 จากนั้นพิมพ์ด้วยฟังก์ชัน Print((unsigned char*)”This application is”); สุดท้ายปิดด้วย loop ไม่รู้จบ while (1) แล้วรอผู้ใช้กดคีย์ GetKey(&key); ตอนรันโปรแกรมถ้าผู้ใช้กดคีย์ “MENU” บนเครื่องคิดเลขก็จะเข้าสู่โหมด “MAIN MENU” แต่โปรแกรมก็ยังรันค้างอยู่ในสถานะเดิม

Compile & Build

มาลองคอมไพล์และบิวด์ดู ที่เมนู “Project” คลิก “Rebuild all” ถ้าไม่มีอะไรผิดพลาดที่กรอบ “Builds” จะแสดงผลดังนี้

Executing Hitachi SH C/C++ Compiler/Assembler phase

set SHC_INC=C:\fx-9860G-SDK\OS\SH\include
set PATH=C:\fx-9860G-SDK\OS\SH\bin
set SHC_LIB=C:\fx-9860G-SDK\OS\SH\bin
set SHC_TMP=D:\FX9860GIISD\Test\Debug
"C:\fx-9860G-SDK\OS\SH\bin\shc.exe" -subcommand=C:\Users\priabroy\AppData\Local\Temp\hmk9A12.tmp

Executing Hitachi OptLinker04 phase

"C:\fx-9860G-SDK\OS\SH\bin\Optlnk.exe" -subcommand=C:\Users\priabroy\AppData\Local\Temp\hmk9D5F.tmp

Optimizing Linkage Editor Completed

HMAKE MAKE UTILITY Ver. 1.1
Copyright (C) Hitachi Micro Systems Europe Ltd. 1998
Copyright (C) Hitachi Ltd. 1998 


	Make proccess completed

"D:\FX9860GIISD\Test\HELLO.G1A" was created.

Build has completed.

ถ้าสำเร็จจะเห็น HELLO.G1A ถูกสร้างขึ้นมา จากนั้นที่เมนู “Run” คลิกที่เมนูย่อย “Run” อีมูเลเตอร์จะเริ่มทำงาน

ที่ “Keyboard” กดปุ่มเลื่อนลูกศรลงที่ไอคอน “Debug” กดคีย์ “EXE” จะเห็นหน้าตาโปรแกรมดังนี้

ดาวน์โหลดซอร์สโค้ด (Source code) โปรแกรมแปลงพิกัดภูมิศาสตร์ (Geographic Calc)

เพื่อการลดระยะเวลาการเรียนรู้สำหรับคนที่เพิ่งจะมาศึกษาการใช้งานเครื่องมือ SDK ผมจะขอแสดงโครงการที่ผมทำไว้แล้ว ไปที่หน้าดาวน์โหลด (Download) มองหาซอร์สโค้ด (Source code) โปรแกรมสำหรับเครื่องคิดเลข Casio fx-9860G II SD จะได้ไฟล์ zip ชื่อ “UTM-Geo.zip” ทำการแตกไฟล์ zip จะเห็นไฟล์ดังนี้

ที่ผมวงสีแดงไว้คือโค้ดที่ได้จาก ไลบรารีที่ผมดาวน์โหลดมาใช้จาก githubพัฒนาโดย Howard Butler ส่วนที่ผมวงสีฟ้าไว้คือโค๊ดของผมเอง แตกไฟล์ zip ไปไว้ที่โฟลเดอร์สำหรับผมเองอยู่ที่ไดรว์ D:\FX9860GIISD

เปิดโครงการโปรแกรมแปลงพิกัดภูมิศาสตร์

ใช้เมนู “Project” > “Open…” เปิดไฟล์ชื่อ “UTMGeo.g1w”  ที่ช่อง panel  ด้านซ้ายสุดของเครื่องมือพัฒนาจะเห็นเป็นชื่อไฟล์ที่ผมได้ add มาไว้ จะสังเกตเห็นชื่อไฟล์บางตัวมี นามสกุล extension *.h, *.c ที่เป็นปกติของไฟล header และซอร์สโค้ดของภาษาซีและ *.hpp, *.cpp ของภาษาซีพลัสพลัส ซึ่งไฟล์ extension ที่เราตั้งชื่อไว้ถ้ามีโค้ดภาษาซีพลัสพลัส จะต้องใช้ extension เป็น hpp และ cpp

แก้ไขโครงการ

ใช้เมนู “Project” > “Edit…” เพื่อแก้ไขชื่อโครงการ หรือเพิ่มลดไฟล์ header และ source file

แก้ไขไอคอนสำหรับปรากฎที่หน้า “MAIN MENU” ของเครื่องคิดเลข ไอคอนขนาดกว้าง 39 pixel และสูงขนาด 19 pixel ฟอร์แม็ตเป็น bitmap แบบขาวดำ 1 bit ผมออกแบบในโปรแกรม Paint.net จัดเก็บเป็นไฟล์ bitmap (bmp) แต่ไม่สามารถจัดเก็บเป็นไฟล์ขาวดำแบบ 1 bit ได้ต้องไปเปิดต่อใน Gimp แล้วเลือกเมนู Image>Mode>Indexed ตรง Color map เลือก Use black and white (1 bit) palette จากนั้น save จึงจะได้ไฟล์รูปที่มีฟอร์แม็ต bitmap แบบ 1 bit ได้ ทั้งสองโปรแกรมนี้ฟรี

ตรง Edit icon เวลาจะคลิกต้องระวังเพราะถ้าเราออกแบบไอคอนมาแล้ว จะโดนทับด้วยรูปไอคอนดีฟอลท์ทันที ผมไม่ใช้เลยเพราะพลาดหลายทีแล้ว

การจัดการโค้ด C++

ถ้ามีโค้ดซีพลัสพลัสมาผสมด้วย ที่ด้านบนสุดไฟล์จะต้องกำหนดดังนี้

#ifdef __cplusplus
  extern "C" {
#endif

ด้านล่างสุดปิดท้ายด้วย

  #ifdef __cplusplus
}
  #endif

โค้ดจัดการการป้อนข้อมูล

การป้อนข้อมูลของงานสำรวจในเครื่องคิดเลขส่วนใหญ่จะเป็นเลขทศนิยมเช่นระยะทางหรือค่าพิกัด แต่ถ้าเป็นมุมทีแยกองศา ลิปดาและฟิลิปดา ปกติในเครื่องคิดเลขเช่น  fx-4500, fx-5800P จะมีคีย์ให้กดสะดวก แต่เครื่อง  fx-9860G II SD กลับเอาไปไว้ลึกมากต้องกดหลายครั้งจากคีย์ “OPTN” แต่ที่ผมช็อคคือใน SDK กลับไม่มีฟังก์ชันให้เรียกใช้งานได้เลย จะต้องเขียนฟังก์ชันขึ้นมาเองทั้งหมด ผมอาศัยไปอ่านตามฟอรั่มที่มีคนแฮ็คไว้ พบว่าสามารถเรียกใช้ฟังก์ชัน EditExpress ที่ทาง Casio ไม่ได้เปิดเผยเอกสารไว้ (สังเกตว่าใช้เป็น function pointer อ่านรายละเอียดการใช้งานได้ที่ลิ๊งค์นี้)

#define SCA 0xD201D002
#define SCB 0x422B0009
#define SCE 0x80010070

const unsigned int sc08DB[] = {SCA, SCB, SCE, 0x08DB};
typedef int(*sc_EE)(int, short, int, char*, char*, short, char*, int );
#define EditExpression (*(sc_EE)sc08DB)

วิธีใช้งานก็ประมาณนี้

#define MAXEDITBUFFER 21
//
void AddIn_main(int isAppli, unsigned short OptionNum){
int key;
char vBCD[24];
unsigned char sBCD[MAXEDITBUFFER];

memset( vBCD, 0, sizeof(vBCD));
memset( sBCD, 0, sizeof(sBCD));
key = EditExpression(0, KEY_CTRL_RIGHT, 1, vBCD, (char*)sBCD, MAXEDITBUFFER - 1, "Input:", 0x04);

locate(1, 4);
PrintLine(sBCD, 21);

GetKey(&key);
return 1;
}

ข้อมูลที่ผู้ใช้ป้อนจะกลับมาที่ข้อมูลสตริง sBCD ผมพบว่าเรียกใช้ EditExpression มี  mode ให้ป้อนเลือกว่าจะป้อนข้อมูลเป็น double ไหม ซึ่งตอนป้อนข้อมูลจะรับแค่ตัวเลข แต่ในโหมดตามตัวอย่างด้านบน (พารามิเตอร์สุดท้าย 0x04) นั้นจะรับค่าทั่วๆไปทั้งตัวอักษรผสมกับตัวเลข แต่ผมสังเกตว่าเวลาเราป้อนข้อมูลแล้วกดคีย์ “EXE” จะมีการประมวลผล expression ด้วย ปัญหาที่ผมพบคือถ้าป้อนสตริงค่าพิกัดแบบ MGRS เช่น “18SVK8588822509” กลับ error ผมเลยต้องเขียนฟังก์ชันขึ้นมาเองคือ inputMGRSString() โดยเฉพาะ

การป้อนข้อมูลมุม

เนื่องจากเครื่องมือพัฒนา SDK ไม่สนับสนุนการป้อนมุมแบบใช้งานปกติ ผมเลยกำหนดว่ามุมองศา ลิปดาและฟิลิปดาให้คั่นด้วยเครื่องหมายลบ (-) เช่นแลตติจูด (Latitude) จะใช้ตัวอักษร “N” หรือ “S” มากำกับว่าอยู่ในซึกโลกเหนือหรือซีกโลกใต้ คั่นด้วยเส้นศูนย์สูตร ส่วนมุมลองจิจูด (Longitude) แบ่งเป็นซีกโลกตะวันออก (“E” ปิดท้าย) และซีกโลกตะวันตก (“W” ปิดท้าย

ตัวอย่าง แลตติจูด 45-5-32.525N ลองจิจูด 98-45-38.587W

โค้ดอ่านเขียนข้อมูลเข้าตัวแปร  Alpha

ส่วนใหญ่เวลาเราป้อนข้อมูลเข้าไปในการคำนวณงาน เราต้องการให้โปรแกรมจดจำค่านั้นไว้ ดังนั้นโปรแกรมต้องมีการจัดเก็บเอาไว้แล้วเรียกมาใช้ทีหลังเมื่อเรียกโปรแกรมมาใช้งานอีกที เพราะไม่ต้องป้อนบ่อย ถ้ายังใช้ค่าเดิมแค่กดคีย์ “EXE” ผ่านได้เลย และก็เหมือนเดิมครับ SDK ไม่ได้เปิดเผยฟังก์ชันนี้ไว้ทั้งที่สำคัญมาก ผมไปค้นหาตามฟอรั่มพบว่าการจัดเก็บข้อมูลเข้าตัวแปรตัวอักษร A-Z ใช้โครงสร้างข้อมูลเฉพาะ อ่านได้ตามลิ๊งค์นี้ การแปลงข้อมูลสามารถใช้ class TBCD (โค๊ดเดิมมีบั๊กผมแก้ไปนิดหน่อย)

#define SCA 0xD201D002
#define SCB 0x422B0009
#define SCE 0x80010070

const unsigned int sc04E0[] = {SCA, SCB, SCE, 0x04E0};
const unsigned int sc04DF[] = {SCA, SCB, SCE, 0x04DF};

typedef void(*sc_agd)(char, TBCDvalue*);
typedef void(*sc_asd)(char, TBCDvalue*);

#define Alpha_GetData (*(sc_agd)sc04DF)
#define Alpha_SetData (*(sc_asd)sc04E0)

void GetAlphaDoubleData(char alpha, double *dval);
void SetAlphaDoubleData(char alpha, double val);

typedef struct{
  unsigned char hnibble:4;
  unsigned char lnibble:4;
} TBCDbyte;

typedef struct{
  unsigned short exponent:12;
  unsigned short mantissa0:4;
  TBCDbyte mantissa[7];
  char flags;
  short info;
} TBCDvalue;

typedef struct{
  int exponent;
  int sign;
  int unknown; 
  char mantissa[15];
} TBCDInternal;

//Implement of class TBCD please see utilities.cpp
class TBCD{
  public:
        TBCDvalue*PValue();
        int Get( TBCDvalue&value );
        int Set( TBCDvalue&value );
        int Set( double&value );
        int Get( double&value );
        int SetError( int error );
        int GetError();
        void Swap();
  protected:
  private:
        TBCDvalue FValue[2];
};

void GetAlphaDoubleData(char alpha, double *dval){
  TBCD *bcd;
  TBCDvalue *bval;
  int i, ii;

  bval = (TBCDvalue *)malloc(sizeof(TBCDvalue));
  Alpha_GetData(alpha, bval);
  bcd = new TBCD;
  i = bcd->Set(*bval);
  ii = bcd->Get(*dval);
  delete bcd;
  free(bval);
}

void SetAlphaDoubleData(char alpha, double dval){
  TBCD *bcd;
  TBCDvalue bval;
  int i, ii;

  bcd = new TBCD;
  i = bcd->Set(dval);
  ii = bcd->Get(bval);
  Alpha_SetData(alpha, &bval);
  delete bcd;
}
วิธีการใช้งานก็ง่ายๆ

      SetAlphaDoubleData('B', 123.456); //เอาค่า 123.456 เข้าเก็บที่ตัวแปรอักษร "B"
      GetAlphaDoubleData('B', &x); //ดึงค่าที่เก็บในตัวอักษร "B" ออกมาเข้าตัวแปร x
      Locate(1, 4);
      sprintf((char*) str, (char*) "Value = %.3lf", x); 
      Print((unsigned char*)str); //Value = 123.456 

วิธีการใช้งานก็ง่ายๆ

      SetAlphaDoubleData('B', 123.456); //เอาค่า 123.456 เข้าเก็บที่ตัวแปรอักษร "B"
      GetAlphaDoubleData('B', &x); //ดึงค่าที่เก็บในตัวอักษร "B" ออกมาเข้าตัวแปร x
      Locate(1, 4);
      sprintf((char*) str, (char*) "Value = %.3lf", x); 
      Print((unsigned char*)str); //Value = 123.456  

โค้ดสำหรับเมนูหลัก

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

โค้ดก็ง่ายๆดังนี้

      memset(s, '-', 21);
      Bdisp_AllClr_DDVRAM(); 
      locate(0, 1);
      Print((unsigned char *)"Geographic Calc");
      locate(0, 2);
      PrintLine((unsigned char*)s, 21);
      locate(0, 3);
      PrintLine("[1]:UTM to Geo", 21);
      locate(0, 4);
      PrintLine("[2]:Geo to UTM", 21);
      locate(1, 5);
      PrintLine("[3]:MGRS to Geo", 21);
      locate(0, 6);
      PrintLine("[4]:Geo to MGRS", 21);
      locate(1, 8);
      PrintLine("Select 1,2,3 or 4", 21);
      while (!((key1 >= 0x31) && (key1 <= 0x34)){ 
        GetKey(&key1);

จะมีลูป while ดักการกดคีย์อีก loop ถ้าพบว่ากดคีย์เลข “1” (character code = 0x31) ถึงเลข “4” (Character code = 0x34) เงื่อนไขจริงจะออกจาก loop เข้าเงื่อนไข if (มีหลายชั้นใช้ case แทนได้)

โค้ดงานคำนวณหลัก

ถ้าผู้ใช้กดคีย์ “1” จะเป็นการคำนวณค่าพิกัดจาก UTM ไปยัง Geographic

    if (key1 == 0x31) { //UTM to Geo
      Bdisp_AllClr_DDVRAM(); 
      locate(0, 1);
      Print((unsigned char *)"UTM to Geo");
      locate(0, 2);
      PrintLine((unsigned char*)s, 21);

      GetAlphaDoubleData('A', A);//ดึงค่าข้อมูลเดิมจากหน่วยความจำตัวอักษร "A" 
      sprintf((char*)sBCD, (char*) "%.3lf", *A); //เตรียมรูปแบบข้อมูลทศนิยมสามตำแหน่งไว้ใน sBCD
      //เรียกฟังก์ชันป้อนข้อมูล โดยส่ง sBCD ไปให้
      key2 = EditExpression(0, KEY_CTRL_RIGHT, 3, vBCD, (char*)sBCD, MAXEDITBUFFER - 1, "N? : ", 0x04);
      y = atof((char*)sBCD); //แปลงข้อมูลที่ป้อนมาเป็นตัวเลขทศนิยม
      SetAlphaDoubleData('A', y);//เก็บค่าที่ป้อนไว้ในหน่วยความจำตัวอักษร "A"

      GetAlphaDoubleData('B', B);
      sprintf((char*)sBCD, (char*) "%.3lf", *B); 
      key2 = EditExpression(0, KEY_CTRL_RIGHT, 4, vBCD, (char*)sBCD, MAXEDITBUFFER - 1, "E? : ", 0x04);
      x = atof((char*)sBCD);
      SetAlphaDoubleData('B', x);
     
      GetAlphaDoubleData('Z', Z);
      sprintf((char*)sBCD, (char*) "%.3lf", *Z); 
      key2 = EditExpression(0, KEY_CTRL_RIGHT, 5, vBCD, (char*)sBCD, MAXEDITBUFFER - 1, "UTM Zone No.? : ", 0x04);
      if (strchr((char*)sBCD, 0x87) != NULL) {
        removechar((char*)sBCD, 0x87);
        hemi = 'S';
      } else if(strchr((char*)sBCD, 0x2D) != NULL) {
        removechar((char*)sBCD, 0x2D);
        hemi = 'S';
      } else
        hemi = 'N';

      zn = atoi((char*)sBCD);
      SetAlphaDoubleData('Z', zn);
      //เรียกไลบรารีแปลงพิกัดที่ประกาศใน utm.c
      err = Convert_UTM_To_Geodetic(zn, hemi, x, y, &lat, &lng);
      if (!err){
	      Lat = lat * RAD2DEG; 
	      Lng = lng * RAD2DEG;
              //แยกทศนิยมจัดรูปแบบที่ที่องศา ลิปดาและฟิลิปดาคั่นด้วยเครื่องหมายลบ
	      slat = degreetodms(fabs(Lat), NUMDECIMAL, 0x2D);
	      if(Lat >= 0)
	        sprintf(str, "Lat= %s N", slat);
	      else
	        sprintf(str, "Lat= %s S", slat);
	      locate(1, 6);
	      Print((unsigned char*)str); 
	      locate(1, 7);
	      slong = degreetodms(fabs(Lng), NUMDECIMAL, 0x2D);
	      if (Lng >= 0)
	        sprintf(str, "Lon= %s E", slong);
	      else
	        sprintf(str, "Lon= %s W", slong);
	      Print((unsigned char*)str);      
	      free(slat);
	      free(slong);
       }

ลองดูโค้ดแปลงพิกัดจากค่าพิกัดฉากยูทีเอ็มไปค่าพิกัดภูมิศาสตร์ โดยใช้ไลบรารี

err = Convert_UTM_To_Geodetic(zn, hemi, x, y, &lat, &lng);

zn คือโซนยูทีเอ็ม hemi คือซีกโลก x และ y คือค่าพิกัดฉากยูทีเอ็ม ส่วนค่าที่คำนวณแล้วจะส่งกลับมาที่ตัวแปร lat, lng ส่วนการแปลงพิกัดอย่างอื่นก็ลองดูได้ตามโค้ด

คอมไพล์และบิวด์ (compile & build)

ทดสอบโดยเมนู “Project” > “Rebuild all” ถ้าเมนูนี้ไม่ขึ้นให้คลิก “Project” > “Reload” ก่อน จากทำการรันโปรแกรมโดย “Run” > “Run” จะเห็นอีมูเลเตอร์ “Display” เลื่อนกดคีย์ ไปที่ไอคอนโปรแกรม จากอีมูเลเตอร์ “Keyboard” แล้วกดคีย์ “EXE”

จะเห็นโปรแกรม “Geographic Cacl” ขึ้นเมนูมา

การใช้งานโปรแกรมนี้พร้อมตัวอย่างไปดูที่ลิ๊งค์นี้ได้

ซอร์สโค้ดหลัก

ท้ายสุดผมเอาซอร์สโค้ดของไฟล์ “main.cpp” มาลงให้ดูเต็มๆ จะเห็นว่าไม่มีอะไรสลับซับซ้อนมาก ง่ายๆครับ

#ifdef __cplusplus
  extern "C" {
#endif

#include "fxlib.h"
#include "string.h"
#include "stdlib.h"
#include "stdio.h"
#include "math.h"
#include "mgrs.h"
#include "tranmerc.h"
#include "utm.h"
#include "utilities.hpp"


#define MAXEDITBUFFER 21
#define NUMDECIMAL 5


void removechar(char* s, const char toremove)
{
  while(s=strchr(s, toremove))   
    memmove(s, s+1, 1+strlen(s+1));
}

int InputMGRSString(int x, int y, unsigned char*prompt, unsigned char*buffer, int bufferSize ){
   unsigned int key, edit_key, return_key = 0;
   int pos, len, len2;
   Cursor_SetFlashMode(1);   // set cursor visibility on
   pos = strlen((char*)buffer);// initially set the cursor to the end of the string

   locate(x, y);
   Print(prompt);
   len2 = strlen((char*)prompt);
   while (!return_key){
      locate(x + len2, y);
      PrintLine(buffer, 22-x);
      locate (x + pos + len2, y);
      GetKey( &key );
      edit_key = 0;
      if ((key >= 0x30 ) && (key <= 0x39)){
         edit_key = key;
      } else if ((key >= 0x41) && (key <= 0x5A)){
         edit_key = key;
      } else{
         switch (key){
            case KEY_CTRL_DEL :
               if ( pos > 0 ) pos--;
               len = strlen( (char*)buffer );   // get the current length of the string
               memmove( buffer+pos, buffer+pos+1, len-pos);   // shift the memory: XXYDYYY -> XXXYYY
               break;
            case KEY_CTRL_RIGHT :
               len = strlen( (char*)buffer );   // get the current length of the string
               if ( pos < len ) pos++;
               break;
            case KEY_CTRL_LEFT :
               if ( pos > 0 ) pos--;
               break;
            case KEY_CTRL_UP :
            case KEY_CTRL_DOWN :
               return_key = key;
               break;
            case KEY_CTRL_EXE :
            case KEY_CTRL_EXIT :
               return_key = key;
               break;
            default :
               break;
         };
      }
      if ( edit_key ){
         if ( pos < bufferSize-1 ){
            buffer[ pos ] = edit_key;
            pos++;
         }
      }
   }
   Cursor_SetFlashMode( 0 );   // set cursor visibility off
   return ( return_key );
}

int AddIn_main(int isAppli, unsigned short OptionNum)
{
  unsigned char buffer[21];
  char str[21], s[21];
  int editresult;
  unsigned int key1, key2;
  double x, y, lat, lng, Lat, Lng;
  char *slat, *slong, *sangle;
  long zn; 
  char hemi;//North hemi is 'N', South hemi is 'S'
  char vBCD[24];
  unsigned char sBCD[MAXEDITBUFFER] = ""; 
  int err;
  unsigned char mgrs[15];
  long precision;
  double *A, *B, *C, *D, *E, *F, *H, *I, *Z;

  A = (double *)malloc(sizeof(double));
  B = (double *)malloc(sizeof(double));
  C = (double *)malloc(sizeof(double));
  D = (double *)malloc(sizeof(double));
  E = (double *)malloc(sizeof(double));
  F = (double *)malloc(sizeof(double));
  H = (double *)malloc(sizeof(double));
  I = (double *)malloc(sizeof(double));
  Z = (double *)malloc(sizeof(double));
  memset(s, '-', 21);
  while(1){
	Bdisp_AllClr_DDVRAM(); 
	locate(0, 1);
	Print((unsigned char *)"Geographic Calc");
	locate(0, 2);
	PrintLine((unsigned char*)s, 21);
	locate(0, 3);
	PrintLine("[1]:UTM to Geo", 21);
	locate(0, 4);
	PrintLine("[2]:Geo to UTM", 21);
	locate(1, 5);
	PrintLine("[3]:MGRS to Geo", 21);
	locate(0, 6);
	PrintLine("[4]:Geo to MGRS", 21);
	locate(1, 8);
	PrintLine("Select 1,2,3 or 4", 21);
    while (!(key1 >= 0x31) && (key1 <= 0x34)){ 
      GetKey(&key1);
    }
    if (key1 == 0x31) { //UTM to Geo
      Bdisp_AllClr_DDVRAM(); 
      locate(0, 1);
      Print((unsigned char *)"UTM to Geo");
      locate(0, 2);
      PrintLine((unsigned char*)s, 21);

      GetAlphaDoubleData('A', A);
      sprintf((char*)sBCD, (char*) "%.3lf", *A); 
      key2 = EditExpression(0, KEY_CTRL_RIGHT, 3, vBCD, (char*)sBCD, MAXEDITBUFFER - 1, "N? : ", 0x04);
      y = atof((char*)sBCD);
      SetAlphaDoubleData('A', y);

      GetAlphaDoubleData('B', B);
      sprintf((char*)sBCD, (char*) "%.3lf", *B); 
      key2 = EditExpression(0, KEY_CTRL_RIGHT, 4, vBCD, (char*)sBCD, MAXEDITBUFFER - 1, "E? : ", 0x04);
      x = atof((char*)sBCD);
      SetAlphaDoubleData('B', x);
     
      GetAlphaDoubleData('Z', Z);
      sprintf((char*)sBCD, (char*) "%.3lf", *Z); 
      key2 = EditExpression(0, KEY_CTRL_RIGHT, 5, vBCD, (char*)sBCD, MAXEDITBUFFER - 1, "UTM Zone No.? : ", 0x04);
      if (strchr((char*)sBCD, 0x87) != NULL) {
        removechar((char*)sBCD, 0x87);
        hemi = 'S';
      } else if(strchr((char*)sBCD, 0x2D) != NULL) {
        removechar((char*)sBCD, 0x2D);
        hemi = 'S';
      } else
        hemi = 'N';

      zn = atoi((char*)sBCD);
      SetAlphaDoubleData('Z', zn);
      //Call function declared in utm.c
      err = Convert_UTM_To_Geodetic(zn, hemi, x, y, &lat, &lng);
      if (!err){
	      Lat = lat * RAD2DEG; 
	      Lng = lng * RAD2DEG;
	      slat = degreetodms(fabs(Lat), NUMDECIMAL, 0x2D);
	      if(Lat >= 0)
	        sprintf(str, "Lat= %s N", slat);
	      else
	        sprintf(str, "Lat= %s S", slat);
	      locate(1, 6);
	      Print((unsigned char*)str); 
	      locate(1, 7);
	      slong = degreetodms(fabs(Lng), NUMDECIMAL, 0x2D);
	      if (Lng >= 0)
	        sprintf(str, "Lon= %s E", slong);
	      else
	        sprintf(str, "Lon= %s W", slong);
	      Print((unsigned char*)str);      
	      free(slat);
	      free(slong);
       }
    } else if (key1 == 0x32) { //Geographic to UTM
      Bdisp_AllClr_DDVRAM(); 
      locate(0, 1);
      Print((unsigned char *)"Geo To UTM");
      locate(0, 2);
      PrintLine((unsigned char*)s, 21);

      GetAlphaDoubleData('H', H);
      sangle = degreetodms(*H, NUMDECIMAL, 0x99); 
      if(*H >= 0)
        sprintf(str, "%sN", sangle);
      else
        sprintf(str, "%sS", sangle);
      key2 = EditExpression(0, KEY_CTRL_RIGHT, 3, vBCD, (char*)str, MAXEDITBUFFER - 1, "Lat?: ", 0x04);
      lat = parsedms((char*)str);
      SetAlphaDoubleData('H', lat);
      free(sangle);

      GetAlphaDoubleData('I', I);
      sangle = degreetodms(*I, NUMDECIMAL, 0x99); 
      if(*I >= 0)
        sprintf(str, "%sE", sangle);
      else
        sprintf(str, "%sW", sangle);
      key2 = EditExpression(0, KEY_CTRL_RIGHT, 4, vBCD, (char*)str, MAXEDITBUFFER - 1, "Lon?: ", 0x04);
      lng = parsedms((char*)str);
      SetAlphaDoubleData('I', lng);
      free(sangle);

      Lat = lat * DEG2RAD;
      slat = degreetodms(fabs(Lat), NUMDECIMAL, 0x2D);
      Lng = lng * DEG2RAD;
      slong = degreetodms(fabs(Lng), NUMDECIMAL, 0x2D);
      //Call function declared in utm.c
      err = Convert_Geodetic_To_UTM(Lat, Lng, &zn, &hemi, &x, &y);
      if (!err) {
	      sprintf(str, "North= %11.3lf", y);
	      locate(1, 5);
	      Print((unsigned char*)str);      
	      sprintf(str, "East= %10.3lf", x);
	      locate(1, 6);
	      Print((unsigned char*)str);      
	      sprintf(str, "UTM Zone No= %d%c", zn, hemi);
	      locate(1, 7);
	      Print((unsigned char*)str); 
      }
      free(slat);
      free(slong);      
    } else if (key1 == 0x33) { //MGRS to Geo
      Bdisp_AllClr_DDVRAM(); 
      locate(0, 1);
      Print((unsigned char *)"MGRS to Geo");
      locate(0, 2);
      PrintLine((unsigned char*)s, 21);

      memset(mgrs, 0, 16);
      editresult = InputMGRSString( 1, 3, "MGRS?", buffer, sizeof(buffer) );
      if (editresult) {
        memcpy(mgrs, buffer, 15);
        //Call function that declared in mgrs.c
        err = Convert_MGRS_To_Geodetic((char*)mgrs, &lat, &lng);
        if (!err){
          Lat = lat * 180.0 / PI;
          slat = degreetodms(fabs(Lat), NUMDECIMAL, 0x2D);
          Lng = lng * 180.0 / PI;
          slong = degreetodms(fabs(Lng), NUMDECIMAL, 0x2D);
          if(Lat >= 0)
            sprintf(str, "Lat= %s N", slat);
          else
            sprintf(str, "Lat= %s S", slat);
          locate(1, 4);
          Print((unsigned char*)str);      
          if(Lng >= 0)
            sprintf(str, "Lon= %s E", slong);
          else
            sprintf(str, "Lon= %s W", slong);
          locate(1, 5);
          Print((unsigned char*)str);        
          free(slat);
          free(slong);
        } 
      }
    } else if (key1 == 0x34) { //Geographic to MGRS
      Bdisp_AllClr_DDVRAM(); 
      locate(0, 1);
      Print((unsigned char *)"MGRS to Geo");
      locate(0, 2);
      PrintLine((unsigned char *)s, 21);

      memset(sBCD, 0, MAXEDITBUFFER);
      key2 = EditExpression(0, KEY_CTRL_RIGHT, 3, vBCD, (char*)sBCD, MAXEDITBUFFER - 1, "Lat?:", 0x04);
      lat = parsedms((char*)sBCD);
      lat = lat * DEG2RAD;
      memset(sBCD, 0, MAXEDITBUFFER);
      key2 = EditExpression(0, KEY_CTRL_RIGHT, 4, vBCD, (char*)sBCD, MAXEDITBUFFER - 1, "Lon?:", 0x04);
      lng = parsedms((char*)sBCD);
      lng = lng * DEG2RAD;

      //Call function that declared in mgrs.c
      err = Convert_Geodetic_To_MGRS(lat, lng, 5, (char*)mgrs);
      if (!err){
        sprintf(str, "MGRS= %s", mgrs);
        locate(1, 5);
        Print((unsigned char*)str);        
      }     
    }
    GetKey(&key1);
  } //while (1)
  free(A);
  free(B);
  free(C);
  free(D);
  free(E);
  free(F);
  free(H);
  free(I);
  free(Z);
  return 1;
}

#pragma section _BR_Size
unsigned long BR_Size;
#pragma section

#pragma section _TOP

int InitializeSystem(int isAppli, unsigned short OptionNum)
{
    return INIT_ADDIN_APPLICATION(isAppli, OptionNum);
}

#pragma section

#ifdef __cplusplus
}
#endif

การนำโปรแกรมไปติดตั้งบนเครื่องคิดเลข

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

สรุป

การพัฒนาโปรแกรมสำหรับเครื่องคิดเลขไม่มีอะไรยากเย็นมากนัก มีความรู้ภาษาซีขั้นพื้นฐานแบบผมก็ทำได้ เพียงแต่ถ้าต้องการคำนวณอะไรซับซ้อนอาจจะต้องหาไลบรารีที่ท่านอื่นได้พัฒนาเขียนไว้ แต่เครื่องมือพัฒนานี้มีข้อจำกัดอยู่บ้างเช่นคอมไพเลอร์นี้สนับสนุน c standard library ได้ไม่ครบทุกอย่าง ตัวอย่างเช่นฟังก์ชั่น strtok() ที่ตัดสตริงออกตามตัวคั่นก็ไม่สนับสนุน ดังนั้นการเลือกไลบรารีที่คนอื่นได้ทำไว้ก็ต้องพิจารณาส่วนนี้ด้วย ส่วนการรับข้อมูลจากผู้ใช้ สำหรับผมแล้วที่มีอยู่ตอนนี้เกือบจะเพียงพอ เพราะงานสำรวจก็จะมีแค่ป้อนมุม ระยะทาง เป็นหลัก ติดตามกันต่อไปครับ

คอมไพล์ Python Script เป็นไฟล์ Executable ด้วย PyInstaller

PyInstaller

คือเครื่องมือที่ช่วยการแปลงโปรแกรมที่เขียนด้วยไพทอนเป็น execute binary file  ที่สามารถนำไปรันได้โดยที่เครื่องคอมพิวเตอร์ปลายทางไม่ต้องติดตั้งไพทอน สำหรับ PyInstaller เป็น cross-platform สามารถใช้งานได้บนวินโดส์ แมค และลีนุกซ์ สนับสนุนไพทอนรุ่น 2.7 และ ไพทอน รุ่น 3.3 ถึง 3.6 จุดมุ่งหมายของ PyInstaller คือต้องการช่วยผู้ใช้ในการแปลงโปรแกรมไพทอน ที่ใช้โมดูลไลบรารีภายนอกเช่น Matplotlib, DJango, wxPython, PyQt เป็นต้น ให้สามารถทำได้ง่ายสะดวก

ติดตั้ง PyInstaller

ติดตั้งง่ายๆด้วยคำสั่ง pip ใน command prompt

pip install pyinstaller

ใช้งาน PyInstaller

การใช้งานสามารถใช้งานผ่าน command line ได้ แต่สำหรับโปรแกรมที่เรียกใช้โมดูลไลบรารีข้างนอกและต้องขนข้อมูล (data) ที่โมดูลไลบรารีนั้นๆต้องการใช้  ผมแนะนำให้ใช้ไฟล์สคริปท์ (Spec file) มาช่วยจะดีกว่า ปรับแต่งได้มากกว่า ตัว spec file จริงๆก็คือไฟล์สคริปท์ของไพทอนนั่นเอง กรณีที่ต้องใช้ Spec file อีกกรณีหนึ่งคือต้องการขนรันไทม์ไลบรารีเช่น .dll หรือ .so ไปแบบแมนวล กรณีที่ผมเจอคือผมใช้ PySide2 ที่รุ่นทางการจริงๆยังไม่ออกมา แต่ hook file ก็มีมาให้แล้วพร้อมกับ PyInstaller รุ่นใหม่ 3.3 แต่ผมใช้งานแล้วยังไม่สำเร็จ ดังนั้นจึงต้องใช้ Spec file นี้เป็นตัวช่วยในการขนรันไทม์ไลบรารีไป ส่วนเรื่อง hook file คืออะไรค่อยว่าอีกที

กรณีศึกษาด้วย Surveyor Pocket Tools บนวินโดส์

โปรแกรม Surveyor Pocket Tools พัฒนาด้วยไพทอน ปัจจุบันใช้ไพทอน รุ่น 3.6 ใช้โมดูลไลบรารีข้างนอกคือ openpyxl, pyproj, geographiclib, gmplot, simplekml, pyshp และที่ขาดไม่ได้คือ PySide2 ซึ่งสำหรับ openpyxl และ pyproj จะมีการขนข้อมูลไปด้วย ส่วน PySide2 ผมจะขนไฟล์ dll  ที่ต้องการด้วยมือล้วนๆ

Spec file ของ Surveyor Pocket Tools

มาดูไฟล์สคริปนี้ ผมตั้งชื่อว่า “setup.spec”

# -*- mode: python -*-
import sys
import PySide2
import os
block_cipher = None

dirname = os.path.dirname(PySide2.__file__)
print("dirname=", dirname)

pyside2_resources_path = os.path.join(dirname, 'resources', '')
print("resources path=", pyside2_resources_path)
pyside2_resources = [(pyside2_resources_path + '*', '.'),
                     (dirname + '/qtwebengineprocess.exe', '.')]

added_files = [
         ( 'markers/*', 'markers/' ),
         ( 'geoidgrids/tgm2017.gtx', 'pyproj/data/' ),
         ( 'geoidgrids/egm96_15.gtx', 'pyproj/data/' ),
         ( 'geoidgrids/egm08_25.gtx', 'pyproj/data/' ),		 
		 ( 'database/*', 'database/'),
		 ( 'example data/*', 'example data/'),
		 ( 'qt.conf', '.'),
         ( 'settings.xml', '.'),
		 ( 'crs.xml', '.')]

a = Analysis(['main.py'],
             pathex=['D:\\sourcecodes\\python\\surveyor pocket tools'],
             binaries=None,
             datas=added_files + pyside2_resources,
             hiddenimports=[],
             hookspath=[],
             runtime_hooks=[],
             excludes=[],
             win_no_prefer_redirects=False,
             win_private_assemblies=False,
             cipher=block_cipher)
pyz = PYZ(a.pure, a.zipped_data,
             cipher=block_cipher)
exe = EXE(pyz,
          a.scripts,
          exclude_binaries=True,
          name='surveyor pocket tools',
		  icon='Land Surveying-64.ico',
          debug=False,
          strip=False,
          upx=False,
          console=False )
coll = COLLECT(exe,
               a.binaries,
               a.zipfiles,
               a.datas,
               strip=False,
               upx=True,
               name='setup')

ลองมาดูโค้ดกัน เริ่มจาก import PySide2 เข้ามาเพื่อจะตรวจสอบว่า PySide2 ที่เราใช้งานเป็น 32 บิตหรือ 64 บิต เพื่อจะได้ขนไฟล์ qtwebengineprocess.exe ที่ผมใช้แสดงผล Google Maps ตอนปักหมุด จากนั้นเก็บไดเรคทอรีของ PySide2 ไปเก็บไว้ในตัวแปร pyside2_resources_path ส่วนที่อยู่ของไฟล์เก็บที่ pyside2_resources

# -*- mode: python -*-
import sys
import PySide2
import os
block_cipher = None

dirname = os.path.dirname(PySide2.__file__)
print("dirname=", dirname)

pyside2_resources_path = os.path.join(dirname, 'resources', '')
print("resources path=", pyside2_resources_path)
pyside2_resources = [(pyside2_resources_path + '*', '.'),
                     (dirname + '/qtwebengineprocess.exe', '.')]

ต่อไปจะขนไฟล์ที่โปรแกรม Surveyor Pocket Tools ต้องการใช้ ให้ใส่ไว้ที่ตัวแปร added_files โครงสร้างเป็น tuple เหมือนกัน และขนไฟล์ชื่อ qt.conf ที่ PySide2 ต้องการไปด้วย

added_files = [
         ( 'markers/*', 'markers/' ),
         ( 'geoidgrids/tgm2017.gtx', 'pyproj/data/' ),
         ( 'geoidgrids/egm96_15.gtx', 'pyproj/data/' ),
         ( 'geoidgrids/egm08_25.gtx', 'pyproj/data/' ),		 
		 ( 'database/*', 'database/'),
		 ( 'example data/*', 'example data/'),
		 ( 'qt.conf', '.'),
         ( 'settings.xml', '.'),
		 ( 'crs.xml', '.')]

มาดูไดเรคทอรีที่โปรแกรมต้องการดังนี้

ต่อไปมาดูโค้ดส่วนที่สำคัญมาก ‘main.py’ คือไฟล์สคริปท์หลักของโปรแกรม Surveyor Pocket Tools ต่อไปคือ pathex เป็นไดเรคทอรีของไฟล์ไพทอนสคริปท์ และ datas ที่ผมจัดการรวม added_files และ pyside2_resources เข้าด้วยกัน สุดท้าย hookspath คือไดเรคทอรีที่เก็บไฟล์ hook ไว้ สำหรับไฟล์ hook นี้ PyInstaller จะอ่านสคริปท์นี้ทีละไฟล์มาตัดสินใจว่าจะขนข้อมูลไดเรคทอรีไหนไป ผมเลือกใช้ดีฟอลท์ครับคือปล่อยว่าง

a = Analysis(['main.py'],
             pathex=['D:\\sourcecodes\\python\\surveyor pocket tools'],
             binaries=None,
             datas=added_files + pyside2_resources,
             hiddenimports=[],
             hookspath=[],
             runtime_hooks=[],
             excludes=[],
             win_no_prefer_redirects=False,
             win_private_assemblies=False,
             cipher=block_cipher)

ใช้ PyInstaller คอมไพล์ไฟล์ setup.spec

ผมใช้ Minoconda เมื่อจะคอมไพล์ก็เรียก command prompt มาดังนี้ ใช้คำสั่ง cd เข้ามาที่พาทของสคริปท์ของไพทอน ใช้คำสั่ง dir ดูไฟล์ setup.spec

ต่อไปทำการคอมไพล์ ด้วยคำสั่ง

pyinstaller setup.spec

ผลลัพธ์ของ PyInstaller

ไฟล์ของไลบรารี Proj4 (pyproj) ที่ต้องใช้

เมื่อคอมไพล์เสร็จแล้ว ไม่มี error จะได้ไดเรคทอรีมาสองคือ “build” และ “dist” เมื่อเข้าไปดูใน “dist” จะเห็นไดเรคทอรีย่อยช “setup” ชื่อไดเรคทอรีนี้ PyInstaller จะสร้างตามชื่อหน้าของไฟล์ setup.spec เมื่อเข้าดูที่ไดเรคทอรี “setup” จะเห็นไฟล์ต่างๆที่โปรแกรมต้องการ

ผมลองดับเบิ้ลคลิกไฟล์ “surveyor pocket tools.exe” ก็สามารถเปิดมาและทำงานได้ตามปกติ ลองดูชื่อไดเรคทอรีจะเห็นสองไดเรคทอรี ที่ได้จากไฟล์ hooks คือ openpyxl และ pyproj ลองเข้าไปดูในไดเรคทอรี จะเห็นข้อมูลที่ pyproj ขนไปใช้ หมายเหตุว่าข้อมูลนี้ pyproj จะนำไปเป็นฐานข้อมูลในการแปลงพิกัดตาม datum และ projection

ทำไฟล์ Setup ด้วย Inno Setup

จากนั้นผมจะ copy ไดเรคทอรีที่อยู่ใน “setup” ไปไว้อีกที่หนึ่ง พื้นที่นี้สำหรับใช้ Inno Setup มาทำไฟล์ติดตั้ง ลองดูไดเรคทอรี

ในไดเรคทอรีนี้ผมจะมีไฟล์ “surveyorpockettools64.iss” เป็นไฟล์สคริปท์ของ Inno Setup เพื่อสร้างไฟล์ติดตั้ง setup สำหรับวินโดส์ 64 บิต

#define MyAppName "Surveyor Pocket Tools"
#define MyAppEXE "Surveyor Pocket Tools.exe"
#define MyShortAppName "SurveyorPocketTools"
#define MyMainRoot "Survey Suite"
#define Developer "Prajuab Riabroy"
#define Version "1.02"
#define Build "641"

[Setup]
AppName={#MyAppName}
AppVerName={#MyAppName} V{#Version}
DefaultDirName={pf}\{#MyMainRoot}\{#MyAppName}
DefaultGroupName={#MyMainRoot}\{#MyAppName}
UseSetupLdr=yes
UninstallDisplayIcon={app}\{#MyAppEXE}
VersionInfoProductName={#MyAppName}
VersionInfoCompany=priabroy
VersionInfoCopyright=Copyright 2000-2017 by {#Developer}
VersionInfoDescription={#MyAppName}
VersionInfoProductVersion={#Version}
VersionInfoVersion={#Version}
OutputDir=Setup
OutputBaseFilename={#MyShortAppName}V{#Version}Build{#Build}Setup64
;OutputDir=TraverseProV250Setup64
; "ArchitecturesAllowed=x64" specifies that Setup cannot run on
; anything but x64.
ArchitecturesAllowed=x64
; "ArchitecturesInstallIn64BitMode=x64" requests that the install be
; done in "64-bit mode" on x64, meaning it should use the native
; 64-bit Program Files directory and the 64-bit view of the registry.
;ArchitecturesInstallIn64BitMode=x64
ArchitecturesInstallIn64BitMode=x64
AppPublisher={#Developer}
AppPublisherURL=https://www.surveyorpockettools.org
AppVersion={#Version}.{#Build}
LicenseFile = eula.txt
ChangesEnvironment=yes
SolidCompression=yes
Compression=lzma2/ultra64
LZMAUseSeparateProcess=yes
LZMADictionarySize=1048576
LZMANumFastBytes=273

[Files]
Source: "{#MyAppName}.exe"; DestDir: "{app}"
Source: "qtwebengineprocess.exe"; DestDir: "{app}"
Source: "base_library.zip"; DestDir: "{app}" ;
Source: "qt5_plugins\*"; DestDir: "{app}\qt5_plugins\"; Flags: ignoreversion recursesubdirs
Source: "certifi\*"; DestDir: "{app}\certifi\"; Flags: ignoreversion recursesubdirs
Source: "PySide2\*"; DestDir: "{app}\PySide2\"; Flags: ignoreversion recursesubdirs
Source: "database\*"; DestDir: "{userappdata}\{#MyAppName}\database\";
Source: "markers\*"; DestDir: "{userappdata}\{#MyAppName}\markers\";
Source: "example data\*"; DestDir:"{userappdata}\{#MyAppName}\example data\";
Source: "pyproj\data\*"; DestDir: "{userappdata}\pyproj\data\";
Source: "pyproj\*.pyd"; DestDir: "{app}\pyproj\";
Source: "lxml\*"; DestDir: "{app}\lxml\";
Source: "shiboken2\*"; DestDir: "{app}\shiboken2\";
Source: "*.pak"; DestDir: "{app}";
Source: "*.dll"; DestDir: "{app}"
Source: "*.pyd"; DestDir: "{app}"
Source: "*.xml"; DestDir: "{userappdata}\{#MyAppName}\";
Source: "qt.conf"; DestDir: "{app}"
Source: "icudtl.dat"; DestDir: "{app}"

[Icons]
;create icon at start menu group
Name: "{group}\{#MyAppName}"; Filename: "{app}\{#MyAppExe}"
;create icon at desktop
Name: "{commondesktop}\{#MyAppName}"; FileName:"{app}\{#MyAppExe}"

[Registry]
; Start "Software\My Company\My Program" keys under HKEY_CURRENT_USER
; and HKEY_LOCAL_MACHINE. The flags tell it to always delete the
; "My Program" keys upon uninstall, and delete the "My Company" keys
; if there is nothing left in them.
Root: HKCU; Subkey: "Software\{#MyMainRoot}"; Flags: uninsdeletekeyifempty
Root: HKCU; Subkey: "Software\{#MyMainRoot}\{#MyAppName}"; Flags: uninsdeletekey
Root: HKLM; Subkey: "Software\{#MyMainRoot}"; Flags: uninsdeletekeyifempty
Root: HKLM; Subkey: "Software\{#MyMainRoot}\{#MyAppName}"; Flags: uninsdeletekey
Root: HKLM; Subkey: "Software\{#MyMainRoot}\{#MyAppName}\Settings"; ValueType: string; ValueName: "InstalledPath"; ValueData: "{app}"
Root: HKLM; Subkey: "Software\{#MyMainRoot}\{#MyAppName}\Settings"; ValueType: string; ValueName: "DevelopedBy"; ValueData: "{#Developer}"
Root: HKLM; Subkey: "Software\{#MyMainRoot}\{#MyAppName}\Settings"; ValueType: string; ValueName: "ApplicationName"; ValueData: "{#MyAppName}"
;Root: HKCU; Subkey: "Environment"; ValueType:string; ValueName:"PROJ_LIB"; ValueData:"{userappdata}\{#MyAppName}\geoidgrids\" ; Flags: preservestringtype ;

ถ้าเป็นไฟล์สำหรับ Surveyor Pocket Tools รุ่น 32 บิตเพียงใส่คอมเมนต์หน้า ;ArchitecturesInstallIn64BitMode=x64 ก็พอ สำหรับรายละเอียดสคริปท์ของ Inno Setup ผมจะไม่กล่าวถึงรายละเอียดในที่นี้ผู้อ่านที่สนใจสามารถศึกษาได้ครับ จากนั้นก็ใช้ Inno Setup ทำการ build ก็จะได้ไฟล์ Exe เดี่ยวๆ ที่สามารถ zip ไปให้ผู้ใช้ได้ download ต่อไป พบกันตอนหน้าครับ