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

ทดสอบข้อมูลค่าพิกัดรูปปิดบนพื้นหลักฐาน Indian 1975

  • ขอพูดเกี่ยวกับ work flow สักนิด เมื่ออ่านไฟล์รูปปิดจากไฟล์  CSV มาแล้ว โปรแกรมจะคำนวณหาจุด centroid หรือจุดศูนย์ถ่วงของพื้นที่ แล้วจะแปลงพิกัดเป็นค่าพิกัดในพื้นหลักฐาน WGS84 ทั้งค่าพิกัดของรูปปิดและจุดศูนย์ถ่วงด้วย จากนั้นโปรแกรมจะสร้างเส้นโครงแผนที่ Lambert Azimuthal Equal Area โดยการใช้จุดศูนย์ถ่วงเป็น latitude of origin, central meridian แล้วเรียกใช้ไลบรารี Proj4 เพื่อทำการแปลงพิกัดไปยังระบบพิกัดฉากของเส้นโครงแผนที่ LAEA สุดท้ายใช้สูตรผูกเชือกรองเท้า ทำการคำนวณหาพื้นที่ จะได้พื้นที่ออกมา แต่ยังเป็นพื้นที่บนทรงรีอยู่
  • ที่หน้าต่างโปรแกรมคลิก “Import” แล้ว browse… เพื่อนำไฟล์ CSV เข้า

compute_area_import2

  • จะใช้ข้อมูลที่ผมเตรียมไว้ให้ อยู่ในโฟลเดอร์ “C:\Users\ชื่อผู้ใช้\AppData\Roaming\Surveyor Pocket Tools\example data” เลือกไฟล์ “boundary2-utm47n-indian1975.csv

surveyor-pocket-tools_2016-12-26_17-06-30

  • เมื่อคลิก “ok” จะเห็น preview ไฟล์นี้ไม่มีชื่อหัวคอลัมน์ เราต้องไประบุทีหลัง

compute_area_import3_noheaders

  • เมื่ออ่านไฟล์เข้าตารางข้อมูล จะเห็นตรงกรอบสีแดงด้านบนเป็น “Col1” ทั้งหมด

compute_area_noheaders

  • ตั้งหัวข้อคอลัมน์ให้ตรงกับ ชื่อจุด Northing Easting

compute_area_set_headers

  • ต่อไปก็ตั้งระบบพิกัดเป็น “Indian 1975 / UTM Zone 47N”

compute_area_set_crs

  • คลิกคำนวณที่ไอคอนเครื่องคิดเลข โปรแกรมจะคำนวณหาพื้นที่ให้ มาลองปักหมุดดูกันครับ บน google maps คลิดที่ทูลส์บาร์ด้านขวา

nvidia-share_2016-12-26_20-27-43

  • ปักหมุดที่ google earth

2016-12-26_20-30-02

การส่งออกไฟล์เป็น Shape file

  • ใครที่ไม่ได้ใช้งานจำพวก GIS ก็ข้ามหัวข้อนี้ไปนะครับ อย่างที่ผมบอกไว้ว่าเส้นโครงแผนที่ Lambert Azimuthal Equal Area รักษาพื้นที่แต่ไม่รักษารูปร่าง แต่ข้อดีคือการจุด origin ไปแปะอยู่ที่จุดศูนย์ถ่วง ทำให้รูปร่างไม่น่าจะเพี้ยนมากนัก ต่อไปคลิกที่ “Export…”

compute_area_export_to_shape_file

  • จะเห็นไดอะล๊อก เลือกปลายทางไฟล์ที่จะเก็บ เลือกรูปแบบเป็น “ESRI Shape file

surveyor-pocket-tools_2016-12-26_21-23-27

  • ผมตั้งชื่อว่า “Boundary2” โปรแกรมจะสร้างไฟล์ให้ มีทั้งหมด 4 ไฟล์คือ  Boundary2.dbf, Boundary2.shx, Boundary.shp และ Boundary.prj ผมจะลองเปิด QGIS แล้วเปิดไฟล์ชุดนี้เข้าไป ที่ QGIS ที่เมนู “Project => New” เพื่อสร้าง project ใหม่ จากนั้นคลิกที่เมนู “Project => Project Proterites” เลือกหน้า CRS (Coordinate Reference System) ตั้งค่าตามรูปด้านล่าง คือเราจะให้ QGIS แปลงพิกัดเป็น WGS84 ในขณะนำเข้า

