QGIS航拍地形模型整理

發表於 分類於 閱讀次數：

專案背景

經由航拍測量後的DEM資料經常因地表植物以及行水中等問題導致經WebODM產製出的地形模型失真，如欲進行後續的地形剖面(橫斷面、縱斷面)時會有需要手動修正的困擾。

本專案當初採用全測站進行收測，但年代已久地形已有部分變動，僅可認定在結構體及山坡地等未開發區域的高程為正確值，其餘部分則先僱工進行雜樹木清除，以航拍測量套疊正射影像後的屬於裸露地的高程視為正確值，結合上述兩項測量成果進行重繪以取得較為精準的地形模型。

重繪方法

經航拍測量後所得到的結果如下:

  • dtm.tif
  • odm_orthophoto.tif

匯入GIS後將裸露地部分的範圍做一個範圍，並且用這個範圍作為遮罩對原始dtm做修剪，之後將dtm的像素格放大，取得像素格中心點進行dem取樣，做為不規則三角網的高程內插來源。

前人測量資料如為結構體部分則須將該線段的頂點及線段加繪內插後的點資料繪出成CSV檔，之後匯入至GIS中做為不規則三角網的高程內插來源，最後與裸露地部分一併進行不規則三角網的高程內插(TIN Interpolation)即可完成，作為後續剖面用的地形模型。

操作步驟

QGIS

  1. 將航照測量產製出的dtm.tif匯入至QGIS
  2. 將航照測量產製出的odm_orthophoto.tif匯入至QGIS
  3. 繪製要留存的測量範圍(裸露地部分)，得到AREA.shp
  4. 以AREA.shp為遮罩範圍，使用工具 Clip raster by mask layer針對dtm.tif進行初剪，得dtm_area.tif。
  5. Raster>Projection>Wrap(reproject)...將dtm_area.tif處理成2M*2M的像素格，得到dtm_area_2m.tif。
  6. 透過工具 Raster pixels to points將dtm_area_2m.tif各像素格生成中心點資料，得到dtm_area_2m_points.shp
  7. 透過工具 Sample raster value將dtm_area_2m_points.shp針對dtm.tif進行採樣，得到 dtm_area_2m_z.shp，帶有高程資料。

目前所得到的航拍高程資料已經初步完成，重投影的部分是將原本過於密集的像素格放大，後續要進行手動編修時才不會那麼痛苦。

CAD部分

Fig1. 原測量成果圖
  1. 將平面圖上既有點資料收進Excel
  2. 確保線段各頂點高程皆可對應到既有點資料
  3. 將線段複製貼上原位置到新檔案
Fig2. VBA控制畫面
  1. 點選按鈕 "頂點加密"
Fig3. 頂點加密成果圖
  1. 點選按鈕 "匯入點位"
  2. 點選按鈕 "輸出CSV"
Fig4. 渠道模型點資料

CSV檔案可以直接匯入QGIS，如果不夠密則容易導致模型失真，不規則三角網會去抓最鄰近的點資料作為高程內插來源，密一點模型比較不會抓錯。

QGIS

經上述兩個步驟後，點資料初步整理內容如下，裸露地部分參考航照測量回傳內容，渠道部分參考平面測量成果，進行TIN Interpolation，將上述兩個資料來源一併計入高程內插即可得到一較不失真的高程模型了!

Fig5. 初步整理後QGIS點資料
Fig6. TIN Interpolation畫面
Fig7. 整理後地形模型成果

頂點加密原始碼

main module

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
'TODO:
'1.因應地形點在QGIS建模時如不夠密的話會導致建模內容出現不合理現象
'  建議從原線段之前後頂點將點資料內插至該線段。
'2.UAV掃描內容如遇水面或植栽容易出現不合理高程，需後續修整不合理點後重新建立
'====

