Tag: คำนวณ

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

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

ไม่นานมานี้มีผมดาวน์โหลดไฟล์โปรแกรมและข้อมูลของ TGM2017 เรียกเต็มๆคือ Thailand Geoid Model 2017 ที่เป็นโครงการร่วมมือจากหลายๆฝ่ายของทางราชการ ผมยังไม่มีโอกาสได้นำไปใช้งาน โดยเฉพาะจะนำมาประยุกต์ใช้ในงานรังวัด GNSS ถือว่าเป็นสิ่งที่พวกเรารอคอยมานานที่จะได้มี local geoid model มาใช้งานกัน

โดยเฉพาะงานรังวัด GNSS เมื่อคำนวณแล้วจะได้ค่าพิกัดทางราบ และทางดิ่งจะได้ความสูง (h – Ellipsoid height) ที่เทียบกับทรงรี WGS84 ถ้าเรามีแบบจำลองความสูงจีออยด์สามารถคำนวณหาความสูง (H – Orthometric height) หรือระดับน้ำทะเลปานกลางได้ จากสูตร h = H + N โดยที่ไม่ต้องมีการเดินระดับไปหาจุดที่รังวัด GNSS แต่อย่างใด โดยที่โปรแกรมคำนวณหาค่าความสูงจีออยด์ (N – Geoid Height หรือ Geoid Separation) ที่ผมจะมาลองเขียนโปรแกรมเล็กๆเพื่อหาค่า N นี้ด้วยภาษาไพทอนในลำดับต่อไป

ในตอนนี้ผมจะใช้ไลบรารี Proj.4 เรียกไลบรารีในเวอร์ชั่นภาษาไพทอนนี้ว่า PyProj ตัวไลบรารี Proj.4 เองก็เพิ่งสนับสนุนการคำนวณ Vertical Datum ไปไม่นานนี้เอง สมมติว่าถ้าไม่ได้ใช้ไลบรารีตัวนี้เป็นตัวช่วยคงเป็นงานหนักหนาสาหัสเหมือนกัน

แปลงรูปแบบไฟล์ของแบบจำลองความสูงจีออยด์ TGM2017

ไฟล์แบบจำลองความสูงจีออยด์ TGM2017 ที่ผมดาวน์โหลดได้มานั้นเป็นแอสกี้ (Text file) ขนาดประมาณ 10.7 MB และมีโปรแกรมคำนวณมาด้วยชื่อไฟล์ execute คือ “TGM2017.EXE” สำหรับไฟล์แอสกี้ในคู่มือแสดงรายละเอียดดังนี้

  • ค่าละติจูด (Latitude) เริ่มต้น : 3 องศา
  • ค่าลองจิจูด (Longitude) เริ่มต้น : 95 องศา
  • ความละเอียดเชิงพื้นที่ : 1 ลิปดา
  • จานวนหลัก (Column) : 780 หลัก
  • จานวนแถว (Row) : 1200 แถว
TGM2017 layout

สำหรับไลบรารี PyProj นั้นถ้าต้องการคำนวณ Vertical Datum ต้องการไฟล์ในรูปแบบ gtx ที่กำหนดมาตรฐานโดย NOAA เพื่อให้การใช้งานมีประสิทธิภาพมาก ทางนี้แนะนำให้ใช้ไฟล์ไบนารีเพื่อประหยัดเมมโมรี สำหรับรูปแบบไฟล์ gtx นั้นไม่มีอะไรมาก กำหนดให้มีไฟล์ header เพื่อบอกจุดเริ่มต้นขอบเขตการใช้งาน ระยะห่างของละติจูดและลองจิจูดของแต่ละจุดในแนวตั้งและแนวนอน จำนวนแถวและจำนวนคอลัมน์ ที่เหลือจะเป็นความสูงจีออยด์ในแต่ละจุดมีขนาดเท่ากับ จำนวณแถวคูณจำนวนคอลัมน์ ตัวหัวไฟล์ผมสามารถเขียนได้ดังนี้

3.0 95.0 0.01666666666 0.01666666666 1200 780

อธิบาย : ละติจูดของมุมล่างซ้าย ลองจิจูดของมุมล่างซ้าย ระยะห่างของจุดในแนวละติจูด ระยะห่างของจุดในแนวลองจิจูด จำนวนแถว จำนวนคอลัมน์

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

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

ผมเขียนสคริปต์ภาษาไพทอน เผื่อมีใครสนใจ ผมตึ้งชื่อไฟล์ว่า “rtsd2gtx.py” โปรแกรมมีขนาดเล็กมากๆ ไม่กี่บรรทัด หลักการคืออ่านไฟล์แอสกี้ “tgm2017.asc” มาแล้วเขียนไฟล์ไบนารีในรูปแบบ gtx ข้อควรระวังคือไฟล์ไบนารี ต้องเป็นแบบ Big Endian ในโมดูล struct ผมจึงใส่เครื่องหมาย “>” ไปด้วย

ตัวอย่างเขียนขนิดข้อมูลเป็น double ใช้คำสั่ง f.write(struct.pack(“>d”, 3.0))
เขียนข้อมูลเป็น integer ใช้คำสั่ง f.write(struct.pack(“>i”, 1200))

<import struct
def readRTSDAsciiFile():
    geoid = []
    with open('tgm2017.asc', 'r') as f:
        geoid = [[float(n) for n in line.split()] for line in f] #2D list, one list contain one row.
    #TGM2017 was organized already from lowest row (min latitude) to highest row (max latitude).
    #No need to reverse.	
    return geoid

def writeGTXBinaryFile(geoid):
    with open('tgm2017.gtx', 'wb') as f:
        # pack with ">" for big endian.
		# write Header of file
        f.write(struct.pack(">d", 3.0))  #lower left Latitude
        f.write(struct.pack(">d", 95.0)) #lower left Longtitude		
        f.write(struct.pack(">d", 1.0/60.0)) #delta Latitude = 1 second
        f.write(struct.pack(">d", 1.0/60.0)) #delta Longitude = 1 second	
        f.write(struct.pack(">i", 1200)) #number of rows
        f.write(struct.pack(">i", 780))	#number of columns
		# write height values.
        for row in geoid:
            for c in row:
                f.write(struct.pack(">f", c))

# main program	
g = readRTSDAsciiFile()
writeGTXBinaryFile(g)

ณ ตอนนี้ผมสร้างโฟลเดอร์ เอาไฟล์ tgm2017.asc มาวางไว้ และมีโปรแกรมสคริปต์ “rtsd2gtx.py” อยู่ในโฟลเดอร์เดียวกัน

โฟลเดอร์ที่เก็บไฟล์

จากนั้นทำการรันโปรแกรมด้วยคำสั่งที่คอมมานด์ไลน์ python rtsd2gtx.py รอแป๊ปหนึ่งไฟล์ gtx ที่เก็บแบบจำลองความสูงจีออยด์จะถูกสร้าง และมีขนาดเล็กลงจากไฟล์แอสกี้ดั้งเดิมประมาณ 3 เท่า

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

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

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

