Tag: ellipsoid

ติดปีกเครื่องคิดเลขเทพ 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 มองหาโปรแกรมสำหรับเครื่องคิดเลข 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 เช่นเคยค่าต่างกันที่ทศนิยมแปด สามารถนำไปใช้งานได้

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

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

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

สนุกกับโปรแกรมเครื่องคิดเลขสำหรับงานสำรวจ ตอนที่ 3 โปรแกรมคำนวณระยะทางบนทรงรี สำหรับเครื่องคิดเลข Casio FX 5800P

โปรแกรมคำนวณระยะทางที่สั้นที่สุดบนทรงรี (Geodesic Distance)

  •  สวัสดีครับผู้อ่านกลับมาพบกันอีกครั้ง ครั้งนี้หอบเอาสูตรยาวๆมาฝากกัน เรื่องระยะทางบนทรงรีความจริงจัดอยู่ในหมวด  Geodesy ที่ถือว่าเป็นเรื่องยากสำหรับนักศึกษาสมัยก่อน เพราะคำนวณทีต้องเปิดตารางล็อก สมัยนี้ถ้าทำความเข้าใจก็ไม่ได้ยากแล้วครับ มีตัวช่วยมากมายเช่นเครื่องคอมพิวเตอร์
  • ระยะทางที่สั้นที่สุดบนทรงรีเรียกว่า Geodesic Distance ถือว่าเป็นะระยะทางที่สั้นที่สุดบนทรงรี จัดเป็นสูตรที่มีมานานนมแล้ว นักคณิตศาสตร์สมัยก่อนคิดค้นขึ้นมาเพื่อช่วยในการคำนวณระยะทางที่สั้นที่สุดเช่นสำหรับการเดินเรือเดินข้ามมหาสมุทร เพราะแผนที่ที่เราใช้กัน ถ้าเมืองท่าที่ต้องการเดินทางนั้นอยู่กันไกลหลายพันไมล์ทะเล จะเอาดินสอมาขีดตรงๆเชื่อมกันบนแผนที่เดินเรือ จะได้ระยะทางที่ไม่ใช่ใกล้ที่สุด จะกลายอ้อมไป
  • เคยสังเกตเวลานั่งเครื่องบินไปต่างประเทศไกลๆไหมครับ จะมีรูปเส้นทางบินบนจอมอนิเตอร์ให้ดู จะเห็นเส้นทางบินเป็นเส้นโค้งๆ ไม่ใช่ตรงๆ นั่นเป็นเพราะเส้นโครงแผนที่ ความจริงแนวเครื่องบินก็ไม่ได้เฉไปไหน แต่เมื่อเอาแนวบินมาวาดบนแผนที่ที่ใช้เส้นโครงแผนที่แบบรักษารูปร่าง(Conformal) ที่เราใช้กันในปัจจุบันเช่น UTM เส้นทางบินจะกลายเป็นเส้นโค้งไป

  • ก่อนจะไปต่อ ผมอยากจะพูดถึงระยะทางอีกระยะทางหนึ่งบนทรงกลม เรียกว่า Great Circle Distance คือระยะทางที่สั้นที่สุดบนทรงกลม คนโบราณนำระยะทางนี้มาคำนวณเส้นทางการเดินเรือ เส้นทางการบิน แต่เนื่องจากโลกเราไมได้กลมแต่ทรงคล้ายออกมาทางทรงรีมากกว่า Geodesic distance เลยแม่นกว่า Great circle distance แต่คนก็ยังนิยมใช้ Greate circle distance มากกว่าอยู่ดีถึงจะคลาดเคลื่อนจาก Geodesice distance ไปประมาณ 0.1% แต่ก็คำนวณง่ายกว่า ในสูตรไม่มีการวนลูป เรียกสูตรการคำนวณบนทรงกลมนี้ว่า “Haversine
  • ถ้าสนใจเรื่อง Geodesic distance และ Greate circle distance ในรายละเอียดและโปรแกรมคำนวณบนคอมพิวเตอร์ติดตามได้ในโปรแกรม Survey Pocket Tools ของผมได้ที่บล็อก priabroy.name

การประยุกต์ใช้งาน

  • บางคนอาจจะว่าไกลตัว แต่เวลาเราค้นหาเส้นทางสำหรับโปรแกรมนำทาง (Navigator) บนโทรศัพท์มือถือทั้งหลาย ทราบไหมว่าเขาคำนวณระยะทางได้อย่างไร อย่างของ google maps ไม่มีเอกสารเปิดเผย แต่หลักการคือเส้นทางต่างๆของกูเกิ้ลจะต่อกันเป็นจุดๆ เรียกว่า node แต่ละโหนดจะมีค่าพิกัดแลตติจูด/ลองจิจูด กำกับอยู่ โปรแกรมจะคำนวณหาระยะทางระหว่างแต่ละโหนดด้วยการนำค่าพิกัดดังกล่าวมาเข้าสูตรคำนวณระยะทางที่สั้นที่สุดบนทรงรี (geodesic distance) หรือไม่ก็ระยะทางบนทรงกลม (great circle distance) เมื่อนำระยะทางแต่ละเส้นมารวมๆกันก็จะได้เป็นระยะทางแสดงให้เราดูบนหน้าจอโทรศัพท์มือถือ
  • ส่วนการเลือกเส้นทางจากต้นทางไปจุดหมายปลายทางแต่ละสถานที่ ว่าไปตามถนนเส้นไหนจะสั้นที่สุดอยูในเรื่อง Travelling Salesman Problem (TSP) อันนี้ลึกซึ้งมากครับ ไม่ขอกล่าวถึง
  • ในฐานะช่างสำรวจ เราใช้กันมันเกือบทุกวันโดยที่ไม่รู้ตัว ก็ลองมาศึกษากันหน่อยว่ามันทำงาน คำนวณมาให้ได้อย่างไร แต่ก่อนจะไปกันต่อ มาลองซักซ้อมเรื่องความรู้เซอร์เวย์กันเล็กน้อยก่อน

