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

คำนวณพื้นที่ polygon จากไฟล์ในระบบพิกัดฉาก UTM

  • ตอนแรกได้สาธยายเรื่อง Scale factor ของเส้นโครงแผนที่ที่มีผลต่อการรังวัดของช่างสำรวจไปแบบพอหอมปากหอมคอ และท้ายสุดได้ใช้โปรแกรมเปิดไฟล์ CSV ที่เก็บค่าพิกัดแบบภูมิศาสตร์ของรูปปิด polygon คำนวณหาพื้นที่บนระนาบ UTM และพื้นที่บนทรงรี ก็เห็นไปแล้วว่าสองพื้นที่นี้มีความแตกต่างกัน เพราะสืบเนื่องจาก  Scale factor เบื้องหลังการทำงานของโปรแกรมก็คือ เมื่ออ่านพิกัดภูมิศาสตร์เข้ามาก็ใช้ไลบรารี  Proj4 ทำการแปลงย้อนกลับไปค่าพิกัดฉาก UTM แล้วทำการคำนวณหาพิกัดทั้งสองแบบ
  • ตัวอย่างต่อไปจะอ่านไฟล์ที่เก็บในรูประบบพิกัดฉาก UTM บน WGS84 ที่โปรแกรมที่ทูลส์บาร์ด้านขวา คลิก “Import” เพื่อเปิดนำเข้าไฟล์ CSV

compute_area_import2

 

  • ไปที่โฟลเดอร์โปรแกรมที่ผมกล่าวไปแล้วในตอนที่ 1 เลือกไฟล์ “boundary1-utm47n-wgs84.csv

 

surveyor-pocket-tools_2016-12-25_10-44-50

  • โปรแกรมจะ  preview ให้ดูก่อนว่าที่เลือกมานั้นใช่หรือไม่  ถ้าใช่ก็คลิก “OK” ไปต่อ

surveyor-pocket-tools_2016-12-25_10-56-05

  • โปรแกรมจะนำค่าพิกัดของรูป polygon มาใส่ไว้ในตารางข้อมูล พร้อมกับจัดคอลัมน์ Northing/Easting ให้ อย่าลืมตรวจสอบด้วยว่าถูกต้องไหม

compute_area_example2

  • ไฟล์นี้ทราบดีว่าระบบพิกัดฉากเป็น UTM Zone 47N WGS84 ตั้งระบบพิกัดที่ Reference Coordinate System ดังรูป เลือก “WGS84 / UTM Zone 47N” จากนั้นคลิกที่ไอคอนรูปเครื่องคิดเลข

compute_area_calc2

  • จะได้ผลลัพธ์ จะเห็นว่าพื้นที่บนทรงรี = 85129.846 ตร.ม. พื้นที่บน UTM = 85061.867 ตร.ม.ต่างกัน 67.979 ตร.ม. ขอพักผลลัพธ์ตรงนี้ไว้ก่อน

compute_area_result2

เบื้องหลังงานคำนวณ

  • การคำนวณพื้นที่ของระบบพิกัดฉาก จะใช้สูตรเชือกรองเท้า (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 กม.) น่าจะยังรักษารูปร่างได้ดี
image4
ภาพจากเวป http://ayresriverblog.com/2011/05/19/the-world-is-flat/
  • จากรูปด้านบนรูปซ้ายสุด จะเป็นการเอา 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 มีขนาดใหญ่โตมโหฬาร ทั้งๆที่ขนาดเล็กกว่าทวีปอาฟริกาหลายเท่า ดูรูปด้านล่างจะเข้าใจดี
ภาพจาก http://geokov.com/education/map-projection.aspx
ภาพจาก http://geokov.com/education/map-projection.aspx
  • ส่วนเส้นโครงแผนที่ที่รักษาพื้นที่ (Equal Area) ไม่ว่าจะขยับไปอยู่ตรงไหน ขนาดพื้นที่ของมันยังเท่าเดิม แต่รูปร่างจะเพี้ยนไป วงกลมจะกลายเป็นวงรี ดูรูปด้านล่างจะเห็นว่าเกาะ Greenland มีขนาดนิดเดียว
450px-tissot_indicatrix_world_map_lambert_cyl_equal-area_proj-svg
ภาพจาก https://en.wikipedia.org/wiki/Lambert_cylindrical_equal-area_projection

การส่งออกไฟล์ผลลัพธ์ (Export file)

  • เมื่อคำนวณแล้วต้องการผลลัพธ์ที่เป็นรูปธรรม นอกจากปักหมุดที่ google maps และ google earth แล้ว โปรแกรมเตรียมส่งออกไฟล์ผลลัพธ์ไปได้สามรูปแบบคือ CSV, Excel และ Shape file
  • มาดูการส่งออกไฟล์ Excel คลิกที่ไอคอน “Export”

compute_area_export

  • ไดอะล๊อกจะถามชื่อไฟล์ แต่ก่อนอื่นเลือกชนิดไฟล์ที่จะส่งออกก่อน ในที่นี้เลือก Excel

surveyor-pocket-tools_2016-12-25_15-21-16

  • เมื่อจัดเก็บไฟล์ excel แล้ว ก็ถึงเวลาเปิดมาดูไฟล์ผลลัพธ์ ผมใช้ Libreoffice Calc เปิดไฟล์ excel ตรงด้านล่างที่แสดง sheet จะมีอยู่สาม sheet คือ “area”, “projection” และ “ellipsoid” อันดับแรกดูที่ sheet “area” ก่อน

soffice-bin_2016-12-25_15-40-56

  • ลองเลื่อนลงมาที่ด้านล่าง

soffice-bin_2016-12-25_15-48-19

  • สองคอลัมน์ด้านขวา “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” กัน

soffice-bin_2016-12-25_19-47-01

  • ก็เป็นค่าพารามิเตอร์มาตรฐานของ WGS84 ได้แก่ค่า a, b, f, e, e’ ดู  sheet ถัดไปคือ “projection”

soffice-bin_2016-12-25_19-49-47

  • สำหรับเส้นโครงแผนที่ 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

compute_area_excel1

  • วิธีการคำนวณ ESF แบบละเอียดผมได้กล่าวไปในตอนที่ 1 แล้ว  เนื่องจาก ESF ต้องการค่าระดับ เราต้องรู้ค่าระดับเฉลี่ย ในพื้นที่ตัวอย่างนี้ผมทราบว่าค่าระดับเทียบกับน้ำทะเลปานกลางประมาณ 180 เมตร ผมจัดการป้อนไปในช่องสีเหลืองดังรูป จากนั้นผมเปิดโปรแกรม “EGM” เพื่อคำนวณความสูงจีออยด์

surveyor-pocket-tools_2016-12-25_20-09-37

  • นำค่า -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)

soffice-bin_2016-12-25_20-30-14

  • หรือถ้าทราบ ESF แบะ GSF ก็รวบรัดดังนี้

soffice-bin_2016-12-25_20-35-31

  • ช่างรังวัดกรมที่ดินกับช่างรังวัดเอกชน เห็นสูตรนี้แล้ว คงคุ้นๆกัน ESF ก็คือ ค่าคูนสัมประสิทธ์ (C) นั่นเอง ส่วน GSF คือค่าคูนมาตราส่วน(K)

soffice-bin_2016-12-25_20-54-14

ตรวจสอบความถูกต้องของพื้นที่บนระบบพิกัดฉาก 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 ดู

Leave a Reply

Your email address will not be published.