โค้ดโปรแกรมเริ่มต้นจากสร้าง projection ต้นทาง (ตัวแปร wgs84 egm) โดยที่กำหนด +geoidgrids=tgm2017.gtx และสร้าง projection ปลายทาง (ตัวแปร wgs84) โดยที่เว้นไม่ใส่แบบจำลองความสูงจีออยด์

import pyproj # Import the pyproj module
# Define source coordinate system with TGM2017.
wgs84egm = pyproj.Proj("+init=EPSG:4326 +geoidgrids=tgm2017.gtx") 
# Define target coordinate system without geoid model.
wgs84 = pyproj.Proj("+init=EPSG:4326")

lat, lon = 18.3353398540, 99.371215293000
x, y, z = pyproj.transform(wgs84egm, wgs84, lon, lat, 0.0)
print("1) Latitude = %.8f, Longitude = %.8f, Geoid Height = %.5f" % (y, x, z))

lat, lon = 7.75906518100, 98.303594340000
x, y, z = pyproj.transform(wgs84egm, wgs84, lon, lat, 0.0)
print("2) Latitude = %.8f, Longitude = %.8f, Geoid Height = %.5f" % (y, x, z))

lat, lon = 9.18558870000, 99.843706600000
x, y, z = pyproj.transform(wgs84egm, wgs84, lon, lat, 0.0)
print("3) Latitude = %.8f, Longitude = %.8f, Geoid Height = %.5f" % (y, x, z))

lat, lon = 6.96243169, 99.77398865
x, y, z = pyproj.transform(wgs84egm, wgs84, lon, lat, 0.0)
print("4) Latitude = %.8f, Longitude = %.8f, Geoid Height = %.5f" % (y, x, z))

lat, lon = 7.20852844, 99.72353542
x, y, z = pyproj.transform(wgs84egm, wgs84, lon, lat, 0.0)
print("5) Latitude = %.8f, Longitude = %.8f, Geoid Height = %.5f" % (y, x, z))

lat, lon = 7.35631365, 100.12797020
x, y, z = pyproj.transform(wgs84egm, wgs84, lon, lat, 0.0)
print("6) Latitude = %.8f, Longitude = %.8f, Geoid Height = %.5f" % (y, x, z))

lat, lon = 7.23806984, 100.55255020
x, y, z = pyproj.transform(wgs84egm, wgs84, lon, lat, 0.0)
print("7) Latitude = %.8f, Longitude = %.8f, Geoid Height = %.5f" % (y, x, z))

lat, lon = 7.09316728, 100.56368280
x, y, z = pyproj.transform(wgs84egm, wgs84, lon, lat, 0.0)
print("8) Latitude = %.8f, Longitude = %.8f, Geoid Height = %.5f" % (y, x, z))

lat, lon = 6.99643503, 100.43034140
x, y, z = pyproj.transform(wgs84egm, wgs84, lon, lat, 0.0)
print("9) Latitude = %.8f, Longitude = %.8f, Geoid Height = %.5f" % (y, x, z))

lat, lon = 13.0+42.0/60+38.64/3600, 100+29.0/60+13.02/3600
x, y, z = pyproj.transform(wgs84egm, wgs84, lon, lat, 0.0)
print("10) Latitude = %.8f, Longitude = %.8f, Geoid Height = %.5f" % (y, x, z))

เมื่อจัดเก็บไฟล์เรียบร้อย ในตอนนี้ ถ้าดูที่โฟล์เดอร์จะเห็นไฟล์ดังนี้

จุดทดสอบ

ถ้าดูในโค้ดจะเห็นจุดทดสอบทั้งหมด 10 จุด โดยที่ 9 จุดแรกผมนำมาจากไฟล์ INPUT.DAT ที่ติดมากับโปรแกรม TGM2017.EXE ถ้าเปิดดูด้วย Notepad จะเห็นดังนี้

จุดทดสอบ

ส่วนจุดที่ 10 ได้มาจากสไลด์ ชื่อ BM11SP Latitude = 13° 42′ 38.64″ N Longitude = 100° 29′ 13.02″E

คำนวณ

จากนั้นทำการรันโปรแกรมบนคอมมานด์ไลน์ ใช้้คำสั่ง python test.py จะได้ผลลัพธ์ดังนี้

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

ผมนำผลลัพธ์มาเปรียบเทียบกับที่คำนวณด้วย TGM2017.EXE พบว่ามีความต่างกันเล็กน้อยที่หลักมิลลิเมตร สาเหตุน่าจะเกิดตอน interpolation โดยการนำจุดรอบๆจากแบบจำลองความสูงจีออยด์ ซึ่งโปรแกรมอาจจะเลือกใช้แบบ BiCubic หรือ BiLinear ผมไม่แน่ใจว่า Proj.4 ใช้แบบไหน ตอนนี้ยังไม่มีโอกาสไปดูโค้ดต้นฉบับ

เปรียบเทียบผลลัพธ์การคำนวณ

การพัฒนาโปรแกรมสำหรับการใช้งานในอนาคต

ผมวางแผนว่าจะนำแบบจำลองความสูงจีออยด์ TGM2017 ไปคำนวณในโปรแกรมรวมเครื่องมือฉบับกระเปาสำหรับช่างสำรวจ Surveyor Pocket Tools ได้แก่การคำนวณความสูงจีออยด์ การคำนวณหาสเกลแฟคเตอร์

สรุป

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


Print Preview ใน Excel

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

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

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

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

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

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

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

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

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

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

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

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

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

 ตั้งสมการ Observation Equation

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

v + Bx = f

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Iterations = 2

Number of Stations = 7

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

 

Adjusted Coordinates (Meters)

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

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

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

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

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

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

การวัดแบบ 2D

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

การวัดแบบ 3D

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

แนะนำตำรา

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

 

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

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

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

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

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

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

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

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

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

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

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

Kaestner-Burkhardt Method

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

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

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

 

Cassini Method

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

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

 

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

Collins Method

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

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

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

 

 

Font-Llagunes

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

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

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

 

 

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

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

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

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

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

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

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

 

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

 

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

 

ติดปีกเครื่องคิดเลขเทพ Casio fx 9860G II SD ด้วยโปรแกรมภาษาซีบน AddIn ตอนที่ 9 โปรแกรมคำนวณหาจุดตัด (Intersection)

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

โปรแกรมคำนวณหาจุดตัด (Intersection)

โปรแกรมเขียนไว้นานแล้ว แต่มีโอกาสมาปรับแก้ใหม่ให้สามารถคำนวณจุดสัมผัส เช่นเส้นตรงกับวงกลมโดยสามารถตั้ง Tolerance ได้ ผมตั้งไว้ที่ 1 มม. ขอขยายความว่าเราสามารถหาจุดตัดเส้นตรงกับวงกลมได้เป๊ะๆ สองจุด แต่ในกรณีที่เป็นจุดสัมผัส (Tangent) นั้นยากโอกาสน้อยมากที่จะคำนวณหาจุดสัมผัสจากโจทย์ที่กำหนดเส้นตรงที่ผ่านจุดและมีทิศทางตามอะซิมัท มาสัมผัสกับวงกลมที่กำหนดจุดศูนย์กลางและรัศมี  ผมจึงตั้ง tolerance ไว้ 1 มม. ถ้าเส้นตรงมาเฉียดจุดสัมผัสนี้ไม่ว่าจะด้านนอกหรือด้านในวงกลมถ้าระยะห่างเส้นตรงห่างจากจุดสัมผัสจริงน้อยกว่า 1 มม. ก็ให้ถือว่าเส้นตรงเส้นนี้สัมผัสกับวงกลม สำหรับไอคอนโปรแกรมเมื่อติดตั้งแล้ว ดังรูปด้านล่าง

