คำนวณพื้นที่ polygon จากไฟล์ในระบบพิกัดฉาก UTM
- ตอนแรกได้สาธยายเรื่อง Scale factor ของเส้นโครงแผนที่ที่มีผลต่อการรังวัดของช่างสำรวจไปแบบพอหอมปากหอมคอ และท้ายสุดได้ใช้โปรแกรมเปิดไฟล์ CSV ที่เก็บค่าพิกัดแบบภูมิศาสตร์ของรูปปิด polygon คำนวณหาพื้นที่บนระนาบ UTM และพื้นที่บนทรงรี ก็เห็นไปแล้วว่าสองพื้นที่นี้มีความแตกต่างกัน เพราะสืบเนื่องจาก Scale factor เบื้องหลังการทำงานของโปรแกรมก็คือ เมื่ออ่านพิกัดภูมิศาสตร์เข้ามาก็ใช้ไลบรารี Proj4 ทำการแปลงย้อนกลับไปค่าพิกัดฉาก UTM แล้วทำการคำนวณหาพิกัดทั้งสองแบบ
- ตัวอย่างต่อไปจะอ่านไฟล์ที่เก็บในรูประบบพิกัดฉาก UTM บน WGS84 ที่โปรแกรมที่ทูลส์บาร์ด้านขวา คลิก “Import” เพื่อเปิดนำเข้าไฟล์ CSV
- ไปที่โฟลเดอร์โปรแกรมที่ผมกล่าวไปแล้วในตอนที่ 1 เลือกไฟล์ “boundary1-utm47n-wgs84.csv“
- โปรแกรมจะ preview ให้ดูก่อนว่าที่เลือกมานั้นใช่หรือไม่ ถ้าใช่ก็คลิก “OK” ไปต่อ
- โปรแกรมจะนำค่าพิกัดของรูป polygon มาใส่ไว้ในตารางข้อมูล พร้อมกับจัดคอลัมน์ Northing/Easting ให้ อย่าลืมตรวจสอบด้วยว่าถูกต้องไหม
- ไฟล์นี้ทราบดีว่าระบบพิกัดฉากเป็น UTM Zone 47N WGS84 ตั้งระบบพิกัดที่ Reference Coordinate System ดังรูป เลือก “WGS84 / UTM Zone 47N” จากนั้นคลิกที่ไอคอนรูปเครื่องคิดเลข
- จะได้ผลลัพธ์ จะเห็นว่าพื้นที่บนทรงรี = 85129.846 ตร.ม. พื้นที่บน UTM = 85061.867 ตร.ม.ต่างกัน 67.979 ตร.ม. ขอพักผลลัพธ์ตรงนี้ไว้ก่อน
เบื้องหลังงานคำนวณ
- การคำนวณพื้นที่ของระบบพิกัดฉาก จะใช้สูตรเชือกรองเท้า (Shoelace’s formula) ที่ผมกล่าวไปในตอนที่ 1 ง่ายๆครับ ยิ่งถ้าใช้ excel แล้วยิ่งง่าย จากข้อมูลข้างต้นผมแปลงค่าพิกัดรูปปิดทั้งหมดไปยังค่าพิกัดภูมิศาสตร์ด้วยไลบรารี Proj4 พอได้ค่าพิกัดภูมิศาสตร์ (แลตติจูด ลองจิจูด)มาแล้ว จะคำนวณหาพืนที่ได้อย่างไร ใช้สูตรเชื่อกผูกรองเท้าไม่ได้แล้ว และพื้นผิวไม่ใช่ระนาบราบแบบ UTM แต่เป็นพื้นผิวโค้ง
- สำหรับสูตรที่ใช้คำนวณพื้นที่บนทรงรี ค่อนข้างซับซ้อนเป็นเรื่อง math หนักๆ ถ้าสนใจลองอ่านของ Charles F. F. Karney, Geodesics on an ellipsoid of revolution,Feb 2011. http://arxiv.org/abs/1102.1215 ตำราสามารถดาวน์โหลดได้ที่หน้าเวปเพจนี้ แต่สำหรับผมไม่ไหวครับหนักเกิน Charles F.F. Karney เป็นผู้พัฒนา GeographicLib มีโมดูล PolygonArea ผมจะใช้ฟังก์ชันนี้ลองคำนวณหาพื้นที่เพื่อเปรียบเทียบกับวิธีที่ผมจะใช้
- วิธีการที่ผมจะใช้คำนวณหาพื้นที่จะใช้เส้นโครงแผนที่ (Map Projection) มาช่วย แต่เป็นเส้นโครงแผนที่ที่รักษาพื้นที่ (ไม่รักษารูปร่าง) เลือกใช้ “Lambert Azimuthal Equal Area” ดูชื่อแล้วเหมือนจะเป็น 2 in 1 คือรักษาทิศทาง (Azimuthal) จากจุดศูนย์กลางที่สัมผัสออกไปยังจุดต่างและยังรักษาพื้นที่ (Equal Area) ด้วย คือได้พื้นที่จริงๆมาแน่นอนบนทรงรี แต่ความรู้สึกผมถ้าพื้นที่ประมาณไม่เกิน 100 ตร.กม. (ประมาณ 10 กม. x 10 กม.) น่าจะยังรักษารูปร่างได้ดี
- จากรูปด้านบนรูปซ้ายสุด จะเป็นการเอา plane ราบๆ ไปสัมผัสที่ขั้วโลกเหนือ การสัมผัสเลือกสัมผัสได้ตามที่ต้องการ ดูสามรูปด้านขวา แต่ที่ผมจะใช้งานจะเลือกเอา plane มาสัมผัสที่จุด centroid คือจุดศูนย์ถ่วงของพื้นที่ จากนั้นตั้ง latitude of origin, central of meridian ที่ค่าพิกัดภูมิศาสตร์ของจุดสัมผัสนี้ และตั้ง False Northing, False Easting เป็น 0.0
- จากนั้น จะใช้ไลบรารี Proj4 แปลงพิกัดจากค่าพิกัดภูมิศาสตร์บนทรงรีเป็นค่าพิกัดบนระบบพิกัดฉากของ Lambert Azimuthal Equal Area (LAEA) เมื่อได้ค่าพิกัดฉากแล้วก็สามารถนำไปคำนวณด้วยสูตรเชือกผูกรองเท้าได้
ข้ามการคำนวณ Grid Scale Factor
- ปกติถ้าเป็นระยะทางการจะทอนจากระบบพิกัดฉากไปยังระยะทางบนทรงรี จะต้องนำระยะราบมาหารด้วย GSF แต่ถ้าเป็นพื้นที่บนระบบพิกัดฉากต้องการทอนไปเป็นพืนที่บนทรงรีจะต้องหารด้วยค่า GSF ยกกำลังสอง
- พท.บนทรงรี = พื้นที่บนระบบพิกัดฉาก UTM / GSF²
- แต่เมื่อเลือกใช้เส้นโครงแผนที่ LAEA แล้ว ไม่ต้องใช้ GSF แต่เราจะคำนวณหา GSF ในภายหลังเพื่อเปรียบเทียบว่าผลลัพธ์พื้นที่ที่คำนวณด้วยสองวิธีนี้แตกต่างกันอย่างไรบ้าง
เพิ่มเติม
- เส้นโครงแผนที่ที่เราใช้กันอยู่บ่อยๆคือ UTM จะเป็นเส้นโครงแผนที่รักษารูปร่าง (Conformal) รูปร่างจะเหมือนเดิมไม่ว่าขยับไปอยู่ตรงไหน แต่ข้อเสียคือขนาดจะเป็นเปลี่ยนไปคืออาจจะใหญ่ขึ้นหรือเล็กลง เอาง่ายๆคือถ้ารูปร่างเป็นวงกลมไม่ว่าจะขยับไปบนลงล่าง รูปร่างจะยังเป็นวงกลม แต่อาจมีขนาดที่ใหญ่ขึ้น ด้วยเหตุนี้แผนที่บนเส้นโครงแผนที่นี้จึงดูเกาะ Greenland มีขนาดใหญ่โตมโหฬาร ทั้งๆที่ขนาดเล็กกว่าทวีปอาฟริกาหลายเท่า ดูรูปด้านล่างจะเข้าใจดี
- ส่วนเส้นโครงแผนที่ที่รักษาพื้นที่ (Equal Area) ไม่ว่าจะขยับไปอยู่ตรงไหน ขนาดพื้นที่ของมันยังเท่าเดิม แต่รูปร่างจะเพี้ยนไป วงกลมจะกลายเป็นวงรี ดูรูปด้านล่างจะเห็นว่าเกาะ Greenland มีขนาดนิดเดียว
การส่งออกไฟล์ผลลัพธ์ (Export file)
- เมื่อคำนวณแล้วต้องการผลลัพธ์ที่เป็นรูปธรรม นอกจากปักหมุดที่ google maps และ google earth แล้ว โปรแกรมเตรียมส่งออกไฟล์ผลลัพธ์ไปได้สามรูปแบบคือ CSV, Excel และ Shape file
- มาดูการส่งออกไฟล์ Excel คลิกที่ไอคอน “Export”
- ไดอะล๊อกจะถามชื่อไฟล์ แต่ก่อนอื่นเลือกชนิดไฟล์ที่จะส่งออกก่อน ในที่นี้เลือก Excel
- เมื่อจัดเก็บไฟล์ excel แล้ว ก็ถึงเวลาเปิดมาดูไฟล์ผลลัพธ์ ผมใช้ Libreoffice Calc เปิดไฟล์ excel ตรงด้านล่างที่แสดง sheet จะมีอยู่สาม sheet คือ “area”, “projection” และ “ellipsoid” อันดับแรกดูที่ sheet “area” ก่อน
- ลองเลื่อนลงมาที่ด้านล่าง
- สองคอลัมน์ด้านขวา “Cross-Multiplying” จะเป็นคอลัมน์สำหรับคูนค่า X,Y ไขว์ คลองคลิกดูที่เซลล์จะเห็นสูตรที่โปรแกรมเขียนฝังไว้ในไฟล์ จากนั้นจะทำการบวกผลรวมแต่ละคอลัมน์ จับมาลบกันแล้วคูนด้วย 0.5 จะได้พื้นที่เป็นตร.ม. ออกมา
- พื้นที่ค่าที่คำนวณได้ตัวนี้เป็นพื้นที่บนเส้นโครงแผนที่ Lambert Azimuthal Equal Area ซึ่งเป็นระบบพิกัดฉากแบบหนึ่ง จะสังเกตว่าค่าที่คำนวณในโปรแกรม excel หรือที่ผมใช้ Libreoffice Calc จะเป็นค่าที่คำนวณจากในโปรแกรมเสปรตชีตในนี้นี้เอง ไม่ได้เอาค่าที่คำนวณจากโปรแกรม “Area” มาใส่แต่อย่างใด
- มีอีกสอง sheet คือ “projection” เป็นรายละเอียดของเส้นโครงแผนที่ที่เราใช้ และอีก sheet “ellipsoid” เป็นค่าพารามิเตอร์ทรงรี WGS84 ผมเขียนไว้เพื่อต้องการนำค่า flattening (f) ไปใช้ใน sheet “area” เพื่อคำนวณหารัศมีของทรงรี ณ ค่าแลตติจูดที่กำหนด มาดู sheet “ellipsoid” กัน
- ก็เป็นค่าพารามิเตอร์มาตรฐานของ WGS84 ได้แก่ค่า a, b, f, e, e’ ดู sheet ถัดไปคือ “projection”
- สำหรับเส้นโครงแผนที่ Lambert Azimuthal Equal Area (LAEA) ผมใช้ Plane มาสัมผัสที่จุดศูนย์ถ่วงของพื้นที่ จากนั้นจะโปรเจคจุดต่างๆบนพื้นทรงรีมาที่ plane นี้
- จุดศูนย์ถ่วงนี้คำนวนในโปรแกรม “Area” ดังนั้นจะเห็นค่า Lattitude of origin = 13.943866348 ค่า central meridian = 99.067227167 และจะเห็น ESRI WKT ซึ่งข้อความนี้จะถูกเขียนไปตอนส่งออกไฟล์ Shape file นามสกุล .prj ซึ่งไฟล์ .prj คน GIS คงทราบดีว่าเป็นไฟล์กำหนดระบบพิกัด สามารถอ่านได้โดยโปรแกรม GIS ทั่วๆไป และผมแถม Proj4 string ด้วยค่านี้ผมใช้งานในโปรแกรมตอนแปลงพิกัด สังเกตดูว่ากระชับกว่า ESRI WKT
คำนวณหาพื้นที่บนพื้นผิวภูมิประเทศ (Surface Area)
- มาดูขั้นตอนการหาพื้นที่ Surface Area จากพื้นที่บนทรงรี (Ellipsoidal area) ขั้นตอนนี้คือการหาค่า Elevation Scale Factor(ESF) นั่นเอง แต่ที่นี้เราย้ายการคำนวณมายัง Excel
- วิธีการคำนวณ ESF แบบละเอียดผมได้กล่าวไปในตอนที่ 1 แล้ว เนื่องจาก ESF ต้องการค่าระดับ เราต้องรู้ค่าระดับเฉลี่ย ในพื้นที่ตัวอย่างนี้ผมทราบว่าค่าระดับเทียบกับน้ำทะเลปานกลางประมาณ 180 เมตร ผมจัดการป้อนไปในช่องสีเหลืองดังรูป จากนั้นผมเปิดโปรแกรม “EGM” เพื่อคำนวณความสูงจีออยด์
- นำค่า -34.7246 มาป้อนที่ช่องสีเหลือง และสูตรที่ฝังในสเปรตชึตจะคำนวณค่า Ellipsoid Height = 180 + (-32.7246) = 145.275 เมตร สุดท้ายสูตรจะคำนวณหาค่ารัศมีทรงรี ณ แลตติจูด = 13.943866348 ได้เท่าก้บ 6376911 เมตร
- คำนวณหาค่า ESF = R/(R+h) = 6376911/(6376911+145.275) = 0.999977219
- คำนวณหาพื้นทีจริง Ground-base area หรือ Surface area = Ellipsoidal Area / ESF² (พืนที่บนทรงรีหารด้วย ESF ยกกำลังสอง)
- พื้นที่จริง = 85129.846 / (0.999977219²) = 85133.725 ตร.ม.
ทบทวนสูตรการแปลงพื้นที่
- แปลงจากพื้นที่บนระบบพิกัดฉาก (Grid-based Area) ไปยังพื้นที่บนทรงรี และจากพื้นที่ทรงรี (Ellipsoidal Area) ไปยังพื้นที่บนพื้นที่ผิวภูมิประเทศบนพื้นโลก (Surface Area)
- หรือถ้าทราบ ESF แบะ GSF ก็รวบรัดดังนี้
- ช่างรังวัดกรมที่ดินกับช่างรังวัดเอกชน เห็นสูตรนี้แล้ว คงคุ้นๆกัน ESF ก็คือ ค่าคูนสัมประสิทธ์ (C) นั่นเอง ส่วน GSF คือค่าคูนมาตราส่วน(K)
ตรวจสอบความถูกต้องของพื้นที่บนระบบพิกัดฉาก Lambert Azimuthal Equal Area
- ตอนที่สั่งคำนวณจะเห็นผลลัพธ์บนหน้าต่างของโปรแกรม ขอทบทวนอีกครั้ง พื้นที่บนทรงรี = 85129.846 ตร.ม. และ พื้นที่บน UTM = 85061.867 ตร.ม. ต่างกัน 67.979 ตร.ม.
- ค่าพิกัดจุดศูนย์ถ่วง แลตติจูด = 13.943866348 ค่าลองจิจูด = 99.067227167 ผมหา GSF ได้ 0.99960065 ลองสูตรแรก Ellipsoidal Area = Grid-based Area / GSF² = 85061.867/0.99960065² = 85129.847 ตร.ม. ต่างจากค่าที่คำนวณด้วยเส้นโครงแผนที่ LAEA 85129.846 ตร.ม. น้อยมากตำแหน่งทศนิยมที่ 3 สรุปว่าวิธีการใช้เส้นโครงแผนที่ Lambert Azimuthal Equal Area ใช้ได้ดี
- ขอจบตอนนี้ พบกับตอนที่ 3 ตอนสุดท้าย จะมาเก็บตกกัน ถ้าพื้นที่อยู่ในเขตแบ่งโซน UTM จะทำยังไง และลองค่าพิกัดยูทีเอ็มของรูปปิดในระบบของ Indian 1975 ดู