เกร็ดความรู้เล็กน้อยสำหรับช่างสำรวจ

  • วันนี้มาว่ากันต่อเรื่องมุมอะซิมัทและฟังก์ชั่น Rec() บนเครื่องคิดเลข Casio FX-5800P ฟังก์ชัน Rec() เอาไว้คำนวณหาค่าพิกัดปลายทาง เมื่อทราบค่าพิกัดต้นทาง, ระยะทางและอะซิมัท
  • Rec(r,θ) เขียนใหม่เป็น Rec(ระยะทาง, อะซิมัท) ก่อนจะเข้าวิธีการใช้งาน มาเท้าความกันหน่อย

  • ลองมาดูรูปด้านบน กำหนดมุมอะซิมัท (θ) และระยะทาง (R) ต้องการหาระยะทางไปตามแกน X จากจุด A ไปหาจุด B หา ได้เท่ากับ  ระยะทาง x Sine (มุมอะซิมัท)
    • ΔX = R x Sin(θ)  เรียกว่า “Departure” คุ้นๆไหม
  • ระยะทางไปตามแกน Y จากจุด A ไปจุด B
    • ΔY= R x Cos(θ)  เรียกว่า “Latitude”
  • ดังนั้นค่าพิกัดจุด A หาได้จาก
    • XB  = XA + ΔX = XA + R x Sin(θ)
    • YB  = YA + ΔY = YA + R x Cos(θ) 
  • ถ้าจำสูตรยากหรือขี้เกียจจำ ฟังก์ชั่น Rec() ช่วยได้ มาดูตัวอย่าง

  • ต้องการทราบค่าพิกัด A เมื่อทราบค่า (X,Y) ของจุด B ดังรูปด้านบน ระยะทาง 100 เมตร อะซิมัทจาก B ไป A = 240º
    • แทนค่าลงไปในสูตร Rec(ระยะทาง, อะซิมัท) = > กดเครื่องคิดเลขเลย Rec(100,240) ผลการคำนวณได้ผลลัพธ์มาสองอย่างคือ ΔX เก็บไว้ในตัวแปร “J” ส่วนค่า ΔY เก็บไว้ในตัวแปร “I”
    •  ผู้อ่านจะสังเกตว่าทำไม  ΔX ไม่เก็บไว้ในตัวแปร “I” และ ΔY ไม่เเก็บไว้ในตัวแปร “J” อย่างที่มันควรจะเป็น ต้องไม่ลืมว่ามุมที่เขาเขียนฟังก์ชันตัวนี้มาคือมุมกวาดจากแกน X แต่เราใช้มุมกวาดจากแกน Y มันเลยสลับร่างสร้างรัก ด้วยประการฉะนี้
  • ลอง Recall (RCL) ค่ามาดูจะได้ I = -50 และ J = -86.603
    • XA  = XB + ΔX = 586.603 – 86.603 = 500
    • YA  = YB + ΔY =  550 – 50 = 500
    • คำตอบคือค่าพิกัดจุด A (500,500)
  • มุมอะซิมัทเป็นมุมที่มหัศจรรย์ครับ เวลาคูณ Cos กับ Sin มันจะไปตกควอดแรนท์ที่ถูกต้องให้ โดยที่หัวสมองเราไม่ต้องไปนึกตาม ลองคิดถึงุมุม Bearing สิครับ สยดสยอง แต่โชคดีบ้านเราไม่ได้ใช้
  • ยังมีเกร็ดเล็กเกร็ดน้อยเกี่ยวกับอะซิมัทอีกนิด มาว่ากันต่อตอนต่อไป

โปรแกรมคำนวณระยะทางที่ส้ั้นที่สุดบนทรงรี (Geodesic Distance)

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

  • การคำนวณจะมีการวนลูปด้วยคำสั่ง For … Next ทุกคร้้งจะมีการเปรียบเทียบค่าใหม่กับค่าเริ่มต้นถ้าค่าต่างกันน้อยกว่า สิบยกกำลังลบสิบสอง (10-12) การคำนวณจะหยุดทันทีแสดงว่าค่าใช้ได้ ผมลองจับจำนวนครั้งที่วนลูปดูประมาณ 3-4 ครั้งแค่นั้น ในตัวโปรแกรมตัวแปรไปดึงตัวแปรอนุกรม Z มาใช้ทั้งหมด 10 ค่า Z[1], Z[2], Z[3], …. Z[10] ตอนประกาศใช้ใช้คำสั่ง 10→DimZ ตอนเลิกใช้ตอนท้ายโปรแกรมประกาศ 0→DimZ เพื่อคืนเมมโมรีให้กับเครื่องคิดเลข
  • เวลาคีย์ชื่อโปรแกรมเนื่องจำกเครื่องคิดเลขจำกัดความยาวตั้งชื่อสั้นว่า “GEODESIC

วิธีการใช้งานโปรแกรม

  • ต้องการระยะทาง Geodesic distance
    • จากจุดที่ 1 ค่า latitude = 14°27’27.0″N longitude = 100°54’12.57″E
    • ไปหาจุดที่ 2 ค่า latitude = 14°36’49.53″N longitude = 98°1’39.63″E
  • เรียกโปรแกรมก่อน กดที่ “File” เลื่อนไปที่ “Geodesic

  • ป้อนค่าพิกัดจุดที่ 1 ค่า latitude = 14°27’27.0″N longitude = 100°54’12.57″E

  • ป้อนค่าพิกัดจุดที่ 2 ค่า latitude = 14°36’49.53″N longitude = 98°1’39.63″E

ผลลัพธ์ของการคำนวณ

  • ผู้อ่านสังเกตไหมเวลากด “EXE” เพื่อคำนวณขั้นสุดท้าย เครื่องจะวูบไปหลายวินาที แต่ถ้าโปรแกรมนี้ถ้าย้ายไปเครื่อง Casio FX-9860G II คงใช้เวลาน้อยกว่าเพราะตามสเป็คแล้วแรงกว่า FX-5800P ประมาณสองเท่า มาดูผลลัพธ์กัน ได้ระยะทาง = 310432.516 เมตร หรือ 310.433 กม.

  • เทียบกับโปรแกรม Geodesic ใน Surveyor Pocket Tools ได้ค่าเท่ากัน ซึ่งหมายเหตุอีกนิดว่าอัลกอริทึ่มที่ใช้คำนวณในโปรแกรม Surveyor Pocket Tools ผมใช้ไลบรารีจาก GeographicLib ซึ่งผลลัพธ์การคำนวณตรงกับสูตรที่เราป้อนในเครื่องคิดเลข ก็ทำให้เรามั่นใจครับว่ามาถูกทาง สูตรถูกต้องแล้ว

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

Surveyor Pocket Tools – คำนวณพื้นที่ เรื่องธรรมดาที่ไม่ธรรมดา (ตอนที่ 1)

  • Surveyor Pocket Tools เป็นโปรแกรมที่รวบรวมเครื่องมือเล็กเครื่องมือน้อยสำหรับช่างสำรวจ ตอนนี้เพิ่มการคำนวณหาพื้นที่ หลากหลายประเด็นที่จะมาคุยกันว่ามันควรจะง่ายธรรมดาแต่มันไม่ธรรมดาอย่างไร
โปรแกรมคำนวณพื้นที่
โปรแกรมคำนวณพื้นที่

ย้อนรอยสูตรคำนวณหาพื้นที่

  • Surveyor Pocket Tools รุ่นก่อนหน้านี้ส่วนใหญ่จะเป็นงานคำนวณในงานเซอร์เวย์เรื่อง Advance เช่นงานคำนวณหาระยะทางบน Ellipsoid งานคำนวณหาความสูงของจีออยด์ การแปลงพิกัด ตอนนี้จะกลับมาเรื่องพื้นฐานคือการคำนวณหาพื้นที่
  • การคำนวณหาพื้นที่ ปกติที่เราคุ้นหูคุ้นตาคือ การเอาค่าพิกัด X,Y ในระบบพิกัดฉากมาวางเรียงกันในแนวนอนหรือแนวตั้งแล้วจัดการเอาค่า X คูนค่า Y สลับแบบผูกเชือกรองเท้า คูนขึ้นเป็นบวก คูนลงเป็นลบ เรียกสูตรนี้ว่า  สูตรเชือกผูกรองเท้า (Shoelace’s formula) แต่เนื่องจากช่างสำรวจคงใช้สูตรนี้กันมากกว่าเพื่อนเลยเรียกอีกชื่อหนึ่งว่า Surveyor’s formula
soffice-bin_2016-12-22_11-12-12
Shoelace Formula (รูปดัดแปลงจากรูปต้นฉบับ wikipedia.org)
  • โดยที่ข้อแม้สูตรนี้คือจุดจะต้องเรียงกันไปไม่มีการไขว้กันเด็ดขาดแบบรูปด้านล่างsoffice-bin_2016-12-22_11-34-25
  • งานสำรวจแบบโลกเก่าคือออกพิกัดลอยของใครของมัน ใช้กล้อง Total station หรือโบราณกว่านั้นดึงเทปวัดมุม เก็บมุมพื้นที่ของแปลงแล้วนำค่าพิกัดรูปปิดมาคำนวณ ข้อดีคือได้พื้นที่จริงๆมาใช้งาน แต่ข้อเสียคือแผนที่ต่างคนต่างออกจากศูนย์พิกัดต่างกัน ถ้าแปลงติดๆกันไม่สามารถเอามาต่อกันได้ เอาไปเขียนลง google earth ก็ไม่ได้ ไม่รู้ว่ามันอยู่ตรงไหนของโลกนี้