ดาวน์โหลดและติตตั้งบนเครื่องคิดเลข

ไปที่หน้าดาวน์โหลด (Download) มองหาโปรแกรมบนเครื่องคิดเลข Casio fx-9860G II SD  ชื่อโปรแกรม Intersection จากนั้นทำการดาวน์โหลดมาจะได้ไฟล์ชื่อ “INSCTEX.G1A” แล้วทำการ copy ไฟล์ตัวนี้ไปยังเครื่องคิดเลขด้วยโปรแกรม Casio FA-124  หรือ copy ผ่านทางตัว SD Card ที่มากับเครื่องคิดเลข

เริ่มต้นใช้งานโปรแกรมหาจุดตัด (Intersection)

กดคีย์ “Main Menu” ของเครื่องคิดเลขใช้คีย์ลูกศรไล่ไปหาไอคอนโปรแกรม Intersection ดังรูป แล้วกดคีย์ “EXE”

เมนูหลักของโปรแกรม

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

F1 – Set สำหรับเลือกการคำนวณว่าจะเป็นเส้นตรงตัดกับเส้นตรง เส้นตรงกับวงกลม หรือวงกลมกับวงกลม

F2 – IN (Input) จะเป็นการป้อนค่าพิกัดหรืออะซิมัท ตามตัวเลือกที่เลือกไว้ในเมนูแรก (F1 – Set)

F3 – Calc คำนวณหาจุดตัด ซึ่งในเมนูย่อยสามารถแสดงรูปเส้นตรง (Plot) วงกลมพร้อมจุดตัดได้

F5 – Info แสดงเครดิตไลบรารี ที่โปรแกรมนำมาใช้งาน

F6-Exit ออกจากโปรแกรม

ตัวอย่างที่ 1 คำนวณหาจุดตัดระหว่างเส้นตรงกับเส้นตรง (Line and Line Intersection)

เป็นการคำนวณหาจุดตัดเมื่อกำหนดเส้นตรงสองเส้นที่มีค่าพิกัดหัวและท้าย ดังรูปด้านล่างเป็นแปลงที่ดิน DEFG ต้องการหาค่าพิกัดจุด “X” ซึ่งเป็นพิกัดศุนย์กลางพื้นที่ที่เกิดจากเส้นตรงจากมุมของแปลงที่ดินลากเป็นเส้นทแยงมุมตัดกัน

เลือกรายการคำนวณ

ที่เมนูหลักของเครื่องคิดเลขกดคีย์ F1-Set จะเห็นตัวเลือกรูปแบบการคำนวณจุดตัด (Intersection Type) เมื่อกดคีย์ “EXE” จะเห็นรายการให้เลือก 4 อย่าง ให้เลือก “Line X Line” ดังรูป กดคีย์ “EXE” จากนั้นกดคีย์ F6-OK เพื่อออก

ป้อนค่าพิกัด

กลับมาที่เมนูหลักของโปรแกรม กดคีย์ F2-IN ป้อนค่าพิกัดของเส้นตรงสองเส้น ตามโจทย์ข้างต้น กดคีย์ F6 – OK เพื่อออก

คำนวณหาจุดตัด

กลับมาที่เมนูหลักของโปรแกรม จากนั้นกดคีย์ F3 – Calc เพื่อคำนวณหาจุดตัด “X” โปรแกรมจะทวนค่าพิกัดของเส้นตรงเส้นที่ป้อนไว้ก่อนจะแสดงค่าพิกัดจุดตัด ถ้าจุดนี้อยู่บนเส้นตรงทั้ง 2 เส้น กดคีย์ F2 – PgDn เพื่อเลื่อนไปหน้าต่อไป

หน้าสุดท้ายจะเห็นค่าพิกัดจุดตัด “X” N = 1652.560, E = 1739.142

แสดงรูปตัด (Plot)

กดคีย์ F5-Plot เพื่อแสดงเส้นตรงสองเส้นและจุดตัด จอภาพเป็น dot pixel หยาบๆก็ได้ขนาดนี้พอให้เกิดจินตนาการว่าเส้นตรง 2 เส้นวางตัวในลักษณะใดและตัดกันที่ตรงไหน กดคีย์ F6-Done เพื่อออก

ตัวอย่างที่ 2 คำนวณหาจุดตัดระหว่างเส้นตรงกับเส้นตรง (Azimuth and Azimuth Intersection)

การกำหนดเส้นตรงโดยกำหนดจุดและอะซิมัทให้ ในกรณีนี้เส้นตรงจะยาวไม่สิ้นสุด ดังนั้นถ้าเส้นตรงที่กำหนดไม่ขนานกัน ก็มั่นใจได้เลยว่าหาจุดตัดได้แน่ (ต่างจากตัวอย่างแรกที่หาจุดตัดบนเส้นตรงทั้งสองเท่านั้น)

กำหนดเส้นตรงที่เป็น Alignment ของงานถนน เส้นที่ 1 เป็น  ผ่านจุด N: 2641990.928 , E = 231848.514  Azimuth: 35°28′ 1.7433″ เส้นที่ 2 ผ่านจุด N: 2641812.446 E: 231753.041 Azimuth: 9°12′ 20.4212″ คำนวณหาจุดตัด PI (Point of Intersection)

เลือกรายการคำนวณ

ที่เมนูหลักของโปรแกรม กดคีย์ F1-Set เลือกรายการคำนวณ “Azi X Azi” กดคีย์ F6-OK เพื่อออก

ป้อนค่าพิกัด

กลับมาที่เมนูหลักโปรแกรม กดคีย์ F2-IN เพื่อป้อนค่าพิกัด ไม่ลืมว่าป้อนอะซิมัทคั่นด้วยเครื่องหมายลบ  เมื่อเสร็จแล้วกดคีย์ F6-OK เพื่อออก

คำนวณหาจุดตัดและแสดงรูป

ที่เมนูหลักของโปรแกรม กดคีย์ F3-Calc เพื่อคำนวณหาจุดตัด โปรแกรมจะทวนค่าพิกัดที่ป้อนเข้าไปก่อน กดคีย์ F2-PgDn เพื่อเลื่อนไปหน้าแสดงผลถัดไป