qgis-bin_2016-12-27_05-36-35

  • แต่หลายครั้งผมพบว่า QGIS รู้สึกจะเอ๋อๆ ไม่ยอมแปลงพิกัดแบบทันทีทันใดตอนนำเข้า จากนั้นเลือกเมนู “Layer => Add Layer => Add Vector Layer…” เลือกไฟล์ชือ “Boundary2.shp”
  • เข้ามาแล้ว ผมตรวจดู QGIS แปลงพิกัดให้เรียบร้อย ไปที่เมนู Setting => Custom CRS… จะเห็นร่องรอยระบบพิกัดที่นำไฟล์เข้าไป QGIS

qgis_user_defined_crs

  • ข้อควรจำ ระบบพิกัดของเส้นโครงแผนที่ LAEA  ที่ผมทำขึ้นมาช่วยหาพื้นที่ ตัวนี้ควรใช้ชั่วคราวเท่านั้นนะครับ เพราะไม่ได้เป็นมาตรฐาน ถ้านำเข้าโปรแกรมด้าน GIS ควรรีบแปลงระบบพิกัดไปหาตัวมาตรฐานอื่นทันที
  • แต่ถ้าผู้อ่านต้องการส่งออกไฟล์ Shape file บนระบบค่าพิกัดเป็น UTM แต่จำกัดเฉพาะบน WGS84 เพราะว่า Indian 1975 ที่อยู่ตามโปรแกรม GIS ทั้งหลายค่าพารามิเตอร์ Transformation ไม่ตรงกับประเทศไทย ที่เราใช้กัน ค่าที่เราใช้กันคือ dx = 204.5, dy = 837.9, dz = 294.8 ตามกรมแผนที่ทหารปี 2551 ดังนั้นถ้ามีการแปลงพิกัดบนโปรแกรมเหล่านี้ไปหาระบบพิกัดอื่นจะไม่ถูกต้อง เว้นเสียแต่ว่าเราสามารถกำหนดตัวแปรพารามิเตอร์เองได้
  • ที่โปรแกรมคลิก “Export…” ด้านขวา เลือกโฟลเดอร์ปลายทางแล้วป้อนขื่อไฟล์ แล้วลองเปิดบนโปรแกรมด้าน GIS ดู

compute_area_export_shapefile_utm

การเปลี่ยนหน่วยพื้นที่

  • ผ่านมาหลายตอนแล้วผมลืมบอกไปว่า สามารถเปลี่ยนหน่วยพื้นที่ได้จาก ตร.ม. ไปยัง หน่วยไร่ของบ้านเรา หรือแม้แต่หน่วย hectare หรือ acre รอบๆบ้านเรายังใชหน่วยพวกนี้อยู่ เช่นเมียนมา ยังใช้หน่วย เอเคอร์อยู่

compute_area_changed_unit

ปัญหาของแปลงที่ดินอยู่คร่อมเส้นแบ่งโซน

  • ปัญหาของระบบพิกัดฉากอีกอย่าง คือตรงบริเวนเส้นแบ่งโซน ที่ศูนย์พิกัดอยู่คนละที่กัน ถ้าเป็นงานก่อสร้างบริเวณช่วงแบ่งโซนนี้ อาจจะทำหมุดไว้อย่างน้อยหนึ่งคู่ พร้อมกับมีพิกัดอิงอยู่ทั้งสองโซน
  • สำหรับแปลงที่ดิน ถ้าไม่ล้ำโซนไปหากันไกลมากนัก ก็น่าจะอนุโลมในการใช้โซนข้างใดข้างหนึ่ง  การคำนวณพื้นที่ถ้าจะต้องมาแบ่งพื้นที่ตามโซนแล้ว นำพื้นที่มารวมกันทีหลัง ค่อนข้างยุ่ง
  • แต่ถ้าแปลงที่ดินเก็บค่าพิกัดแบบภูมิศาสตร์ ก็ง่ายครับ มาดูข้อมูลทดสอบแปลงรูปปิด อยู่บนพื้นที่ระว่างโซน 47 และโซน 48 ค่าพิกัดเป็นแบบภูมิศาสตร์
  • คลิก “Import…” ที่ทูลส์บาร์ด้านขวา เลือกไฟล์ชื่อ “boundary4-crossed-zone47n-zone48n.csv

python_2016-12-27_11-06-49

  • ดู preview ไฟล์มีหัวชื่อคอลัมน์ ค่าพิกัดเป็นแลตติจูด ลองจิจูด

python_2016-12-27_11-07-01

  • ที่ตารางข้อมูล ตั้งระบบพิกัดเป็น “WGS84 / Geographic” คลิกที่ไอคอนเครื่องคิดเลขทำการคำนวณ ที่นี้ผลการคำนวณพื้นที่ในระบบพิกัดฉาก โปรแกรมจะพบว่ามีการข้ามโซน จะคำนวณพื้นที่ให้ทั้งสองโซน ค่าจะต่างกันเพราะว่าค่า scale factor ไม่เท่ากัน พื้นที่แปลงนี้เมื่อคำนวณบนระบบพิกัดฉากโซน 47 จะได้ค่าออกมา = 291 ไร่ 2 งาน 87.47 ตารางวา 