ปัญหาการคำนวณพื้นที่กับค่า Scale factor

  • ในการทำแผนที่พื้นที่ใหญ่ๆ เนื่องพื้นผิวโลกมีความโค้ง ดังนั้นการทำแผนที่บนผิวโลกแล้วนำมาเขียนลงบนระบบแบนราบแบบกระดาษ จะเกิดความเพี้ยน (distortion) ยิ่งไกลมากยิ่งเพี้ยนมาก นักคณิตศาสตร์รุ่นเก่าจึงแก้ปัญหานี้โดยการคิดค้นประดิษฐ์เส้นโครงแผนที่ (Map projection) ที่เป็นระนาบราบ (ระบบพิกัดฉาก) โดยมีการฉาย (project) จุดจากบนผิวโลกลงมาจุดบนระนาบ
  • มีได้ย่อมมีเสีย ข้อดีคือแผนที่ต่างคนต่างทำ แต่อยู่บนระบบพิกัดฉากเดียวกัน ศูนย์กำเนิดเดียวกันสามารถนำมาต่อกันได้เพราะอยู่ในระบบเดียวกัน ข้อเสียคือแต่ละจุดบนเส้นโครงแผนทีจะมี Scale factor ที่แตกต่างกันไปไม่เท่ากัน การทำแผนที่บนเส้นพิกัดฉาก ในทางปฏิบัติจะต้องมีการทอนระยะทางที่วัดได้บนพื้นโลกด้วยการคูนเข้าค่าคูนมาตราส่วน (scale factor) จึงจะได้ระยะทางบนระบบพิกัดฉาก ปัญหาที่ scale factor แต่ละที่ไม่เท่ากันจึงเป็นปัญหาคลาสสิคที่สร้างความเวียนหัวให้กับคนทำแผนที่พอสมควร
  • ดังนั้นการทำแผนที่บนระบบพิกัดฉากเช่น UTM ทุกสิ่งทุกอย่างที่วัดบนพื้นโลกจะถูกทอนลงมาบนระบบพิกัดฉากด้วยค่า scale factor ปัญหาก็คือ ถ้านำค่าพิกัดของรูปปิด (polygon) ระบบพิกัดฉากมาคำนวณจะได้พื้นที่เรียกว่า Grid-based area จะไม่ตรงกับพื้นที่จริงๆ (Surface area) พื้นที่จริงๆนี่แหละครับที่เราต้องการ ที่จะนำไปใช้ในการซื้อขายที่ดิน หรือนำมาบริหารจัดการเกี่ยวกับที่ดินเช่นด้านการเกษตร เช่นการเตรียมปุ๋ยให้สอดคล้องกับพื้นที่จริงๆในกรณีมีที่ดินมากๆ

ทำความเข้าใจเรื่อง Scale Factor

  • ก่อนจะไปต่อเรื่องคำนวณพื้นที่จะมาทำความเข้าใจเรื่อง Scale factor บนแผนที่ระบบพิกัดฉาก Scale factor จะมีผลต่อความยาว ระยะทาง ทำให้มีผลต่อการคำนวณพื้นที่โดยปริยาย
Scale factor
Scale factor
  • จากรูปโมเดลด้านบน จะมีพื้นผิวสี่อย่างคือ
    • พื้นผิวภูมิประเทศบนโลก ช่างสำรวจทำการรังวัดบนนี้ครับ
    • พื้นผิวจีออยด์ (Geoid) คือพื้นผิวที่แทนระดับน้ำทะเลปานกลาง เราวัดความสูงสิ่งต่างๆบนโลกอ้างอิงค่าระดับจากพื้นผิวนี้
    • พื้นผิวทรงรีอาจจะเป็น WGS84 (ช่างรังวัดกรมที่ดินใช้ทรงรี Everest 1830) ค่าพิกัดภูมิศาสตร์ที่ได้จากเครื่อง GPS อ้างอิงที่พื้นผิว WGS84 นี้
    • พื้นระนาบราบ(ระบบพิกัดฉาก) บ้านเราใช้ UTM แผนที่ที่เราใช้กันส่วนใหญ่อยู่บนระบบนี้
  • การคำนวณหาค่าคูนมาตราส่วน (Scale factor) จะมีสองอย่างประกอบกัน โดยที่คำนวนสองขั้นตอน
    • Elevation Scale Factor (ESF) เป็นค่าคูนมาตราส่วนที่คูนกับระยะราบที่รังวัดได้บนพื้นผิวภูมิประเทศของโลกจะได้ระยะราบบนทรงรี หรือ ระยะราบบนทรงรี = ESF x ระยะราบบนพื้นโลก
    • Grid Scale Factor (GSF)  เป็นค่าคูนมาตราส่วนที่คูนกับระยะราบบนทรงรีจะได้ระยะราบบนระบบพิกัดฉาก หรือ ระยะราบบนระบบพิกัดฉาก = GSF x ระยะราบบนทรงรี
  • ดังนั้นถ้ารวบรัดสองขั้นตอนในขั้นเดียวจะได้ ระยะราบบนระบบพิกัดฉาก = ESF x GSF x ระยะราบบนพื้นโลก ค่า ESF x GSF นิยมเรียกว่า Combined Scale Factor (CSF)

ขั้นตอนที่ 1 คำนวณ Elevation Scale Factor (ESF)

  • ชื่อก็บอกมาแล้วว่าเกี่ยวข้องกับค่าระดับ ค่าระดับที่ใช้เป็นความสูงเหนือทรงรี นะครับไม่ใช่ระดับน้ำทะเลปานกลางเหนือจีออยด์ ถ้ากำหนดค่าระดับจากน้ำทะเลปานกลาง (MSL) ก็ใช้โปรแกรม EGM ของผมแปลงลงมาได้จากสูตร h = H + N (ความสูงจากระดับน้ำทะเลปานกลาง + ความสูงจีออยด์)
  • มาดูสูตรกันก่อน R คือรัศมีทรงรี ถ้าใช้ WGS84 ใช้ค่าหยาบประมาณ 6,372,000 เมตร

esf

  • สมมติจุดที่ต้องการหา ESF อยู่แลตติจูด = 12.9725° ลองจิจูด = 99.8329° ค่าระดับ =63.9 เมตร(เหนือระดับน้ำทะเลปานกลาง) คำนวณหา N ได้ด้วยโปรแกรม EGM (อยู่ในชุด Surveyor Pocket Tools อยู่แล้ว) ได้ความสูงเหนือจีออยด์ (N) = -31.6638 เมตร

surveyor-pocket-tools_2016-12-23_11-19-39

  • ความสูงเหนือทรงรี (h) = H + N = 63.9 – 31.6638 = 32.236 เมตร
  • ค่า ESF แบบหยาบๆ ใช้ค่า R โดยประมาณ ESF = 6372000/(6372000+32.236) = 0.999994941
  • ถ้าต้องการค่า ESF ทีละเอียดระดับทศนิยมตัวที่ 8-9 จะต้องคำนวณหาค่า R ละเอียดจากค่าแลตติจูดของบริเวณนั้นก่อน

soffice-bin_2016-12-24_09-21-50

