เป็นความฝันของผมอันหนึ่งตั้งแต่สมัยจบใหม่ๆที่จะเขียนโปรแกรมสร้างชั้นความสูงจากข้อมูลจุดงานสำรวจ x, y, z แต่จนแล้วจนเล่าโครงการนี้ไม่เคยเกิดสักที เนื่องจากความรู้ความสามารถและทักษะไม่เพียงพอ ด้วยความยากในการคิดอัลกอริทึมที่จะสร้างสามเหลี่ยมด้านที่สั้นที่สุดจากจุด (point) งานสำรวจ จนกระทั่งเลยวัยแห่งความฝันอันนั้นมาไกลมากแล้ว
ปัจจุบันในยุค open source มีไลบรารีด้านนี้ที่มีโมดูลส่วนหนึ่งที่มีความสามารถใกล้เคียงที่สามารถนำมาสร้างเส้นชั้นความสูงได้คือ matplotlib
ไลบรารี matplotlib
ตัว matplotlib เองถูกนำไปใช้ร่วมกับ numpy, pandas สุดยอดของไลบรารีไพทอนที่นำไปประมวลผลงานด้าน Data Science & Data Analysis ตลอดจน Machine Learning ทะลุไปถึง AI & Big Data จนเกิดทำให้กระแสความนิยมภาษาไพทอนมาแซงภาษาจาวาที่นำโด่งก่อนหน้านี้มาหลายสิบปี matplotlib ถูกนำไปใช้ในงาน Data Visualization ที่นำเสนองานคำนวณต่างๆในรูปกราฟฟิคในลักษณะกราฟต่างๆ ได้สะอาด สวยงาม ถูกต้อง
Tricontour Smooth Delaunay
โค้ดตัวอย่างอยู่ในเว็บไซต์ของ matplotlib เรียกวิธีการนี้ว่า “Tricontour Smooth Delaunay” เป็นสร้างเส้นชั้นความสูงที่มีโค้งมนสวยงาม ผมดูแล้วเทียบได้กับการลากเส้นด้วยมือในสมัยก่อนทีเดียว รายละเขียดขั้นตอนและวิธีการ เท่าที่ผมอ่านและประมวลได้คือ
- ใช้ฟังก์ชัน Triangulation สร้างสามเหลี่ยมด้านสั้นที่สุดทั้งหมดโดยอาศัยจุดสำรวจ x, y, z ในที่นี้จะรวมสามเหลี่ยมที่รูปทรงยาวลีบตรงขอบงาน (ถ้าผู้อ่านใช้โปรแกรมเชิงพานิชย์เช่น Autocad Civil3D หรือ Trimble Business Center จะคุ้นกับการแก้ไขลบรูปสามเหลี่ยมทรงลีบๆยาวๆตรงขอบงานสำรวจ) ถ้าไม่มีการแก้ไขตรงนี้จะทำให้เส้นชั้นความสูงเพี้ยน เพราะจุดอยู่ไกลกันมากเกินไป
- ใช้ฟังก์ชัน TriAnalyzer เพื่อที่จะกรองสามเหลี่ยมที่ไม่ดีออกได้แก่ สามเหลี่ยมที่มีความสูงเท่ากันทั้งสามด้าน (flat triangles) รวมถึงสามเหลี่ยมที่มีทรงยาวลีบตรงขอบงาน
- ใช้ฟังก์ชัน UniformTriRefiner ในการประมาณค่า (interpolate) โดยการเพิ่มจำนวนสามเหลี่ยมให้มากขึ้น (recursive subdivisions) ทำให้เส้นชั้นความสูงมีความโค้งมนเนียนมากขึ้น
- สร้างเส้นชั้นความสูงจากสามเหลี่ยมด้วยฟังก์ชัน tricontour
รายละเอียดข้อมูลงานสำรวจสำหรับทดสอบ
ส่วนจุดสำรวจ x, y, z ได้จากงานสำรวจของอ่าวพัทยากลาง ที่สำรวจตั้งขอบฟุตบาททางเดินข้างถนน ลงมาที่ชายหาดและก็ลงน้ำไป ใช้เครื่องมือ RTK ในการสำรวจงานบก ส่วนในน้ำได้ข้อมูลจากการใช้เรือสำรวจติด GNSS และ Echo Sounder หรือเครื่องหยั่งน้ำ นำข้อมูลสองส่วนนี้มารวมกัน พื้นที่ในการสำรวจหน้ากว้างตามชายหาดประมาณ 3 กม. ลึกลงไปในทะเลประมาณ 2.5 กม. ไฟล์ที่เก็บข้อมูลชือ “pattaya_topo_hydro_utm47n_wgs84.csv” มีจำนวนจุดทั้งหมด 4704 จุด
โค้ดไพทอน
โค้ดต้นฉบับนั้นใช้การ random จุดขึ้นมาในพื้นที่ 2 x 2 ถ้าผมคิดหน่วยเป็นเมตร ก็ได้ 4 ตารางเมตร ซึ่งงานสำรวจจริงๆจะใหญ่กว่านั้นมาก อย่างเช่นจุดทดสอบที่ผมกล่าวมาข้างต้นประมาณ 7.5 ตร.กม. ดังนั้นผมจึงยำโค้ดต้นฉบับให้สามารถอ่านไฟล์จากงานสำรวจ และตัดทอนบางอย่างออกไป เพิ่มบางอย่างเข้ามาเช่น ให้เลือกพล็อทรูปสามเหลี่ยมที่ถูกแบ่งย่อยออกมาให้เห็น รูปสามเหลี่ยมเมื่อแบ่งแล้วเป็นยังไง
นี่เป็นโค้ดไพทอน ตั้งชื่อไฟล์ “tricontour.py” สามารถคัดลอกและจัดเก็บไฟล์นี้ได้ตามสะดวก สังเกตว่าโค้ดมีไม่กี่สิบบรรทัด นั่นหมายความว่าไลบรารี matplotlib เอางานไปคำนวณเบื้องหลังฉาก
import os
import matplotlib
from matplotlib.tri import Triangulation, TriAnalyzer, UniformTriRefiner
import matplotlib.pyplot as plt
from matplotlib.ticker import FuncFormatter
import matplotlib.cm as cm #โทนสี color map
import numpy as np
# ตั้งค่าสภาพแวดล้อมให้เรียกใช้กราฟฟิค GUI ของ Qt5 เป็น backend
matplotlib.use('Qt5Agg')
matplotlib.rcParams['backend']='Qt5Agg'
# ใส่ 1000 ในกรณีมีจำนวนเส้นและจุดจำนวนมาก รันในวินโดส์ถ้าไม่ใส่โปรแกรมจะแฮ้งค์
plt.rcParams['agg.path.chunksize'] = 1000
subdiv = 3 # ตัวเลขปรับความละเอียด โดยสร้างรูปสามเหลี่ยมเพิ่ม จำนวน = 4**subdiv (ยกกำลัง)
# ถ้าใช้เลข 3 จำนวนสามเหลี่ยม = 4x4x4 = 64 รูป
min_circle_ratio = .01 # ใช้เลข -1 ถ้าต้องการสามเหลี่ยมทั้งหมด แนะนำเลข 0.01
# ดึงข้อมูลงานสำรวจตัวอย่างจากเว็บไซต์ของผม
data = np.genfromtxt(r'https://www.priabroy.name/shared/misc/pattaya_topo_hydro_utm47n_wgs84.csv', delimiter=',')
xs = np.asarray(data[:,2]) # คอลัมน์ที่สาม เป็นค่า Easting
ys = np.asarray(data[:,1]) # คอลัมนที่์สอง เป็นค่า Northing
zs = np.asarray(data[:,3]) # คอลัมน์ที่สี่ เป็นค่าความสูง
print('Data size: {0:d} rows x {1:d} columns'.format(data.shape[0], data.shape[1]))
print('Number of points: ', xs.shape[0])
# ฟังก์ชัน Triangulation สร้างรูปสามเหลี่ยม
tri = Triangulation(xs, ys)
ntri = tri.triangles.shape[0] # จำนวนรูปสามเหลี่ยมพื้นฐานทั้งหมด
print("Number of triangles: ", ntri)
print("Number of subdivision triangles: {0:d}".format(4**subdiv*ntri))
# กรองสามเหลี่ยมที่มีความสูงเท่ากันทัั้งสามจุด (flat triangle)
# กรองสามเหลี่ยมรูปทรงยาวลีบตรงขอบงาน (สามเหลี่ยมไม่ดี)
mask = TriAnalyzer(tri).get_flat_tri_mask(min_circle_ratio)
tri.set_mask(mask)
# เพิ่มจำนวนสามเหลี่ยมให้มากขึ้น โดยการ interpolate แบบ recursive subdivisions
refiner = UniformTriRefiner(tri)
tri_refi, z_test_refi = refiner.refine_field(zs, subdiv=subdiv)
flat_tri = Triangulation(xs, ys) # ดึงรูปสามเหลี่ยมรูปทรงไม่ดีออกมาเพื่อพล็อท
flat_tri.set_mask(~mask)
plot_tri = True # ต้องการแสดงรูปสามเหลี่ยมพื้นฐาน
plot_tri_refi = False # แก้เป็น True ถ้าต้องการเห็นสามเหลี่ยมที่ interpolate แต่กินแรงเครื่องคอมพิวเตอร์
plot_masked_tri = True # ต้องการแสดงสามเหลี่ยมทรงยาวลีบตรงขอบงาน
levels = np.arange(-25, 12.0, 1) # แสดงเส้นชั้นความสูงตั้งแต่ -25m ถึง 12m interval 1m
cmap = cm.get_cmap(name='Blues', lut=None) # ไล่โทนสีเส้นชั้นความสูงในย่านสีน้ำเงิน
fig = plt.figure() # สร้างตัวแปรเกี่ยวกับรูปทรงเพื่อเขียนหัวข้อ เขียนไตเติ้ล
ax = plt.axes() # สร้างตัวแปรที่เกี่ยวข้องกับการพล๊อทในแกน X และแกน Y
ax.set_aspect('equal') # ให้มาตราส่วนแกน X และ แกน Y เท่ากัน
fig.canvas.set_window_title('สร้างเส้นชั้นความสูงแบบความละเอียดสูงด้วยไลบรารี matplotlib')
fig.suptitle("Pattaya Beach\nHigh-resolution contour") #เขียนไตเติ้ล ด้านบนสุด
# จัดรูปแบบให้เขียนค่า Northing และ Easting โดยใช้เครื่องหมายคอมม่าคั่นหลักพัน
ax.yaxis.set_major_formatter(FuncFormatter(lambda x, p: format(int(x), ',')))
ax.xaxis.set_major_formatter(FuncFormatter(lambda x, p: format(int(x), ',')))
#กำหนดขนาดตัวอักษรที่เขียนค่าโคออดิเนทตามแกน X และแกน Y
ax.tick_params(axis = 'both', which = 'major', labelsize = 8)
# เขียนตัวอักษร Northing ด้านซ้ายแกน Y และ Easting ใต้แกน X
plt.xlabel('Easting', fontsize=10, color = 'black')
plt.ylabel('Northing', fontsize=10, color = 'black')
#เขียนเส้นชั้นความสูงแบบละเอียด อาศัยพารามิเตอร์ที่สร้างไว้ก่อนแล้ว
ax.tricontour(tri_refi, z_test_refi, levels=levels, cmap=cmap,
linewidths=[2.0, 0.5, 1.0, 0.5])
if plot_tri_refi: #เขียนรูปสามเหลี่ยมที่แบ่งแล้ว
ax.triplot(tri_refi, color='green', linewidth = 0.075)
if plot_tri: # เขียนรูปสามเหลี่ยมพื้นฐาน
ax.triplot(tri, color='0.97')
if plot_masked_tri: # เขียนสามเหลี่ยมรูปทรงไม่ดี
ax.triplot(flat_tri, color='red')
plt.grid() # สร้างเส้นกริด
plt.show() # เอางานพล็อททั้งหมดแสดงที่หน้าต่างจะไม่อธิบายโค้ดมาก เพราะในโปรแกรมผมได้เขียน comment ไว้หมดแล้ว จะพูดเฉพาะที่สำคัญ ไลบรารี numpy เป็นไลบรารีที่เรานำมาใช้คำนวณพวก Data Science & Data Analysis ที่มีประสิทธิภาพมาก สามารถจัดข้อมูลอยู่ในรูปเมตริกซ์หรืออนุกรม (array) ผมใช้ numpy มาอ่านไฟล์ csv ซึ่งสามารถอ่านได้เร็วมาก ตอนแรกเขียนโค้ดอ่านเอง แต่เมื่อลองใช้ numpy มาอ่าน ความเร็วคนละเรื่อง มันเร็วกว่าหลายสิบเท่า ผมใช้ฟังก์ชัน genfromtxt() อ่านข้อมูลจากเว็บไซต์ผมโดยตรง เข้ามาเก็บในเมตริกซ์ชื่อ data
data = np.genfromtxt(r'https://www.priabroy.name/shared/misc/pattaya_topo_hydro_utm47n_wgs84.csv', delimiter=',')หรือถ้าดาวน์โหลดไฟล์นี้มาแล้ว ก็แก้ไขโค้ดอ่านจากเครื่องได้
data = np.genfromtxt(r'pattaya_topo_hydro_utm47n_wgs84.csv', delimiter=',')รูปแบบข้อมูล csv จะเป็น 5 คอลัมน์ประกอบไปด้วย Point No, Northing, Easting, Elevation, Code จากนั้นทำการเฉือนแบ่งข้อมูลออกมาเป็น x, y, z อย่างเดียว เพราะไลบรารี matplotlib ต้องการแบบนั้น คือแบ่งออกมาเป็นคอลัมน์ใครคอลัมน์มัน จัดเก็บใน array ผมตั้งชื่อ xs, ys, zs (เติม s เข้าไปให้รู้ว่ามีหลายตัว)
จากนั้นจะนำ xs, ys, zs ไปผ่านกระบวนการสร้างรูปสามเหลี่ยม จากนั้นนำไปกรองเอาเฉพาะสามเหลี่ยมที่ดี (สามเหลี่ยมรูปทรงแบบยาวลีบๆตรงขอบงานจะถูกตัดออก) สุดท้ายก็สร้างเส้นคอนทัวร์และนำออกมาแสดงผล ก่อนจะรันโปรแกรมนี้ ต้องติดตั้งไลบรารี numpy และ matplotlib ให้เรียบร้อยก่อน ถ้ายังไม่ได้ติดตั้ง โดยใช้ pip ในการติดตั้งดังนี้
pip install numpy
pip install matplotlibทดสอบโปรแกรม
นำไฟล์ “tricontour.py” มาจัดไว้ในโฟลเดอร์เดียวกับไฟล์ข้อมูลจุดงานสำรวจ จากนั้นใช้ command prompt ในวินโดส์ เข้าไปในโฟลเดอร์นี้แล้วพิมพ์คำสั่งเพื่อเรียกไพทอนมาคอมไพล์และรัน
python tricontour.pyถ้าไม่มีอะไรผิดพลาดจะได้ผลลัพธ์ออกมาดังนี้
ซึ่งเราให้ปรินท์ขนาดจำนวนจุดข้อมูล 4704 จุด รูปสามเหลี่ยมพื้นฐานที่สร้างจำนวน 9342 รูป รูปสามเหลี่ยมรูปเล็กๆที่แบ่งโดยการ interpolate แบบ subdivisions จำนวน 597888 รูป ผมค่อยมาขยายรายละเอียดอีกทีตรงจุดนี้
ซูมเข้าไปเพื่อไปดูรายละเอียด ผู้อ่านจะเห็นรูปสามเหลี่ยมสีเทา นั่นคือรูปสามเหลี่ยมพื้นฐานที่ได้จากการทำ triangulation จากไลบรารี matplotlib
ลองซูมเข้าไปอีกระดับ ก็จะพอเห็นรอยหยักหักของเส้นชั้นความสูงบ้างเล็กน้อย ซึ่งสามารถยอมรับได้
รูปสามเหลี่ยมที่แบ่งย่อย
ผมจะให้ดูรูปสามเหลี่ยมที่แบ่งย่อยเป็นรูปเล็กจำนวน 64 รูป โดยแก้ไขโค้ดให้เห็นรูปสามเหลี่ยมที่แบ่งย่อยแล้ว และจะเขียนรูปสามเหลี่ยมพื้นฐานด้วยสีน้ำเงิน สุดท้ายปิดเส้นกริดเพื่อให้ดูรูปง่าย
plot_tri_refi = True # แก้เป็น True ถ้าต้องการเห็นสามเหลี่ยมที่ interpolate แต่กินแรงเครื่องคอมพิวเตอร์if plot_tri: # เขียนรูปสามเหลี่ยมพื้นฐาน
ax.triplot(tri, color='b')#plt.grid() # สร้างเส้นกริดจะได้ผลลัพธ์ออกมา ถ้าแอบไปดู Task Manager ของวินโดส์ จะเห็นว่าโปรแกรมเขมือบเมมโมรีไป 1.5 GB ผมลองซูมไปที่เก่า ลองเปรียบรูปด้านบน จะเห็นรูปสามเหลี่ยมพื้นฐานสีน้ำเงินที่แบ่งย่อยแล้วเป็นสีเขียวอ่อนๆ จำนวน 64 รูป อยู่ภายในแต่ละรูปสามเหลี่ยมพื้นฐาน จึงทำให้โปรแกรมสามารถสร้างเส้นชั้นความสูงได้โค้งมนเนียนมากยิ่งขึ้นจากสามเหลี่ยมรูปเล็กที่แบ่งย่อยมานี้เอง
แล้วถ้าต้องการนำไปใช้งานจะทำยังไง นี่เป็นปัญหาใหญ่เพราะไลบรารี matplotlib ไม่ได้เตรียมฟังก์ชันให้เขียนแปลงรูปแบบเป็นพวก DXF ผมเลยลองเซฟออกมาเป็นไฟล์ SVG ซึ่งเป็นเวคเตอร์ โดยผมปิดพวกรูปสามเหลี่ยมให้หมด เพราะไม่ได้ใช้งาน เอาแค่เส้นกริดกับเส้นชั้นความสูงไปแค่นั้น
เมื่อได้ไฟล์แล้วผมเอาไฟล์ SVG นี้ไปอัพโหลดเข้าเว็บออนไลน์ที่แปลงไฟล์ SVG => DXF (ลองค้นดูนะครับ ไม่อยากจะระบุจะเป็นการโฆษณา) จะได้ไฟล์ DXF แล้วมาเปิดด้วย CAD ความเป็นพิกัดโลกจะหาย ไม่เป็นไร เรามีเส้นกริดสามารถ Align ให้เข้าตำแหน่งพิกัดโลกได้ดังเดิม จากนั้นสามารถไปใช้งานได้ในแคด ลบพวกตัวหนังสือออกให้หมดเพราะแตกยุ่ยหมดแล้ว
นำเข้า Google Earth
กรณีต้องการนำเข้า Google Earth ผมก็หาเว็บออนไลน์แปลงอีกจาก DXF => KML อัพไฟล์ DXF ที่เซฟจาก CAD แล้วแปลง อาจจะต้องตั้ง UTM Zone ให้ด้วยเป็น UTM Zone 47N เวลาแปลงเป็น KML พิกัดจะได้ไม่หลุดไปอยู่ประเทศอื่น เมื่อแปลงแล้วสามารถนำไฟล์ KML ไปเปิดได้ใน Google Earth แต่ผมพลาดท่าหน่อยตรงสีเส้นชั้นความสูงดันไปกำหนดให้อยู่ในโทนสีน้ำเงิน เวลาไปเจอสีน้ำทะเลแล้วจะดูยากหน่อย
แนวทางพัฒนาต่อยอดโปรแกรม
คงเป็นเพียงแนวคิดแค่นั้นครับ โปรแกรมสามารถนำไปพัฒนาให้มีการติดต่อผู้ใช้โดยตรง สามารถเปิดไฟล์จุดงานสำรวจได้ สามารถตั้งค่าต่างๆได้ กำหนดสีเส้นชั้นความสูง interval แต่การที่จะให้โปรแกรมเขียนเลขหรือตัวอักษรเส้นชั้นความสูงออกมาแสดงผล คงเป็นงานหนักน่าดู ที่อาจจะต้องไปศึกษาโค้ดของ matplotlib ในเชิงลึก ให้สามารถเขียนออกมาเป็นไฟล์ DXF ได้โดยตรง หรือสามารถนำเข้าเส้น Alignment อย่างง่ายๆเข้าไปแล้วให้โปรแกรมสร้างเส้นแนวตั้งฉากเพื่อตัด cross-section ก็ทำได้ เนื่องจากรูปสามเหลี่ยมที่สร้างโดย matplotlib แต่ละรูปมีค่าพิกัดและค่าความสูงแต่ละจุดกำกับไว้อยู่
เอาละครับจากโปรแกรมโค้ดไม่กี่สิบบรรทัด ตอนนี้มาไกลมาก โปรดติดตามกันตอนต่อไปครับ
เป๊ะเวอร์
ขอบคุณครับ