compute_area_crossed_zone_comparison1

  • ลองคลิกเลือกเป็นโซน 48 ตามรูปจะได้พื้นที่ 291 ไร่ 2 งาน 86.61 ตารางวา สองโซนต่างกัน = 87.47 – 86.61 = 0.86 ตารางวา ถือว่าต่างกันน้อยมาก แต่ถ้าเทียบกับพื้นที่บนทรงรีแล้ว ต่างกัน 2 งานกว่าๆ ถือว่ามากพอสมควร

compute_area_crossed_zone_comparison2

ตรวจสอบพื้นที่บนทรงรี

  • ที่จุดศูนย์ถ่วง latitude = 13.9468690 longitude = 102.0021300  ค่าระดับเฉลี่ยเทียบกับรทก. (H) ประมาณ 30 เมตร ความสูงจีออยด์(N) =-26.653 เมตร ความสูงเทียบกับทรงรี (h) = H + N = 30-26.653 = 3.347 เมตร รัศมีทรงรีบริเวณนี้ = 6,376,910 เมตร หาค่า ESF = 6376910/(6376910+3.347) = 0.999999475
  • ค่า บน UTM Zone 47N GSF = 1.00090202 บน UTM Zone 48N GSF = 1.000898320
  • พื้นที่ในระบบพิกัดฉาก UTM Zone 47N = 466749.898 ตร.ม. คิดเป็นพื้นที่บนทรงรี =  466749.898 /  1.00090202² = 465909.000 ตร.ม. พื้นที่บนทรงรีที่ได้จากเส้นโครงแผนที่ LAEA = 465909.007 ตร.ม. ต่างกันที่ทศนิยมตำแหน่งที่ 3 ถือว่าน้อยมาก
  • พื้นที่ในระบบพิกัดฉาก 48N = 466746.458 ตร.ม. คิดเป็นพื้นที่บนทรงรี  = 466746.458/1.000898320 = 465909.011 ตร.ม. ต่างกันที่ทศนิยมตำแหน่งที่ 3 เช่นเดียวกัน อย่างไรก็ตามค่าที่ได้จาก LAEA ผมถือว่าให้ค่าที่ถูกต้องที่สุด
  • คิดเป็นพื้นที่จริง = 465907.007 / 0.999999475² = 465907.496 ตร.ม.
  • ลองปักหมุดดูบน google earth เส้นสีน้ำเงินหนาๆคือเส้นแบ่งโซน ด้านซ้ายคือโซน 47 ด้านขวาคือโซน 48

googleearth_2016-12-27_11-27-45

สรุป

  • การคำนวณพื้นที่บนทรงรีจากเส้นโครงแผนที่ Lambert Azimuthal Equal Area มีความน่าเชื่อถือ ทำให้การคำนวณพื้นที่จากระบบพิกัดฉาก UTM สามารถหาพื้นที่ได้โดยไม่ต้องอาศัยการคำนวณหาค่า Grid Scale Factor
  • ก็หวังว่าทูลส์ตัวนี้จะช่วยช่างสำรวจคำนวณพื้นที่ได้โดยสะดวกและสามารถนำผลลัพธ์ไปใช้งานอื่นๆได้ง่าย

เครดิต

  • ผมได้เพิ่มไลบรารีภาษา python มาช่วยอีกสองไลบรารีคือช่วยในการเขียน Excel และเขียน shape file ทำให้งานยากๆกลายเป็นเรื่องง่ายๆ
  • เขียนไฟล์ Excel ใช้ openpyxl พัฒนาโดย Eric Gazoni, Charlie Clark ใช้งานง่าย มีทุกอย่างที่ต้องการ ไม่ต้องอ่านคู่มือมาก
  • เขียนไฟล์ ESRI Shape file ใช้ pyshp พัฒนาโดย Joel Lawhead ใช้งานง่ายมาก อ่านคูมือไม่กี่บรรทัดก็ใช้งานได้แล้ว