soffice.bin_2017-02-21_11-45-39

 

  • ถ้าทรงรี WGS84 a = 6378137, b = 6356752.314, f=1/298.2572236, e² = 0.00669437998 แลตติจูด (Ɵ) = 12.9725° แทนที่ลงไปสูตร
  • R = 6378137 x √(1-0.00669437998 ) / (1-0.00669437998 x (sin(12.9725))² = 6358897.478
  • ESF = 6358897.478 / (6358897.478 + 32.236) = 0.999994931 (ค่าคำนวณแบบหยาบ 0.999994941)
  • สรุปได้ค่า ESF = 0.999994931 (ชีวิตยุ่งยากไปไหม โอกาสหน้าผมจะทำทูลส์เล็กๆ ไว้คำนวณหา ESF)

ขั้นตอนที่ 2 คำนวณหา Grid Scale Factor (GSF)

  • ต่อไปจะคำนวณหา GSF ค่านี้จะไม่ขึ้นกับค่าระดับ แต่จะขึ้นกับค่าลองจิจูด เช่นเดียวกันจะมีสูตรคำนวณค่าโดยประมาณและค่าละเอียด สูตรค่าละเอียดค่อนข้างจะยาวหน่อย ด้านล่างเป็นสูตรหาค่าโดยประมาณ

soffice-bin_2016-12-24_09-07-23

  • เช่นถ้าจุดมีค่าพิกัด N = 1,434242.632 E = 590,334.084. UTM Zone 47N WGS84 หา GSF โดยประมาณได้โดยการนำตัวเลขไปแทนที่สูตรด้านบน ใช้ค่ารัศมีโลก 6372000 เมตร จะได้ GSF = 0.9996(1 + (590334.084 – 500000)²/(2 * 6372000)) = 0.9997004496
  • ต่อไปมาดูค่าคูนมาตราส่วน Grid Scale Factor แบบละเอียด ก่อนจะไปต่อ ต้องรู้พารามิเตอร์ทรงรีที่ใช้ รู้ค่าพิกัดภูมิศาสตร์ในตำแหน่งที่ต้องการหา

soffice-bin_2016-12-24_09-21-50

  • ค่าที่ต้องเตรียมก่อนจะคำนวณหา Grid Scale Factor (GSF)

soffice-bin_2016-12-26_07-52-18

 

  • สูตร Grid Scale Factor สำหรับเส้นโครงแผนที่ UTM
soffice-bin_2016-12-24_11-09-00
Grid Scale Factor for UTM
  • มาลองคำนวณ กำหนดจุดมีค่าพิกัดดังนี้ N = 1,434242.632 E = 590,334.084 อยูบน UTM Zone 47N WGS84 เราต้องการหา latitude/longitude เปิดโปรแกรม “UTM -Geo converter” แล้วป้อนค่าพิกัด คำนวณดังรูป ได้ค่า latitude(Φ) = 12.9725056694 longitude(λ) = 99.8329042500 central meridian สำหรับโซน 47 คือ λ0=99  สำหรับนำไปคำนวณต่อ

surveyor-pocket-tools_2016-12-24_09-36-29

  • แทนที่สูตรหาค่า T = (tan(12.9725056694))² = 0.053067021
  • แทนที่สูตรหาพารามิเตอร์ของทรงรี WGS84 a = 6378137 f = 1/298.2572236 หาค่า e = √(2f – f²) = 0.08181919084 และค่า e’ = √(e²/(1-e²)) = 0.08209443794
  • หาค่า C = 0.08209443794² x (cos(12.9725056694))² = 0.00639987446
  • หาค่า A = (99.8329042500  – 99) x 3.141592654/180 * (cos (12.9725056694))² = 0.01416590874 (note : แปลงเป็น ดีกรี เป็น radian คูน ¶/180.0 เข้าไป)
  • หาค่า GSF (k) = 0.9996 x 1.000100987 = 0.9997009464 (ค่าโดยประมาณ 0.999704496)
  • สรุปค่า  GSF = 0.9997009464 

คำนวณ Combined Scale Factor (CSF)

  • ตัวอย่างข้างต้นที่ผมคำนวณหาค่า ESF กับ GSF เป็นจุดเดียวกัน ดังนั้น CSF = ESF x GSF
  • CSF = 0.999994931 x 0.9997009464 = 0.9996958785
  • ผมมีหมุด GPS อีกตัวหนึ่ง ค่าพิกัด  N = 1434091.770 E=590362.138   ค่าระดับ(รทก.) =  60.513 m.
  • คำนวณระยะราบได้ 153.448 เมตร ลองทอนระยะราบบนกริดกลับไปบนพื้นภูมิประเทศบนโลก = 153.448 / 0.9996958785 = 153.495 เมตร ถ้าเอากล้อง Total station ลองวัดระยะราบดูควรจะต่างจากนี้ไม่มาก ที่จริง combined scale factor จะต้องหาอีกจุดด้วยก่อนนำมาเฉลี่ยกัน แต่คงต่างกันไม่มากนัก ผมแค่อยากจะบอกว่า ระยะราบของหมุด  GPS สองตัว ที่เป็นระยะราบบนพิกัดฉาก UTM ทำไมเวลาช่างสำรวจไปตั้งกล้อง Total station แล้ววัดระยะราบแล้วทำไมไม่เท่ากัน ไม่เท่ากันเพราะอะไร แล้วจะคำนวณทอนระยะทางไปหากันได้อย่างไร

การประยุกต์ใช้ Combined Scale Factor

  • หลายๆท่านที่เป็นวิศวกรสำรวจ ช่างสำรวจ อาจจะผ่านงานโครงการใหญ่ๆ ยาวๆ เช่นงานถนน งานรถไฟ อาจจะมีโอกาสได้ใช้ CSF ปัจจุบันกล้อง total station ทันสมัย เราสามารถป้อนเอาค่า CSF พวกนี้เข้าไปในกล้องได้ กล้องเกือบๆทุกยี่ห้อจะเรียกตัวนี้ว่า “scale factor” ทำให้ลดความยุ่งยากกว่าสมัยแต่ก่อนที่ดึงเทปวัดระยะ ที่ต้องมาทอนระยะสดๆ เวลาไป setting out ก็อัพโหลดค่าพิกัดก่อน ในสนามก็สามารถ layout out ได้ถูกต้องไม่ต้องกังวลเพราะเรื่องที่เหลือกล้องจัดการให้
  • โครงการยาวๆแบบนี้ จึงต้องมาแบ่งตัดพื้นที่ออกเป็น block แต่ละ block คำนวณหาค่าเฉลี่ย CSF ที่แทนแต่ละพื้นที่ เวลาก่อสร้างช่วงรอยต่ออาจจะมีปัญหาเล็กน้อยที่โครงสร้างอาจจะไม่ต่อกันเป๊ะทีเดียว
  • ข้อควรจำ งานก่อสร้างประเภทนี้ จะออกแบบไว้บนระบบพิกัดฉาก UTM ดังนั้นคนออกแบบควรจะระลึกไว้ตลอดว่ากำลังออกแบบบนกริด ดังนั้นของที่วัดได้บนแบบ 500 เมตร ไม่ใช่ไปวัดในสนามได้ 500 เมตร (ยังมีคนที่เข้าใจแบบนี้อยู่ไม่น้อย) ผมเคยอ่านเจอบางกรณีที่ต่างประเทศ ถ้าบริเวณนั้น  CSF = 0.9995 คนออกแบบจะลง dimension ไว้ 499.75 เมตร ในแบบ เพื่อให้ที่เวลาก่อสร้างจริงจะวัดได้พอดี 500 เมตร
  • มีเอกสารที่เกี่ยวข้องกับการประยุกต์ใช้ Scale factor มากมายบนอินเทอร์เน็ตให้อ่านศึกษา ผมจะลงไว้ท้าย blog
  • เพราะระยะทางมันมี scale factor ที่เปลี่ยนผันไปตามพื้นที่และเป็นระบบ ดังนั้นเรื่องการคำนวณหาพื้นที่จึงต้องคำนึงถึงเรื่องนี้

การใช้งานโปรแกรมคำนวณพื้นที่ (Area) ในชุดโปรแกรม Surveyor Pocket Tools

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

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

  • ดูด้านขวามือดูตรง “ดาวน์โหลด (Download” มองหา “Surveyor Pocket Tools Vxxx build xxxx” คลิกเพื่อดาวน์โหลด สำหรับวินโดส์ 32 บิต ผมยังทำไฟล์ setup ให้อยู่ เมื่อโหลดมาแล้วก็ unzip แล้วก็ทำการติดตั้งได้ง่ายๆ
  • เมื่อติดตั้งแล้วคลิกที่โปรแกรม Surveyor Pocket Tools จะเห็นหน้าต่าง คลิกที่ไอคอนรูปโพลีกอน “Area” จะเห็นโปรแกรมคำนวณหาพื้นที่ขึ้นมาดังนี้
"Compute Area" on Surveyor Pocket Tools
“Compute Area” on Surveyor Pocket Tools

ส่วนประกอบของโปรแกรม

ส่วนประกอบของโปรแกรม
ส่วนประกอบของโปรแกรม

ขั้นตอนใช้งานพอสังเขป

  • การใช้งาน ผมจะกล่าวย่อๆเป็น work flow ดังนี้
    • เริ่มต้นจากเตรียมไฟล์ CSV ที่เก็บค่าพิกัดของรูปปิดที่จะหาพื้นที่ ค่าพิกัดอาจจะเป็นค่าพิกัดระบบพิกัดฉากบน UTM (อาจจะบน WGS84 หรือ Indian 1975) หรือว่าค่าพิกัดภูมิศาสตร์บน WGS84 ก็ได้
    • นำเข้าด้วยการ import โปรแกรมจะอ่านไฟล์แล้วทำการประมวลผลว่ามีตัวคั่นที่เหมือนกันทั้งไฟล์ไหม หรือว่ามีบรรทัดที่ไม่เหมือนบรรทัดอื่น ถ้าอ่านไม่ได้จะแจ้งให้ผู้ใช้ทราบ เมื่ออ่านเข้ามาแล้วจะแสดงผลที่ตารางข้อมูล
    • ทำการกำหนด header หรือหัวคอลัมน์ว่าช่องไหนเป็น Easting/Northing หรือ Latitude/Longitue หรือกำหนดชื่อจุดด้วยถ้ามี แต่ถ้าโปรแกรมมีหัว (header) โปรแกรมพยายามจัดให้อัตโนมัติ
    • ตั้งรูปแบบมุม ถ้าค่าพิกัดทีอ่านมาเป็นค่าพิกัดภูมิศาสตร์
    • คำนวณพื้นที่
    • ปักหมุดบน Google Maps หรือว่า Google Earth
    • ดูผลการคำนวณพื้นที่แบ่งเป็นพื้นที่บนทรงรี และพื้นที่บนระบบพิกัดฉาก
    • แปลงหน่วยพื้นที่เช่นจากตร.ม. เป็นไม่
    • Export ข้อมูลผลลัพธ์ เป็นไฟล์ CSV, Excel หรือ Shape file เพื่อนำไปใช้ในโปรแกรมด้าน GIS

โฟลเดอร์และไฟล์ทดสอบโปรแกรม

  • ผมมีไฟล์ตัวอย่างให้ทดสอบ แต่ไฟล์ตัวอย่างจะอยู่ที่โฟลเดอร์ที่ซ่อนของวินโดส์ จากหน้าโปรแกรมหลัก ดูที่ไอคอนล่างสุด จะเห็น “Example folder” ลองดับเบิ้ลคลิกเข้าไป จะเห็นโฟลเดอร์ย่อย “example data” ดับเบิ้ลคลิกเข้าไปจะเห็นไฟล์ดังรูป ตอนแนะนำการใช้โปรแกรมผมจะใช้ไฟล์ในโฟลเดอร์นี่เป็นหลัก
explorer_2016-12-24_13-41-57
Example folder
  • เพื่อให้ดูง่ายๆว่าในไฟล์ใช้ระบบพิกัดเป็นอะไร ชื่อไฟล์จะใส่ระบบพิกัดไปในชื่อไฟล์ด้วย และมีไฟล์ที่ขึ้นต้นด้วย “error….” เป็นตัวอย่างไฟล์ถ้าโปรแกรมเจอไฟล์แบบนี้แล้วจะอ่านไม่ได้ หรืออ่านได้แต่ไม่ประมวลผลเช่น ไฟล์มีค่าพิกัดแค่สองจุด กลายเป็นเส้นตรงไม่สามารถคำนวณหาพท.ได้

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

  • ลอง import ไฟล์แรก “boundary1-geo-wgs84.csv” ค่าพิกัดเป็นค่าพิกัดภูมิศาสตร์ ตอน browse เลือกไฟล์เสร็จจะเห็น preview สำหรับไฟล์นี้จะมีชื่อคอลัมน์ “No,Name,Latitude,Longitude”
Import CSV file
Import CSV file
  • ลองเปิดไฟล์นี้มา โปรแกรมจะแสดงข้อมูลที่อ่านมาลงในตารางดังรูปด้านล่าง โปรแกรมจะจัดคอลัมน์ให้ ถ้าไม่ตรงผู้ใช้ต้องมาจัดเอง หลักๆของการออกแบบคือไม่ว่าไฟล์ CSV ผู้ใช้จะเรียงแบบไหนก็ตามเช่น P,N,E หรือ P,E,N หรือ N,E,P ก็ตามโปรแกรมจะอ่านมาก่อน แล้วผู้ใช้ค่อยมากำหนดตามหลัง เพื่อให้ยืดหยุนและสะดวกต่อการใช้งาน

compute_area_header

  • ตอนนี้ยังไม่ได้กำหนดระบบพิกัด แต่ถ้าเผลอไปคลิกไอคอนรูปเครื่องคิดเลขเพื่อคำนวณโปรแกรมก็คำนวณให้ โดยที่นึกว่าค่าแลตติจูดและลองจิจูดเป็นค่าพิกัดฉาก x,y แต่พื้นที่ก็น้อยจนไม่เห็น ต่อไปจะปรับระบบพิกัด (Coordinate Reference System)
  • ตั้งระบบพิกัดเป็น Geographic โดยคลิกเลือกเป็น “WGS84 / Geographic” แล้วคลิกที่ไอคอนเครื่องคิดเลข จะเห็นผลลัพธ์สองอย่าง คือพื้นที่บนทรงรีคำนวณได้ 85129.837 ตร.ม. ส่วนพื้นที่บนระบบพิกัดฉากได้ 85061.858 ตร.ม. ต่างกัน 67.979 ตร.ม.
Area Computation
Area Computation

ปักหมุดรูปแปลงบน Google Maps

  • ตอนนี้เราทราบแล้วว่าทำไมพื้นที่ถึงไม่เท่ากันเพราะ scale factor นั่นเอง ต่อไปจะลองมาปักหมุดบน google maps ที่ไอคอนด้านขวามือคลิกไปที่ไอคอนของ google maps โปรแกรมจะ plot รูปปิดให้พร้อม ปักหมุดรอบรูปแปลง ตรงกลางจะมีหมุดสีแดง แสดงจุดศูนย์ centroid ของพื้นที่

nvidia-share_2016-12-24_15-49-01

ปักหมุดรูปแปลงบน Google Earth

  • คลิกที่ไอคอน google earth โปรแกรมจะถามฃื่อไฟล์ KML ที่จะจัดเก็บ จากนั้นก็พาบินเข้าพื้นที่ จะเห็นแปลงที่ดินดังรูป จะเห็น label แสดงพื้นที่บนทรงรี (Ellipsoidal Area) ที่จริงตอนผมเปิดเรื่องมาเป็นพื้นที่จริงบนผิวโลก ดังนั้นจะมีกระบวนการแปลงพื้นที่จากผิวทรงรีขึ้นไปบนพื้นผิวภูมิประเทศอีกที เอาไว้ในตอนหน้าครับ

explorer_2016-12-24_15-55-16

  • ในตอนหน้ามาดูการใช้งานโปรแกรม ถ้ากำหนดให้รูปแปลงที่ดินเป็นค่าพิกัดฉากในระบบ UTM บน Indian 1975 หรือ WGS84 จะตั้งค่ายังไง และมาดูเบื้องหลังการคำนวณว่าจากค่าพิกัดภูมิศาสตร์จะคำนวณเนื้อที่ได้อย่างไร

เอกสารอ้างอิง

การเขียนโปรแกรมเพื่อคำนวณระยะทางและอะซิมัท (Distance/Azimuth) บน Ellipsoid ด้วย Lazarus (ตอนที่ 1)

  • ตอนก่อนหน้านี้ ผมเขียนโปรแกรมแปลงพิกัดระหว่าง UTM และ Geographic (Lat/Long) และและถ้าไม่เขียนการหาระยะทางและอะซิมัท (เมื่อกำหนดจุด Latitude, Longitude ให้สองจุด) ก็ดูจะขาดอะไรไปอย่าง

Model ที่ใช้ในการคำนวณ

  • สัณฐานหรือรูปทรงที่ใช้แทนโลก ใช้กันอยู่ 2 แบบ คือ ทรงกลม(Spherical)และทรงรี(Ellipsoid) แต่ทรงรีจะใกล้เคียงกับกับโลกเรามากกว่า ดังนั้นการหาระยะทางและทิศทาง(อะซิมัท) บนทรงรีจะให้ความละเอียดถูกต้องมากกว่า
รูปทรงที่ใช้แทนสัณฐานของโลก
รูปทรงที่ใช้แทนสัณฐานของโลก (ภาพอ้างอิงจาก http://www.codeguru.com/cpp/cpp/algorithms/article.php/c5115/)
  • การคำนวณระยะทางและอะซิมัทบนทรงกลมนั้นง่ายกว่าการคำนวณบนทรงรี การคำนวณระยะทางบนทรงกลมเช่นต้องการทราบระยะทางจากจุด A ไปจุด B เราจะใช้ plane คือแผ่นเรียบตัดผ่านจุด A และจุด B และที่สำคัญคือต้องผ่านจุดศูนย์กลางทรงกลมด้วย เส้นที่เกิดจากแผ่น plane มาตัดกับทรงกลมจะเรียกว่า  Great Circle ระยะที่สั้นที่สุดบนผิวทรงกลมก็คือระยะที่ไปตาม Great circle นี่เอง
Great circle
Great circle อ้างอิงจาก http://www.codeguru.com/cpp/cpp/algorithms/article.php/c5115/
  • ส่วนทรงรีที่สัณฐานใกล้เคียงกับโลกมากกว่า จะยุบลงที่ขั้วโลก โป่งออกที่เส้นศูนย์สูตร เราจะหาระยะทางและอะซิมัทโดยอิงสัณฐานของทรงรีโดยใช้สูตรอันลือชื่อของ Thaddeus Vincenty

สูตรที่ใช้อ้างอิงในการคำนวณ (Vincenty’s formulae)

  • Thaddeus Vincenty เกิดที่โปแลนด์ อพยพเข้าสหรัฐอเมริกาเมื่อสงครามโลกครั้งที่สอง ทำงานให้กับ U.S. Air Force และ  National Geodetic Survey ตามลำดับ ที่เป็นที่รู้จักกันมากที่สุดก็คือ Vincenty’s Formulae นี้เอง เขาเผยแพร่สูตรนี้เมื่อปี 1975 ซึ่งให้ความละเอียดสูงมากจนถึงครึ่งมิลลิเมตร
  • สำหรับสูตรของคุณ Thaddeus Vincenty ถ้าเห็นสูตรยาวๆแล้วสารอะดรีนาลีนหลั่งออกมาก็ ไปดาวน์โหลดได้ที่ http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
  • สำหรับสูตรที่จะนำเสนอในที่นี้ผมแปลมาจาก Javascript ดูได้ที่ http://www.movable-type.co.uk/scripts/latlong-vincenty.html สำหรับสูตรของ Vincenty จะมีการอินทิเกรต แต่ในทางโปรแกรมมิ่งจะใช้เทคนิค iteration แทนทำให้สูตรที่จะโปรแกรมดูกระชับมาก
  • สำหรับสูตรที่นำสรุปจาก Vincenty’s formulae คำนวณหาระยะทางและอะซิมัทจากจุด 2 จุด แสดงให้เห็นดังข้างล่าง
a, b = major & minor semiaxes of the ellipsoid
f = flattening (ab)/a
φ1, φ2 = geodetic latitude
L = difference in longitude
U1 = atan((1−f).tanφ1) (U is ‘reduced latitude’)
U2 = atan((1−f).tanφ2)
λ = L (first approximation)
iterate until change in λ is negligible (e.g. 10-12 ≈ 0.06mm) {
sinσ = √[ (cosU2.sinλ)² + (cosU1.sinU2 − sinU1.cosU2.cosλ)² ] (14)
cosσ = sinU1.sinU2 + cosU1.cosU2.cosλ (15)
σ = atan2(sinσ, cosσ) (16)
sinα = cosU1.cosU2.sinλ / sinσ (17)
cos²α = 1 − sin²α (trig identity; §6)
cos2σm = cosσ − 2.sinU1.sinU2/cos²α (18)
C = f/16.cos²α.[4+f.(4−3.cos²α)] (10)
λ′ = L + (1−C).f.sinα.{σ+C.sinσ.[cos2σm+C.cosσ.(−1+2.cos²2σm)]} (11)
}
u² = cos²α.(a²−b²)/b²
A = 1+u²/16384.{4096+u².[−768+u².(320−175.u²)]} (3)
B = u²/1024.{256+u².[−128+u².(74−47.u²)]} (4)
Δσ = B.sinσ.{cos2σm+B/4.[cosσ.(−1+2.cos²2σm) − B/6.cos2σm.(−3+4.sin²σ).(−3+4.cos²2σm)]} (6)
s = b.A.(σ−Δσ) (19)
α1 = atan2(cosU2.sinλ, cosU1.sinU2 − sinU1.cosU2.cosλ) (20)
α2 = atan2(cosU1.sinλ, −sinU1.cosU2 + cosU1.sinU2.cosλ) (21)

where :

    • s is the distance (in the same units as a & b)
    • α1 is the initial bearing, or forward azimuth
    • α2 is the final bearing (in direction p1→p2)

สิ่งที่ต้องเตรียมก่อนจะโปรแกรม

เริ่มต้นโปรแกรม

  • ผมจะสร้างยูนิตที่ทำหน้าที่คำนวณชื่อ GeodesicCompute (อยู่ในไฟล์ geodesiccompute.pas) โค๊ดดังที่กล่าวไปแล้วจาก Javascript แต่ดัดแปลงให้เป็นโค๊ด OOP ในยูนิตนี้มี 2 class คือ TGeodesicInverse ทำหน้าที่คำนวณระยะทางและอะซิมัทเมื่อป้อนค่าพิกัดเป็น Latitude,Longitude 2 จุด เรียกการคำนวณนี้ว่า Inverse
  • และอีกยูนิตหนึ่งคือ TGeodesicDirect(อยู่ในไฟล์ geodesiccompute.pas) ทำหน้าที่คำนวณหาค่าพิกัดจุดที่สอง (Lat/Long) เมื่อกำหนดค่าพิกัดจุดที่หนึ่งให้และ กำหนดอะซิมัทและระยะทาง เรียกการคำนวณนี้สั้นๆว่า Direct

GeodesicCompute Unit

  • โค๊ดที่ผมเขียนดัดแปลงมาจาก Javascript เป็นไปดัง source code ด้านล่าง
unit GeodesicCompute;

{$mode objfpc}{$H+}

interface

uses
Classes, SysUtils, GeoEllipsoids, Math;

type
TZoneHemisphere = (hemEast, hemWest, hemNorth, hemSouth);

TCoordinate = record
X : double;
Y : double;
end;

//
//  Vincenty Inverse Solution of Geodesics on the Ellipsoid.
// Calculate geodesic distance (in m) between two points
// specified by latitude/longitude (in numeric degrees) using
// Vincenty inverse formula for ellipsoids.
//

TGeodesicInverse = class
private
FEllipsoid : TEllipsoid;
FGeoCoor1 : TCoordinate;
FGeoCoor2 : TCoordinate;
FEllDist, FInitAzi, FFinalAzi : double;
procedure SetGeoCoordinate1 (const geocoor1 : TCoordinate);
procedure SetGeoCoordinate2 (const geocoor2 : TCoordinate);//input
procedure SetEllipsoid(const Ellipsoid : TEllipsoid);//input
function GeodesicInverse(lat1, lon1, lat2, lon2 : double) : double;
public
//input
property GeoCoordinate1 : TCoordinate write SetGeoCoordinate1;
property GeoCoordinate2 : TCoordinate write SetGeoCoordinate2;
//output
property Distance : double read FEllDist;
property InitialAzimuth : double read FInitAzi;
property FinalAzimuth : double read FFinalAzi;
//input
property Ellipsoid : TEllipsoid write SetEllipsoid;
//main computation
procedure Compute;

constructor Create;
destructor Destroy; override;
end;

// Calculate destination point given start point lat/long
// (numeric degrees), azimuth or bearing (numeric degrees)
// & distance (in m).
//
// from: Vincenty direct formula - T Vincenty,
// "Direct and Inverse Solutions of Geodesics on the
// Ellipsoid with application of nested equations",
// Survey Review, vol XXII no 176, 1975
// http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf

TGeodesicDirect = class
private
FEllipsoid : TEllipsoid;
FGeoCoor1 : TCoordinate;
FGeoCoor2 : TCoordinate;
FAzimuth : double;
FDistance : double;
procedure SetGeoCoordinate1 (const geocoor1 : TCoordinate);
procedure SetEllipsoid(const Ellipsoid : TEllipsoid);
procedure SetAzimuth(const azi : double);
procedure SetDistance(const dist : double);
function GeodesicDirect(lat, lon, azimuth, distance : double) : TCoordinate;
public
//input
property Ellipsoid : TEllipsoid read FEllipsoid write SetEllipsoid;
property GeoCoordinate1 : TCoordinate write SetGeoCoordinate1;
property Azimuth : double write SetAzimuth;
property Distance : double write SetDistance;
//Output
property GeoCoordinate2 : TCoordinate read FGeoCoor2;
//main computation
procedure Compute;

constructor Create;
destructor Destroy; override;
end;
implementation

// TGeodesicInverse
procedure TGeodesicInverse.SetGeoCoordinate1 (const geocoor1 : TCoordinate);
begin
FGeoCoor1 := geocoor1;
end;

procedure TGeodesicInverse.SetGeoCoordinate2 (const geocoor2 : TCoordinate);
begin
FGeoCoor2 := geocoor2;
end;

procedure TGeodesicInverse.SetEllipsoid(const Ellipsoid : TEllipsoid);
begin
FEllipsoid.EllipsoidName := Ellipsoid.EllipsoidName;
FEllipsoid.MajorAxis := Ellipsoid.MajorAxis;
FEllipsoid.InvFlattening := Ellipsoid.InvFlattening;
end;

procedure TGeodesicInverse.Compute;
begin
FEllDist := GeodesicInverse(FGeoCoor1.Y, FGeoCoor1.X, FGeoCoor2.Y, FGeoCoor2.X);
end;

function TGeodesicInverse.GeodesicInverse(lat1, lon1, lat2, lon2 : double) : double;
var
a, b, f , e : double;
dlat1, dlon1, dlat2, dlon2 : double;
L, U1, U2 : double;
sinU1, cosU1, sinU2, cosU2 : double;
lambda, lambdaP, sinlambda, coslambda, sinsigma : double;
cossigma, sinalpha, sigma, cosSqalpha, cos2sigmaM : double;
cossigmaM, C, uSq, AA, BB, deltaSigma : double;
az1, az2 : double;
iterLimit : integer;
diff : double;
begin
dlat1 := PI/180 * lat1;
dlon1 := PI/180 * lon1;
dlat2 := PI/180 * lat2;
dlon2 := PI/180 * lon2;

a := FEllipsoid.MajorAxis;
f := 1/FEllipsoid.InvFlattening;
b := a*(1 - f);
L := dlon2-dlon1;
U1 := arctan((1-f) * tan(dlat1));
U2 := arctan((1-f) * tan(dlat2));
sinU1 := sin(U1);
cosU1 := cos(U1);
sinU2 := sin(U2);
cosU2 := cos(U2);
lambda := L;
iterLimit := 100;

repeat
sinlambda := sin(lambda);
coslambda := cos(lambda);
sinsigma := sqrt((cosU2*sinLambda) * (cosU2*sinlambda) +
(cosU1*sinU2-sinU1*cosU2*coslambda) *
(cosU1*sinU2-sinU1*cosU2*coslambda));
if (sinsigma = 0) then
begin
result := 0; // co-incident points
exit;
end;
cossigma := sinU1*sinU2 + cosU1*cosU2*coslambda;
sigma := arctan2(sinsigma, cossigma);
sinalpha := cosU1 * cosU2 * sinlambda / sinsigma;
cosSqAlpha := 1 - sinAlpha*sinAlpha;
cos2SigmaM := cosSigma - 2*sinU1*sinU2/cosSqAlpha;
if (IsNaN(cos2sigmaM)) then cosSigmaM := 0; // equatorial line: cosSqAlpha=0 (6)
C := f/16*cosSqalpha*(4+f*(4-3*cosSqalpha));
lambdaP := lambda;
lambda := L + (1-C) * f * sinAlpha *
(sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*
(-1+2*cos2SigmaM*cos2SigmaM)));
dec(iterLimit);
diff := abs(lambda-lambdaP);
until ((diff < 1e-12) and (iterLimit > 0));

if (iterLimit = 0) then
begin
result := NaN; // formula failed to converge
exit;
end;

uSq := cosSqAlpha * (a*a - b*b) / (b*b);
AA := 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
BB := uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));
deltaSigma := BB*sinsigma*(cos2sigmaM+BB/4*
(cossigma*(-1+2*cos2sigmaM*cos2sigmaM)-
BB/6*cos2sigmaM*(-3+4*sinsigma*sinsigma)*
(-3+4*cos2sigmaM*cos2sigmaM)));
result := b*AA*(sigma-deltaSigma); //return Ellipse Distance.
//initial azimuth.
az1 := arctan((cosU2*sinlambda)/(cosU1*sinU2 - sinU1*cosU2*coslambda));
//final azimuth.
az2 := arctan((cosU1*sinlambda)/(-sinU1*cosU2 + cosU1*sinU2*coslambda));