จะได้ค่าพิกัดจุดตัดที่เป็นจุด PI (Point of Intersection N: 2641870.013, E: 231762.371 กดคีย์ F5-Plot จะเห็นรูปร่างเส้นตรงสองเส้น พร้อมจุดตัด “I1” บนหน้าจอภาพกดคีย์ F6-Done สองครั้งเพื่ออก

ตัวอย่างที่ 3 คำนวณหาจุดตัดระหว่างเส้นตรงกับวงกลม (Azimuth and Circle Intersection)

อย่างที่ผมเกริ่นไปข้างต้นว่าเส้นตรงถ้าตัดผ่าเข้าไปในวงกลมจะได้จุดตัดสองเส้น แต่กรณีพิเศษที่เส้นตรงไปสัมผัสกับวงกลมกรณีนี้จะได้จุดตัดคือจุดสัมผัสมาจุดเดียว มาลองดูตัวอย่างแบบนี้ กำหนดให้เส้นตรงผ่านจุด N: 2642178.562, E: 231597.085 Azimuth: 161°8′ 58.2981″ กำหนดวงกลมมีจุดศูนย์กลาง N: 2641772.451, E: 231999.821  รัศมี 249.921 เมตร

เลือกรายการคำนวณ

ที่เมนูหลักของโปรแกรม กดคีย์ F1-Set เลือกรายการคำนวณ “Azi X Cir” กดคีย์ F6-OK เพื่อออก

ป้อนค่าพิกัดและรัศมีวงกลม

ที่เมนูหลักของโปรแกรม กดคีย์ F2-IN เพื่อป้อนค่าพิกัดของเส้นตรงและป้อนค่าอะซิมัท ป้อนค่าพิกัดศูนย์กลางวงกลมพร้อมทั้งรัศมี เสร็จแล้วกดคีย์ F6-OK เพื่อออก

คำนวณจุดตัดและแสดงรูป

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

สังเกตว่าตรงจุดตัดโปรแกรมจะแสดงว่าได้จุดสัมผัส Tangent มา 1 จุดคือ N: 2641691.702, E: 231763.305 กดคีย์ F5-Plot เพื่อแสดงรูปเส้นตรงและวงกลม ก็พอกล้อมแกล้มการวาดรูปร่างวงกลม ไม่มีไลบรารีของเครื่องคิดเลขต้องอาศัยวาดจุดลงไปตามเส้นรอบวงแทน  จุดสัมผัสแสดงด้วยตัวอักษร “I1” จุดศูนย์กลางวงกลมแสดงด้วยตัวอักษร “C1”

ตัวอย่างที่ 4 คำนวณหาจุดตัดระหว่างวงกลมกับวงกลม (Circle and Circle Intersection)

มาถึงตัวอย่างสุดท้าย ผมจะขอรวบรัดแสดงเฉพาะรูปหน้าจอ กำหนดโจทย์ วงกลมวงแรกมีค่าพิกัดศูนย์กลาง N: 2641210.885, E: 232480.916 รัศมี 525 เมตร วงกลมวงที่ 2 N: 2641256.635 E: 233130.568 รัศมี 250 เมตร คำนวณหาจุดตัดระหว่างวงกลมสองวงนี้

จัดเก็บข้อมูลและเรียกมาใช้ภายหลัง

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

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

จากนั้นนำ SDCard มาเสียบบนเครื่องคิดเลขอีกครั้ง เมื่อนำไปใช้งานได้สักพักถ้าเอามาเปิดอีกครั้งจะเห็นไฟล์หลายๆไฟล์ มีนามสกุลเป็น “CFG”  หมายถึง config ตัวอย่างถ้าใช้โปรแกรมคำนวณหาจุดตัดนี้ไฟล์ที่จัดเก็บข้อมุลคือ “INTSCT.CFG

ก็ยังมีหลายซีรี่ย์สำหรับโปรแกรมเครื่องคิดเลขในชุดนี้ ก็ติดตามกันได้ต่อไปครับ

Update: โปรแกรมแปลงค่าพิกัดภูมิศาสตร์ Geographic Calculator (GeoCalc) บนเครื่องคิดเลข Casio fx-9860G II SD

Geographic Calculator

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

เปลี่ยนรูปแบบการติดต่อกับผู้ใช้ด้วยไลบรารี MyLib

ไลบรารี MyLib เป็นไลบรารีภาษาซีเล็กๆที่ผู้พัฒนาใช้นาม hayzel ได้เขียนไว้เพื่อใช้บนเครื่องคิดเลข Casio fx-9860G II SD ผมนำมาใช้และชอบ ทำให้มีแรงใจที่จะเขียนโปรแกรมบนเครื่องคิดเลขรุ่นเทพรุ่นนี้ได้มาหลายโปรแกรม และก็เหมือนเดิมครับว่าโปรแกรมที่ผมเขียนนั้นใช้งานได้ฟรี (Freely Usability) เพื่อใช้ในแวดวงงานสำรวจทำแผนที่ตลอดจนงานสำรวจเพื่อการก่อสร้างก็ตาม

ก็ขอตั้งชื่อโปรแกรมเล็กๆสำหรับแปลงพิกัดบนพื้นหลักฐาน WGS84 นี้ว่า “GeoCalc Extra” ก่อนอื่นสูตรที่ใช้ในการคำนวณผมใช้ไลบรารีชื่อ mgrs สามารถคำนวณแปลงพิกัดในระบบพิกัด UTM, geographic, UPS และ MGRS ได้ ขนาดไม่ใหญ่มากนัก สามารถคอมไพล์และบิวท์มาใส่เครื่องคิดเลขรุ่นนี้ได้ ทั้ง MyLib และ mgrs เป็นโปรแกรมเปิดโค้ด ฟรีทั้งคู่

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

ไปที่หน้าดาวน์โหลด (Download) มองหาโปรแกรมบนเครื่องคิดเลข Casio fx-9860G II SD  ชื่อโปรแกรม GeoCalc Extra จากนั้นทำการดาวน์โหลดมาจะได้ไฟล์ชื่อ “GEOCALC.G1A” แล้วทำการ copy ไฟล์ตัวนี้ไปยังเครื่องคิดเลขด้วยโปรแกรม Casio FA-124  หรือ copy ผ่านทางตัว SD Card ที่มากับเครื่องคิดเลข

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

ที่ “Main Menu” ของเครื่องคิดเลขเลื่อนไปหาไอคอนของโปรแกรมดังรูปด้านล่าง เมื่อกดคีย์ “EXE” แล้วจะเข้าเมนูของโปรแกรม ดังนี้

F1 – Set เลือกรายการคำนวณแปลงพิกัดระหว่าง UTM, Geographic หรือ MGRS

F2 – Src (Source) ป้อนค่าพิกัดที่ต้องการแปลงพิกัด

F3 – Calc คำนวณแปลงพิกัดพร้อมแสดงผลลัพธ์

F5 – Info แสดงเครดิตไลบรารีที่โปรแกรมใช้งาน

F6 – Exit ออกจากโปรแกรม

เลือกรายการคำนวณ (Menu)

ที่เมนูหลักกดคีย์ F1 – Set เพื่อเข้าไปเลือกรายการคำนวณ จะเห็น ระบบพิกัดเริ่มต้น (Source)   และระบบพิกัดปลายทาง (Target) ส่วนด้านล่าง MGRS Precision จะเป็นความละเอียดของระบบพิกัด MGRS (Military Grid Reference System) เลือกได้ 6 ระดับคือ 0, 2, 4, 6, 8, 10

ตัวอย่างที่ 1 แปลงค่าพิกัดจากค่าพิกัดภูมิศาสตร์ไปยังค่าพิกัดยูทีเอ็ม (Geographic to UTM)

ตั้งค่าระบบพิกัดเริ่มต้นและปลายทาง

ตั้งค่าระบบพิกัดต้นทางและปลายทางดังรูป จากนั้นกดคีย์ F6 – OK เพื่อออก

ป้อนค่าพิกัด

กลับมาที่เมนูหลักของโปรแกรมอีกครั้ง กดคีย์ F2 – Src เพื่อป้อนค่าพิกัดภูมิศาสตร์

ป้อนค่าพิกัดละติจูด Latitude 39°57’9.34803″N โดยการป้อน 39-57-9.34803N ค่าลองจิจูด 75°9’54.75490″W ป้อนค่า 75-9-54.75490W เสร็จแล้วกดคีย์ F6 – OK เพื่อออกไปคำนวณ

คำนวณแปลงพิกัด

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

จะแสดงค่าพิกัดเริ่มต้นให้และค่าพิกัดปลายทางคือยูทีเอ็มให้ พร้อมทั้งบอกโซนของยูทีเอ็มให้ กดคีย์ F6 – Done เพื่อออก

ตัวอย่างที่ 2 แปลงค่าพิกัดจากค่าพิกัดยูทีเอ็มไปยังค่าพิกัด MGRS (UTM to MGRS)

กำหนดค่าพิกัดยูทีเอ็ม (UTM) ดังนี้ N: 2642783.110, E: 232030.949 UTM Zone No: 46N กลับมาที่เมนูหลักของโปรแกรม กดคีย์ F1 – Set เพื่อเปลี่ยนรายการคำนวณ ตั้งค่าตามรูปด้านล่าง 

เสร็จแล้วกดคีย์ F6 – OK เพื่อออก

กลับมาเมนูหลักของโปรแกรมกดีย์ F2 – Src เพื่อออกป้อนค่าพิกัดยูทีเอ็มดังนี้ จากนั้นกดคีย์ F6 – OK เพื่อออก

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

ตัวอย่างที่ 3 แปลงค่าพิกัดจากค่าพิกัด MGRS ไปยังค่าพิกัดภูมิศาสตร์ (MGRS to Geographic)

กำหนดค่าพิกัด MGRS: 46QCK0907425049 ส่วนขั้นตอนจะขอรวบรัดแสดงด้วยรูปภาพ

จัดเก็บข้อมูลและเรียกมาใช้ภายหลัง

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

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

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

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

ติดปีกเครื่องคิดเลขเทพ Casio fx 9860G II SD ด้วยโปรแกรมภาษาซีบน AddIn ตอนที่ 8 โปรแกรมคำนวณสเกลแฟคเตอร์ (Scale Factor)

ติดปีกเครื่องคิดเลขเทพ Casio fx 9860G II SD ด้วยโปรแกรมภาษาซีบน AddIn ตอนที่ 8 โปรแกรมคำนวณสเกลแฟคเตอร์ (Scale Factor)

โปรแกรมคำนวณสเกลแฟคเตอร์ (Scale Factor) สำหรับเครื่องคิดเลข Casio fx-9860 G

ช่วงนี้อยู่ในชุดซีรี่ย์โปรแกรมเครื่องคิดเลข Casio fx-9860G  ต่อไปขอนำเสนอโปรแกรมคำนวณสเกลแฟคเตอร์ (Scale Factor) ตัวโปรแกรมพัฒนาด้วยภาษาซี  ใช้เป็นโปรแกรม AddIn ดังรูปด้านล่าง

ปัญหาการคำนวณสเกลแฟคเตอร์

สเกลแฟคเตอร์ในที่นี้ผมจะขอหมายถึง Elevation Scale Factor (ESF), Grid Scale Factor (GSF) ตลอดจน Combined Scale Factor (CSF) สามารถอ่านได้จาก blog ก่อนหน้านี้ของผมได้ที่ใช้โปรแกรม Surveyor Pocket Tools มาช่วยคำนวณ ผู้อ่านอาจจะเคยสังเกตว่าจะไม่เคยเห็นโปรแกรมคำนวณสเกลแฟคเตอร์ในเครื่องคิดเลขเท่าใดนัก ในกรณีนี้คือพวกเราใช้ค่าระดับที่อิงระดับน้ำทะเลปานกลาง (แทนด้วยสัญลักษณ์ตัว “H“) กันเป็นหลัก น้ำทะเลปานกลางนี้สามารถใช้พื้นผิวจีออยด์ (Geoid) แทนได้ แต่การคำนวณหา Elevation Scale Factor เราจะต้องทราบค่าระดับที่อ้างอิงกับ Ellipsoid (แทนด้วยสัญลักษณ์ตัว “h“) ไม่ใช่เทียบกับระดับน้ำทะเลปานกลาง ความต่างระหว่างพื้นผิวจีออยด์และพื้นผิวทรงรี เราเรียกว่า Geoid Separation (แทนด้วยสัญลักษณ์ตัว “N“) ความสัมพันธ์ระหว่างสามสิ่งนี้คือ h = H + N

การหา Geoid Separation (N) จะต้องใช้ Earth Gravity Model ซึ่งก่อนหน้านี้ใช้ EGM96 ปัจจุบันใช้ EGM2008 ที่มีความละเอียดถูกต้องมากกว่า ซึ่งโมเดลต่างๆเหล่านี้จะเป็นไฟล์ที่มีขนาดค่อนข้างใหญ่จนใหญ่มาก ทำให้ไม่สามารถนำมาใช้กับเครื่องคิดเลขในระดับนี้ได้ และถ้าเขียนโปรแกรมด้วย Casio Basic ก็คงหมดสิทธิ์ครับเพราะไม่มีคำสั่งในการอ่านเขียนไฟล์

ปัญหาและทางออก

ผมพยายามใช้โค้ดภาษาซีในการอ่านไฟล์โมเดล EGM96  ด้วย egm96-f477.c โปรแกรมสามารถคอมไพล์และบิวท์ผ่าน แต่ตอนประมวลผลอ่านไฟล์ เนื่องจากต้องใช้เมโมรีเครื่องคิดเลขมาอ่านไฟล์เข้าไปเก็บในหน่วยความจำพบว่า เมโมรีไม่พอ พยายาม optimize ด้วยตัวเองหลายๆทางแต่ไม่สำเร็จ จนมาพบโดยบังเอิญในเว็ปไซต์ ตามลิ๊งค์นี้ จุดประสงค์ทางผู้พัฒนาไลบรานี้เพื่อเขียนแอพลงบนโทรศัพท์มือถือแอนดรอยด์ วิธีการคือเขาได้ Normalize ค่า geoid separation จาก EGM2008 ขนาดเต็มลงมาเหลือขนาดแถวคอลัมน์ 181 x 91 ความถูกต้องจะถูกลดทอนลง ผมสนใจเลยเอามาคอมไพล์บน Casio SDK พบว่าทำงานได้ ใช้เมโมรีเครื่องคิดเลขที่มีอยู่เพียงพอ แต่ก็จวนเจียนจะใช้หมดเหมือนกัน เอาละความคิดที่เคยจะเขียนโปรแกรมคำนวณสเกลแฟคเตอร์บเครื่องคิดเลขก็มาเป็นจริง อาจจะได้ค่าสเกลแฟคเตอร์ที่ไม่ละเอียดเท่าคำนวณจากคอมพิวเตอร์แต่ในระดับทศนิยมหกตำแหน่ง ถือว่าได้ครับ

เครดิตผู้พัฒนาไลบรารี Geoid

ต้องขอบคุณผู้พัฒนาคือ Stefan Röttger (stefan@stereofx.org) มา ณ ที่นี้ ติดตามผลงานได้ที่เว็บไซต์ http://www.open-terrain.org/index.php ไลบรานี้ใช้สัญญาอนุญาตแบบ New BSD License

เตือนสติเรื่องสเกลแฟคเตอร์

ผมชอบวาทะของอ.ดีบุญ เมธากุลชาติ ก็ขออนุญาตเอาจากสไลด์มาแสดง ณ ที่นี้ด้วย อ.ดีบุญเรียก Grid Scale Factor ว่าตัวประกอบมาตราส่วน ส่วน Elevation Scale Factor เรียกว่าตัวประกอบความสูง ได้ใจความครบถ้วนดีครับ และก็เตือนสติคนใช้แผนที่ดังนี้

Grid Scale Factor เมื่อคูณกับ Elevation Scale Factor จะได้ Combined Scale Factor ตัวประกอบตัวนี้คือตัวสำคัญที่จะนำไปคูนกับระยะทางที่วัดบนพื้นผิวโลกเพื่อทอนลงมาบนระนาบราบของแผนที่ (Grid Distance)

ใช้ไลบรารีแปลงพิกัด MGRS

เนื่องจากโปรแกรมนี้ต้องแปลงค่าพิกัดจากระบบพิกัดฉาก  UTM  บนดาตั้ม WGS84 ไปยังค่าพิกัดภูมิศาสตร์ เพื่อส่งไปให้ฟังก์ชั่นคำนวณหาค่าสเกลแฟคเตอร์ที่ต้องการทราบค่าละติจูด ดังนั้นผมใช้ไลบรารี mrgs พัฒนาโดย Howard Butler ไลบรานี้ใช้สัญญาอนุญาตแบบ MIT License

สิ่งที่โปรแกรมต้องการ

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

  1. เลือกโซนของยูทีเอ็ม (UTM  Zone) ตั้งแต่โซน 1 ถึง 60 และเลือกซีกโลก (Hemisphere)  เลือกได้สองอย่างคือเหนือหรือใต้เส้นศูนย์สูตร
  2. ค่าพิกัดและค่าระดับของจุดที่ต้องการทราบค่าสเกลแฟคเตอร์ ต้องการค่าพิกัดเหนือใต้ในระบบยูทีเอ็มและค่าระดับที่อ้างอิงเหนือระดับน้ำทะเลปานกลาง (Orthometric Height above Meas sea level)

ในการคำนวณตัวประกอบความสูง (Elevation Scale Factor) โปรแกรมจะนำค่าความสูงเหนือระดับน้ำทะเลปานกลางมาคำนวณเป็นค่าระดับความสูงเทียบกับระดับทรงรี (Ellipsoid Height) ที่เราไม่สามารถนำค่าระดับความสูงเหนือระดับน้ำทะเลปานกลางมาใช้คำนวณเพราะว่า ตัวอย่างเข่นในประเทศไทย ค่าระดับน้ำทะเลปานกลางจะอยู่ใต้ทรงรีประมาณ 6-42 เมตร ถ้านำไปคำนวณจะได้ค่าตัวประกอบความสูงที่ไม่ถูกต้อง

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

ไปที่หน้าดาวน์โหลด Download มองหาโปรแกรม Scale Factor สำหรับเครื่องคิดเลข Casio fx-9860 G II SD ดาวน์โหลดไฟล์โปรแกรมชื่อ SFACTOR.G1A ติดตั้งลงเครื่องคิดเลข

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

เมื่อเมนูหลัก (Main Menu) เมื่อเลื่อนลูกศรมาที่ไอคอนโปรแกรมกดคีย์ “EXE” จะเข้ามาในโปรแกรมดังรูปด้านล่าง จะเห็นเมนูด้านล่างเรียงรายกันไป ได้แก่

Sett- Setting (ตั้งค่าโซนยูทีเอ็มและเลือกซีกโลก)

PSF – Point Scale Factor ( คำนวณหาค่าสเกลแฟคเตอร์ของจุดเดี่ยว

LSF – Line Scale Factor (คำนวณหาค่าสเกลแฟคเตอร์เฉลี่ยจากจุดสองจุด)

Info – Information (แสดงไลบรารีที่ผมใช้งานอยู่ ให้เครดิตแก่พัฒนา)

Exit – Exit Program (ออกจากโปรแกรม)

ตั้งค่ายูทีเอ็ม (UTM Settings)

ทีเมนูหลักกดคีย์ F1 – Sett จะเห็นหน้าตาจอภาพเครื่องคิดเลขดังนี้

เนื่องจากโปรแกรมบนเครื่องคิดเลขมีข้อจำกัดมากๆ ผมไม่สามารถใส่ระบบพิกัดได้มากเหมือน Surveyor Pocket Tools จึงจำกัดแค่ datum “WGS84” เท่านั้น สำหรับประเทศไทยก็เลือกโซนได้สองโซนคือ 47 และ 48 ตัวอย่างที่จะแสดงต่อไปเลือกเป็นโซน 47 ส่วนซีกโลก (Hemisphere) เลือกเป็น “North” กดคีย์ F6 – OK

 

คำนวณ Point Scale Factor

ตัวอย่างที่คำนวณต่อไปจะยกมาจากตัวอย่างที่แสดงในบล็อกของ Surveyor Pocket Tools เพื่อเปรียบ

เทียบผลลัพธ์กันได้ กลับมาที่เมนูหลักกดคีย์ F1-Sett ป้อนค่าพิกัด N=1499662.173, E=683068.285, H=0.791 

เสร็จแล้วกดคีย์ F1 – Calc เพื่อคำนวณสเกลแฟคเตอร์ของจุดนี้

ผลลัพธ์จะแสดงออกมาดังรูปด้านล่าง โปรแกรมจะแสดงค่าพิกัดและค่าระดับมาให้ จากนั้นแสดงค่าละติจูดและลองจิจูด แสดงความสูงเทียบกับทรงรี = -29.598 m กดคีย์ F2-PgDn เพื่อเลื่อนไปดูหน้าต่อไป

จะเห็นค่า Geoid Separation หรือ Geoid Height = -30.389 m และได้ค่า Elevation Scale Factor = 1.00000465, Grid Scale Factor = 1.00001458 สุดท้ายค่า Combine Scale Factor = 1.00001924 ลองเทียบกับค่าที่คำนวณด้วย Surveyor Pocket Tools ดังนี้

ค่า Combined Scale Factor  = 1.0000192644 ต่างกันที่ทศนิยมแปด ผมถือว่าใช้ได้เพราะเราใช้แค่ทศนิยมที่หกก็เพียงพอแล้ว ความต่างที่แตกต่างจากที่ผมเคยบอกไปแล้วคือมาจากค่า Ellipsoid height ที่คำนวณจากโปรแกรมเครื่องคิดเลขได้ -29.598 เมตร ส่วนค่าความสูงทรงรีที่ได้จาก Surveyor Pocket Tools ได้ค่า -29.776 เมตร ต่างกันประมาณ 20 ซม.

คำนวณหาค่า Line Scale Factor

จะแสดงไดอะแกรมของการคำนวณค่าสเกลแฟคเตอร์เฉลี่ยดังรูปด้านล่าง

ที่โปรแกรมเครื่องคิดเลขกลับมาที่เมนูหลัก

กดคีย์ F3-LSF จะเห็นหน้าตาโปรแกรมดังรูปด้านล่าง ป้อนค่าพิกัดจุดที่ 1 และกดคีย์ F2-Next เพื่อไปป้อนจุดที่ 2

พร้อมแล้วกดคีย์ F1-Calc ต่อไปจะได้ผลลัพธ์ดังนี้ เลื่อนลงไปดูหน้าต่อไปก็กดคีย์ F2-PgDn

สรุปแล้วจะได้ค่าเฉลี่ย Combined Scale Factor = 1.00082320

เปรียบเทียบกับค่าที่ได้จาก Surveyor Pocket Tools เท่ากับ 1.0008232302 เช่นเคยค่าต่างกันที่ทศนิยมแปด สามารถนำไปใช้งานได้

จัดเก็บข้อมูลและเรียกมาใช้ภายหลัง

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

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

จากนั้นนำ SDCard มาเสียบบนเครื่องคิดเลขอีกครั้ง เมื่อนำไปใช้งานได้สักพักถ้าเอามาเปิดอีกครั้งจะเห็นไฟล์หลายๆไฟล์ มีนามสกุลเป็น “CFG”  หมายถึง config ตัวอย่างถ้าใช้โปรแกรมคำนวณหาสเกลแฟคเตอร์นี้ไฟล์ที่จัดเก็บข้อมุลคือ “SFACTOR.CFG

เครดิตไลบรารี

ที่เมนูหลักกดคีย์ F5-Info จะเห็นเครดิตดังรูปด้านล่าง

โปรแกรมนี้มีขนาดใหญ่กว่าโปรแกรมอื่นๆก่อนเกือบสองเท่าเพราะว่าใช้ไลบรารีอื่นๆสามไลบรารีด้วยกัน ก็ติดตามชุดโปรแกรมเครื่องคิดเลข Casio fx-9860 G กันต่อไปครับ

ติดปีกเครื่องคิดเลขเทพ Casio fx 9860G II SD ด้วยโปรแกรมภาษาซีบน AddIn ตอนที่ 7 โปรแกรมคำนวณโค้งราบ (Horizontal Curve)

ติดปีกเครื่องคิดเลขเทพ Casio fx 9860G II SD ด้วยโปรแกรมภาษาซีบน AddIn ตอนที่ 7 โปรแกรมคำนวณโค้งราบ (Horizontal Curve)

โปรแกรมคำนวณโค้งราบ HCurve สำหรับเครื่องคิดเลข Casio fx-9860 G

ช่วงนี้อยู่ในชุดซีรี่ย์โปรแกรมเครื่องคิดเลข Casio fx-9860G ที่ใช้ไลบรารี MyLib ต่อไปขอนำเสนอโปรแกรมคำนวณโค้งราบ (Simple Horizontal Curve) ตัวโปรแกรมพัฒนาด้วยภาษาซี  ใช้เป็นโปรแกรม AddIn ดังรูปด้านล่าง

โปรแกรม HCurve คำนวณโค้งราบ (Simple Horizontal Curve)

องค์ประกอบของโค้งราบ (Elements of Horizontal Curve)

องค์ประกอบของโค้งราบ (Elements of Curve)

คำนิยาม (Abbrivations)

R – Radius คือรัศมีของโค้งราบ รัศมีจะตั้งฉากกับเส้นสัมผัสวงกลมเสมอ

PC –  Point of Curvature คือจุดเริ่มต้นโค้ง บางครั้งเรียกว่า BC (beginning of curve) หรือ TC (tangent to curve)

PI – Point of Intersection คือจุดตัดของ tangent 2 เส้น

PT – Point of Tangency คือจุดสิ้นสุดโค้ง บางครั้งเรียกว่า EC (end of curve) หรือ CT (curve to tangent)

POC – Point of Curve คือจุดบนโค้งในตำแหน่งใดก็ตาม

L – Length of Curve คือความยาวโค้งวัดตามโค้งจากจุด PC ไปจุด PT

T – Tangent Distance หรือ Tangent Length คือเส้นตรงที่สัมผัสโค้งวัดจากจุด PC ไปจุด PI หรือวัดจาก จุด PI ไปจุด PT

คุณสมบัติของโปรแกรม

โค้งราบ (Simple Horizontal Curve) ในงานสำรวจใช้ในงานสำรวจเพื่อการก่อสร้าง (Construction Survey) โดยที่ให้ตำแหน่ง (Setting out) งานก่อสร้างถนน ทางรถไฟ ผู้ใข้ป้อนข้อมูลองค์ประกอบของโค้งราบที่โปรแกรมต้องการจนครบ โปรแกรมสามารถคำนวณองค์ประกอบโค้งที่เหลือและสามารถคำนวณค่าพิกัดบนเส้น Center line หรือกระทั่งสามารถ offset ไปด้านซ้ายหรือด้านขวาก็ได้เช่นเดียวกัน สามารถคำนวณค่าพิกัดแบบทุกช่วงระยะ (interval) ให้ค่าพิกัดออกมาเป็นบัญชีรายการได้ โปรแกรมนี้ออกแบบและพัฒนามาเพื่อช่วยช่างสำรวจให้สามารถนำโปรแกรมไปตรวจสอบข้อมูลโค้งราบได้ด้วยตัวเอง

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

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

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

ไปที่หน้าดาวน์โหลด มองหาโปรแกรมคำนวณโค้งราบ HCurve  จะได้ไฟล์ชื่อ “HCURVEEX.G1A” เมื่อดาวน์โหลดมาแล้วโอนเข้าเครื่องคิดเลขผ่านทางโปรแกรม FA-124 หรือ SD Card

ใช้งานฟรี (Freely Usable)

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

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

ที่เครื่องคิดเลขกดคีย์ “Menu” เพื่อเข้าสู่หน้า AddIn หรือ  Main Menu ของเครื่องคิดเลข เลื่อนลูกศรลงไปด้านล่างๆจะเห็นไอคอนของโปรแกรมดังรูปด้านล่าง กดคีย์ “EXE” เพื่อเข้าสู่โปรแกรม

จะเห็นหน้าเมนูหลักของโปรแกรมดังนี้

ก็เหมือนโปรแกรมคำนวณโค้งดิ่ง (VCurve) ที่ผ่านมาคือมีเมนูย่อยเรียงจากซ้ายไปขวา สัมพันธ์กับการกดคีย์ F1-F6 เรียงกันไป

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

F1 (Know) – Known สำหรับป้อนจุดที่ทราบค่า station และค่าพิกัด

F2 (Angl) – Angle สำหรับป้อนที่เกี่ยวกับมุมเช่น back tangent azimuth และมุมเบี่ยงเบน (Deflection Angle) ตลอดจนถึงทิศทางของโค้ง (Curve Direction)

F3 (Elem) – Elements สำหรับป้อนองค์ประกอบของโค้งเช่นรัศมีหรือความยาวโค้าง

F4 (Info) – Information สำหรับคำนวณหาข้อมูลพื้นฐานของโค้งทั้งหมด

F5 (Calc) – Calculate สำหรับคำนวณหาค่าพิกัดโค้งได้หลายรูปแบบเช่นกำหนดสถานี ระยะ offset ตลอดจนคำนวณจากช่วงระยะทาง (interval) ที่กำหนดให้

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

ตัวอย่างการคำนวณโค้งราบ (Example)

มาดูวิธีการใช้งานจากตัวอย่างจะเข้าใจง่ายที่สุด

เลือกและป้อนสถานีที่ทราบและค่าพิกัด (Known Station and Know Coordinates)

ที่เมนูหลักกดคีย์ F1 (Know) แก้ไขค่าตามโจทย์ตัวอย่างที่ 1 ดังนี้

ในที่นี้สถานีและค่าพิกัดกำหนดที่จุด PI ทั้งคู่ เมื่อป้อนค่าเสร็จแล้วกดคีย์ F6 (OK) เพื่อออก

เลือกป้อนมุมอะซิมัทและมุมเบี่ยงเบน(Known Tangent and Deflection Angle)

ที่เมนูหลักกดคีย์ F2 (Angl) เพื่อเลือกอะซิมัทของเส้นสัมผัสที่ทราบค่าและมุมเบี่ยงเบน ตลอดจนป้อนทิศทางของโค้งราบว่าเลี้ยวซ้ายหรือเลี้ยวขวา จากตัวอย่างเลือกและป้อนค่าดังนี้

เสร็จแล้วกดคีย์ F6 (OK) เพื่อจัดเก็บค่าและออก

ป้อนองค์ประกอบโค้งราบ (Elements of Curve)

ที่เมนูหลักกดคีย์ F3 (Elem)  เลือกว่าจะป้อนค่ารัศมีโค้ง (Radius) หรือว่าความยาวโค้ง (Length of Curve) ในที่นี้เลือกรัศมีโค้งและป้อนค่า 201.950 เมตร

คำนวณหาข้อมูลพื้นฐานโค้งราบ (Curve Information)

ที่เมนูหลักกดคีย์ F4 (Info)

แสดงองค์ประกอบของโค้งและข้อมูลพื้นฐาน ค่าพิกัดของจุด PC, PI และ PT ตลอดจนจุดศูนย์กลางของโค้ง เนื่องจากจอภาพมีขนาดเล็กดังนั้นใช้การกดคีย์ F1 (PgUp) หรือ F2 (PaDn) เพื่อเลื่อนดูหน้าก่อนหน้านี้หรือหน้าถัดไป กดคีย์ F6 (Done) เพื่อออก

การคำนวณโค้งราบ (Horizontal Curve Calculation)

ที่เมนูหลักกดคีย์ F5 (Calc)จะเห็นเมนูย่อยอีกเมนูคือเมนูสำหรับคำนวณโค้งราบ

จะมีเมนูดังนี้

F1 (Sta) – Station คำนวณหาค่าพิกัดเมื่อกำหนดสถานี

F2 (INT) – Interval คำนวณหาค่าพิกัดสถานีเมื่อกำหนดช่วงระยะทาง (Interval) ให้

F4 (Info) – Information คำนวณข้อมูลพื้นฐานโค้งราบ โดยที่ผลลัพธ์เหมือนกับเมนู Info บนเมนูหลัก

F5 (Plot) – Plot Curve วาดรูปร่างโค้งราบ

F6 (Done) ออกจากเมนูคำนวณโค้งราบ

คำนวณหาค่าพิกัดเมื่อกำหนดสถานี (Calculate Coordinates of Station)

ที่เมนูคำนวณโค้งราบกดคีย์ F1 (Sta) จะมีไดอะล็อกให้ป้อนสถานี ตัวอย่างนี้ต้องการทราบค่าพิกัดของสถานี 17+200 โดยที่ offset ไปด้านซ้าย 8 เมตร ป้อนข้อมูลดังรูป ถ้าไม่ต้องการคำนวณหรือเก็บข้อมูลที่ป้อนก็กดคีย์ F5 (Canc) เพื่อ Cancel ออกไป หรือต้องการเก็บค่าแต่ไม่คำนวณก็กดคีย์ F6 (OK) ออกไป ถ้าต้องการคำนวณก็กดคีย์ F1 (Calc) จะได้ผลลัพธ์ดังรูปถัดไป กดคีย์ “EXE” เพื่อออก

คำนวณหาค่าพิกัดสถานีแบบกำหนดช่วงระยะทาง (Interval Calculation)

ที่เมนูคำนวณโค้งกดคีย์ F2 (INT) ในที่นี้ต้องการคำนวณทุกๆระยะ 25 เมตร โดยคำนวณในแนว Center Line (ไม่มีการ offset ไปซ้ายหรือขวา) กดคีย์ F1 (Calc) เพื่อคำนวณ

จะได้ผลลัพธ์ เริ่มตั้งแต่ PC (17+151.314), Sta 17+175, Sta 17+200, Sta 17+225, Sta 17+250, Sta 7+275, Sta 17+300 และสุดท้ายที่ PT (17+313.794) กด F1 (PgUp) เพื่อเลื่อนไปหน้าก่อนหน้านี้ และกด F2 (PgDn) เพื่อไปดูหน้าถัดไป กด F6 (Done)

 

วาดรูปโค้งราบ (Plot Curve)

จากเมนูคำนวณโค้งราบ กดคีย์ F5 (Plot) จะมีเมนูย่อยลงไปอีกสำหรับย่อ F1 (Z-) ขยาย F2 (Z+) ดึงรูปไปด้านซ้าย F3 (Lt) ดึงรูปไปด้านขวา F4 (Rt) กดคีย์ F5 (>) เพื่อไปเมนูย่อยอีกเมนูด้านขวา เมนูด้านขวาจะมีดึงรูปลง F1 (Dn) หรือดึงรูปขึ้น F2 (Up) กดคีย์ F6 (Exit) เพื่อออกจากเมนู

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

จัดเก็บข้อมูลและเรียกมาใช้ภายหลัง

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

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

จากนั้นนำ SDCard มาเสียบบนเครื่องคิดเลขอีกครั้ง เมื่อนำไปใช้งานได้สักพักถ้าเอามาเปิดอีกครั้งจะเห็นไฟล์หลายๆไฟล์ มีนามสกุลเป็น “CFG”  หมายถึง config ตัวอย่างถ้าใช้โปรแกรมคำนวณโค้งราบนี้ไฟล์ที่จัดเก็บข้อมุลคือ “HCVEX.CFG