ก้าวเล็กๆต่อไปของ Surveyor Pocket Tools

  • ในภายภาคหน้า ผมจะเขียนทูลส์เล็กๆมาช่วยคำนวณเรื่อง scale factor ทั้ง ESF, GSF และ CSF  ให้มาใช้งานได้ง่ายๆสะดวก อานิสงส์ของโปรแกรมคำนวณพื้นที่ตัวนี้ ทำให้ผมสามารถเอาสูตรแปลงพิกัด Geographic => UTM และแปลงจาก UTM => Geographic และการคำนวณ ESF & GSF ลงเครื่องคิดเลข Casio Fx5800p เนื่องจากติดตามสูตรเขียน scale factor แต่ไปเจอสูตรทั้งหมดอยู่ด้วยกัน ก็เลยเอามาลงที่เครื่องคิดเลขได้ทั้งหมดแบบนึกไม่ถึง
  • ทุกโปรแกรมบนเครื่อง desktop pc & notebook ของผมยังฟรีเหมือนเดิม ร่วมแบ่งปันกัน โลกนี้จะน่าอยู่มากยิ่งขึ้น” ติดตามตอนต่อไปครับ

17 thoughts on “Surveyor Pocket Tools – คำนวณพื้นที่ เรื่องธรรมดาที่ไม่ธรรมดา (ตอนที่ 3 ตอนจบ)”

  1. ขอบคุณสำหรับโปรแกรมครับ รบกวนนิดหนื่งครับ Surveypockettools ตอนนีใน app.box.com ส่วนโปรแกรมอื่นดาวโดหลดได้ปกติครับ

  2. สุดยอดเลยครับ, ผมเป็นอาจารย์สอนเกี่ยวกับการสํารวจอยู่ วิทยาลัยเทคนิค สับพะวิชา (Polytechnic College) ที่ สปป ลาว, ผมชื่นชอบบลอกของคุณมาก…มีสาระ, ความรู้ และ อธิบายเข้าใจง่ายๆ. แต่ถ้าเป็นไปได้ขอให้คุณเพี่มการคำนวณที่ใช้ของลาวเพี่มด้วยจะเป็นพระคุณอย่างสูง, ผมจะเอามาสอนนักเรียนของผม. ลาวใช้ละบบ Lao National Datum 1997.Lao97_UTM_47N, Lao97_UTM 48N. Lao97 to WGS84: dx=46.012 dy=-127.108 dz=-38.131 ex=0 ey=0 ez=0 m=0. Lao97 use Krassovsky 1940. ขอบคุณครับ.

  3. ยินดีมากมากครับ ถ้าโปรแกรมของผมจะเป็นประโยชน์แก่น้องๆหลานๆบ้านใกล้เรื่อนเคียงกัน ผมจะเพิ่ม Lao International Datum 1997 ให้ทั้ง UTM47N และ UTM48N แต่ขอติดวงเล็บ unofficial ไว้ท้ายเพราะไม่เคยใช้ระบบพิกัดนี้ ถ้าอาจารย์มั่นใจ ว่าใช้เป็นทางการแล้ว ช่วยส่งเอกสารยืนยันมาทางอีเมล์ผม riabroy@gmail.com หรือมีลิ๊งค์มายืนยันก็ได้ ผมจะได้ปลด unofficial ออกให้ครับ สำหรับโปรแกรมอัพเดทรออีกสักพัก

  4. ขออนุญาตสอบถามครับ ผมสร้าง LDP ขึ้นใช้เอง แต่เมื่อนำมาใช้กับเมนู File Transform Coordinate ปรากฎว่ามันไม่คำนวณครับ (จาก UTM ไปยัง LDP)

    1. เดี๋ยวจะตรวจสอบให้ครับ ไลบรารี pyproj ที่ใช้งานเป็นแกนหลักงานคำนวณ มีการเปลี่ยนแปลง API บ่อย บางครั้งผมอัพเดทแต่ไม่ได้ทดสอบ ทุกฟังก์ชัน ทำให้ฟังก์ชันบางอย่างไม่ทำงานครับ

      1. ขอบคุณครับ ตอนนั้นผมย้อนกลับไปเวอร์ชันเก่า ใช้ได้ครับ

      2. อีกประเด็นครับผม เวอร์ชันที่ผมใช้อยู่ V1.20 build 663 เมื่อใช้ File transform Coordinates แล้ว save ผลเป็นไฟล์ excel ไม่แน่ใจว่า Northing กับ Easting สลับกันรึเปล่าครับ หรือผมเข้าใจโปรแกรมผิด

        1. สามารถสลับคอลัมน์ได้ที่โปรแกรม ตามรูปนี้ครับ รูปตัวอย่าง

          1. หมายถึงข้อมูลไม่ตรงช่องน่ะครับ ผมทดลองให้แปลง จาก UTM 47N เป็น 48N ดูเหมือนในไฟล์ excel จะสลับ N, E ครับ
            https://ibb.co/g9d4DLG

Leave a Reply

Your email address will not be published.