FInitAzi := az1 * 180/PI;
if (FInitAzi < 0) then FInitAzi := 360 + FInitAzi;

FFinalAzi := az2 * 180/PI;
if (FFinalAzi < 0) then FFinalAzi := 360 + FFinalAzi;
end;

constructor TGeodesicInverse.Create;
begin
inherited create;
FEllipsoid := TEllipsoid.Create;
end;

destructor TGeodesicInverse.Destroy;
begin
FEllipsoid.Free;
inherited Destroy;
end;

// TGeodesicDirect
procedure TGeodesicDirect.SetGeoCoordinate1 (const geocoor1 : TCoordinate);
begin
FGeoCoor1 := geocoor1;
end;

procedure TGeodesicDirect.SetAzimuth (const azi : double);
begin
FAzimuth := azi;
end;

procedure TGeodesicDirect.SetDistance (const dist : double);
begin
FDistance := dist;
end;

procedure TGeodesicDirect.SetEllipsoid(const Ellipsoid : TEllipsoid);
begin
FEllipsoid.EllipsoidName := Ellipsoid.EllipsoidName;
FEllipsoid.MajorAxis := Ellipsoid.MajorAxis;
FEllipsoid.InvFlattening := Ellipsoid.InvFlattening;
end;

procedure TGeodesicDirect.Compute;
begin
FGeoCoor2 := GeodesicDirect(FGeoCoor1.Y, FGeoCoor1.X, FAzimuth, FDistance);
end;