Sub Polyline_Interpolate_ZPoints(ByVal ent, Optional ByVal interval As Integer = 20)

    Dim pt As New clsPt
    Dim PL_Obj As New clsPL
    Dim x1 As Double
    Dim y1 As Double
    Dim x2 As Double
    Dim y2 As Double

    Call PL_Obj.getPropertiesByPL(ent)

    coords = PL_Obj.getPTS

    ' 2D Polyline Coordinates 格式：x0,y0,x1,y1,x2,y2...
    For i = LBound(coords) To UBound(coords) - 3 Step 2
        
        x1 = coords(i)
        y1 = coords(i + 1)
        x2 = coords(i + 2)
        y2 = coords(i + 3)
        
        z1 = FindZByEN(x1, y1)
        z2 = FindZByEN(x2, y2)
        
        Call pt.getPropertiesByUser("ORI", x1, y1, z1, "ORI")
        Call pt.CreatePoint(0.5)
        
        segLen = Distance2D(x1, y1, x2, y2)
        
        d = interval
        
        Do While d < segLen
            
            ratio = d / segLen
            
            newX = x1 + (x2 - x1) * ratio
            newY = y1 + (y2 - y1) * ratio
            newZ = Round(z1 + (z2 - z1) * ratio, 3)
            
            Call pt.getPropertiesByUser("ADD", newX, newY, newZ, "ADD")
            Call pt.CreatePoint(0.5)
            
            outRow = outRow + 1
            d = d + interval
            
        Loop
        
    Next i
    
    Call pt.getPropertiesByUser("ORI", x2, y2, z2, "ORI")
    Call pt.CreatePoint(0.5)
    
    Debug.Print "內插完成，共產生 " & outRow - 2 & " 個新點"

End Sub

Function Distance2D(x1 As Double, y1 As Double, x2 As Double, y2 As Double) As Double
    Distance2D = Sqr((x2 - x1) ^ 2 + (y2 - y1) ^ 2)
End Function

Function FindZByEN(E As Double, N As Double) As Double

    Dim ws As Worksheet
    Dim lastRow As Long
    Dim r As Long
    Dim tol As Double
    
    tol = 0.001 ' 座標容許誤差
    
    Set ws = ThisWorkbook.Worksheets("FIX_POINT")
    lastRow = ws.Cells(ws.Rows.Count, 1).End(xlUp).Row
    
    For r = 2 To lastRow
        
        If Abs(ws.Cells(r, 2).Value - E) <= tol _
        And Abs(ws.Cells(r, 3).Value - N) <= tol Then
            
            FindZByEN = ws.Cells(r, 4).Value
            Exit Function
            
        End If
        
    Next r
    
    Err.Raise vbObjectError + 1000, , "找不到座標對應的 Z：" & E & "," & N

End Function

clsPt

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

Private PT_NUM As String
Private E As Double
Private N As Double
Private Z As String
Private CD As String

Private CAD As New clsACAD
Private sht_PT As Object

Const NUM_DIGIT As Integer = 4

Sub getPropertiesByUser(ByVal bPT_NUM As String, ByVal bE As Double, ByVal bN As Double, ByVal bZ As Double, ByVal bCD As String)

PT_NUM = bPT_NUM
E = bE
N = bN
Z = bZ
CD = bCD

End Sub

Sub CreatePoint(ByVal txtheight As Double)

If txtheight = 0 Then txtheight = getTxtHeight

Dim txtpt(2) As Double

txtpt(0) = E: txtpt(1) = N

'coll_attr.Add PT_NUM
'coll_attr.Add Z
'coll_attr.Add CD

arr = Array(PT_NUM, Z, CD)

Set ptObj = CAD.InsertBlock(txtpt, "Point", txtheight, arr) ' coll_attr)

Call CAD.setLayer(ptObj, "平面圖-點")

'ptObj.Layer = "平面圖-點"

End Sub

clsPL

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
Private CAD As New clsACAD
Private Math As New clsMath
Private pt As New clsPt
Private myfunc As New clsMyfunction

Public PL As Object
Public PLname As String '不是那麼重要
Private pts As Variant
Private co As Integer
Private totalLen As Double

Private sht_PL As Object

Const NUM_DIGIT As Integer = 4

Sub getPropertiesByPL(bPL) 'get the PL object

Set PL = bPL
pts = CAD.tranIPoints(PL.Coordinates)
totalLen = PL.Length
co = 3: If TypeName(PL) Like "*LWPolyline" Then co = 2
PLname = PL.Layer

End Sub

Function getPTS()

Dim coords_new() As Double
Dim idx As Long

idx = -1

For j = 0 To UBound(pts) Step co

    x1 = pts(j)
    y1 = pts(j + 1)

    idx = idx + 2
    ReDim Preserve coords_new(idx)

    coords_new(idx - 1) = x1
    coords_new(idx) = y1

Next

getPTS = coords_new

End Function