function TGeodesicDirect.GeodesicDirect(lat, lon, azimuth, distance : double) : TCoordinate;
var
a, b, f , e, s : double;
dlat1, dlon1, dlat2, dlon2 : double;
sinU1, alpha1, sinAlpha, sinAlpha1, cosAlpha1, tanU1, cosU1, sigma1 : double;
cos2SigmaM, sinSigma, cosSigma, cosSqAlpha, uSq, deltaSigma : double;
AA, BB, sigma, sigmaP, tmp, lambda, C, L, lat2, revAz : double;

begin
dlat1 := PI/180 * lat;
dlon1 := PI/180 * lon;

a := FEllipsoid.MajorAxis;
f := 1/FEllipsoid.InvFlattening;
b := a*(1 - f);

s := distance;
alpha1 := azimuth * PI/180.0;
sinAlpha1 := sin(alpha1);
cosAlpha1 := cos(alpha1);

tanU1 := (1-f) * tan(dlat1);
cosU1 := 1 / sqrt((1 + tanU1*tanU1));
sinU1 := tanU1*cosU1;
sigma1 := arctan2(tanU1, cosAlpha1);
sinAlpha := cosU1 * sinAlpha1;
cosSqAlpha := 1 - sinAlpha*sinAlpha;
uSq := cosSqAlpha * (a*a - b*b) / (b*b);
AA := 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
BB := uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));

sigma := s / (b*AA);
sigmaP := 2*PI;
while (abs(sigma-sigmaP) > 1e-12) do
begin
cos2SigmaM := cos(2*sigma1 + sigma);
sinSigma := sin(sigma);
cosSigma := cos(sigma);
deltaSigma := BB*sinSigma*(cos2SigmaM+BB/4*
(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-
BB/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*
(-3+4*cos2SigmaM*cos2SigmaM)));
sigmaP := sigma;
sigma := s / (b*AA) + deltaSigma;
end;

tmp := sinU1*sinSigma - cosU1*cosSigma*cosAlpha1;
dlat2 := arctan2(sinU1*cosSigma + cosU1*sinSigma*cosAlpha1,
(1-f)*sqrt(sinAlpha*sinAlpha + tmp*tmp));
lambda := arctan2(sinSigma*sinAlpha1,
cosU1*cosSigma - sinU1*sinSigma*cosAlpha1);
C := f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha));
L := lambda - (1-C) * f * sinAlpha *
(sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*
(-1+2*cos2SigmaM*cos2SigmaM)));

revAz := arctan2(sinAlpha, -tmp); // final bearing
result.X := lon + L * 180/PI;
result.Y := dlat2 * 180/PI;
end;

constructor TGeodesicDirect.Create;
begin
inherited create;
FEllipsoid := TEllipsoid.Create;
end;

destructor TGeodesicDirect.Destroy;
begin
FEllipsoid.Free;
inherited Destroy;
end;

end.

  • สำหรับ source code ที่แสดงดังข้างต้น เป็น class ง่ายๆ เช่น TGeodesicInverse จะมี property สำหรับ input อยู่ 3 ตัวคือ ค่ารูปทรงสัณฐาน Ellipsoid, GeoCoordinate1 และ GeoCoordinate2 เมื่อ input ค่า lat,long 2 ค่าพิกัดให้เสร็จ แล้วก็เรียกเมธอด Compute เพื่อทำการคำนวณ จากนั้นเรียกค่าจาก 2 property คือ Azimuth และ Distance ซึ่งจะเป็นผลลัพธ์การคำนวณ
  • ส่วนคลาส TGeodesicDirect ก็เหมือนๆกัน มี property 4 ตัว คือ Ellipsoid, GeoCoordinate1, Azimuth และ Distance เมื่อรับค่ามาก็จะเรียกเมธอด Compute และให้ผลลัพธ์เป็นค่าพิกัด (lat,long) ของจุดที่ต้องการหา

ปัญหาความเข้ากันของ Widget ในแต่ละ Platform

  • รูปที่จะแสดงให้ดูเป็นปัญหาของ Widget ผมเพิ่งจะประสบความสำเร็จติดตั้ง Widget ของ Qt บน PCLinuxOS เมื่อรันโปรแกรมด้วยชุด source code เดียวกันเปรียบเทียบกับ Win32 บนวินโดส์ และ Gtk2 บน Ubuntu  พบว่า ขนาดของ Lable, EditBox, ComboBox, Groupbox ทั้งหลายบน Gtk2 มีขนาดใหญ่กว่าเพื่อน แต่พบว่าขนาดของ object พวกนี้บน Win32 และ Qt มีขนาดใกล้เคียงกัน เพื่อให้เข้ากับ Gtk2 จึงต้องมาขยับ เช่นตัว Label กับ EditBox จะต้องอยู่ห่างกันพอสมควร ไม่งั้นตัว Label จะล้ำเข้าไปหา EditBox แต่ความสวยงามดูๆ QT จะสวยเนียนกว่าเพื่อน รองลงมาเป็น Gtk2 กับ Win32 ดูรูปผลลัพธ์ของโปรแกรมที่ผมเขียนในตอนนี้
Geodesic computation on Qt widget.
โปรแกรมคำนวณหาระยะทางและอะซิมัทบน Ellipsoid รันด้วย Qt (PCLinuxOS)
โปรแกรมคำนวณหาระยะทางและอะซิมัทบน Ellipsoid รันด้วย Gtk2 (Ubuntu)
โปรแกรมคำนวณหาระยะทางและอะซิมัทบน Ellipsoid รันด้วย Gtk2 (Ubuntu)
โปรแกรมคำนวณหาระยะทางและอะซิมัทบน Ellipsoid รันด้วย Windows
โปรแกรมคำนวณหาระยะทางและอะซิมัทบน Ellipsoid รันด้วย Windows
  • ทิ้งท้ายกันก่อนในตอนนี้ ติดตามกันตอนหน้า