2018-11-01

程式解析公開地形模型資料測試

目前台灣政府開放了相當多的資料上網,
其中有一組資料是台灣附近20公尺網格數值地形資料,
https://data.gov.tw/dataset/35430
也就是工程上常說的DEM或DTM,
當然,20公尺的解析度可能無法拿來做工程用途,
但對資料工程師來說,
資料處理的問題都只有「量」與「質」的差異,
需要精度更高的特定範圍資料,
應該要由專業的測量公司協助提供。

另外一組新版資料 https://data.gov.tw/dataset/88733
在筆者撰寫此文章時測試完畢後部份花蓮資料有缺失,
相信不久後應該會被修正,
所以本文還是以上一組資料進行說明。

本文簡單的說明一下如何用程式技巧來處理 大約2.5G的台灣地形資料,
為了方便展示,筆者使用了Python跟C#來做立體圖表及平面的處理。

首先,
從上述政府開放資料平台下載該資料
解壓縮後有12,508個檔案 19個資料夾 (本測試並沒有納入澎湖縣)
打開內部的檔案來看一下,
除了有資料來源說明之外,
另外就是每個點的二度分帶座標及高程。

如果眼尖一點會發現有一些 -999 的存在,
估計可能是海面等無法量測的區域,
後續處理的時候就直接排除或處理成透明資料。

首先我們先用簡單的Python工具去將這個圖表以Mesh3D的方式畫出來,
工具非常的多,這邊我挑一個容易展示的Plotly來使用,
這裡我以桃園的某個區域進行測試,
範圍約2.5km * 2.5km,圖表不照實際比例。


關鍵Python程式碼如下:

import plotly.plotly as py
import plotly.graph_objs as go
xdata=[]
ydata=[]
zdata=[]
dem = open('E:/TY_DEM/96221012dem.grd','rt')
for x in dem:
    point = x.split()
    xdata.append(float(point[0]))
    ydata.append(float(point[1]))
    zdata.append(float(point[2]))
trace = go.Mesh3d(x=xdata,y=ydata,z=zdata,color='#FFB6C1',opacity=0.50)
fig = go.Figure(data=[trace])
py.plot(fig)

當然,相同的方法也可以套用在整個桃園區:

甚至是整個台灣:

只是點位數量有所限制,
需要做一點前處理,
但因為這個圖表工具本身並非處理地形使用,
在比例上並不正確,
資料處理後的精度也並不足以表現地貌,
這裡僅是用來展示大致的資料內容。

但很明顯的,
圖表工具並不適合來展示複雜地形資料,
單純的平面地圖配合顏色反而是很適合的工具。

先來畫一下整個桃園區:
藍色大約是高程 0~750公尺,
綠色大約是高程 750~1500公尺,
紅色則是高程超過 1500公尺。

進一步我們來看一下整個台灣畫出來的結果:
因為高程變化更大了,
這裡分的更細一點,
強化平原的細節,弱化了高山的部份,
藍色大約是: 0~ 250 m
黃色大約是: 250 ~ 750 m
紅色大約是: 750 ~ 1750 m
青色大約是: > 1750 m

C#關鍵程式碼如下:

Bitmap bmp = new Bitmap(width, height);
string[] filePaths = Directory.GetFiles(@"E:\DEM\", "*.grd", SearchOption.AllDirectories);
foreach (string file in filePaths)
{
    using (StreamReader sr = new StreamReader(file))
    {
        string tmp;
        while ((tmp = sr.ReadLine()) != null)
        {
            string[] point = tmp.Split(' ');
            int xx = Convert.ToInt32((Convert.ToInt32(point[0]) - x_offset) / scale);
            int yy = height - (Convert.ToInt32((Convert.ToInt32(point[1]) - y_offset) / scale));
            int zz = Convert.ToInt32(Convert.ToSingle(point[2]));
            if (xx < 0 || yy < 0 || xx > width || yy > height || zz<0)
                continue;
            if (zz < 255)
                bmp.SetPixel(xx, yy, Color.FromArgb(255 - zz, 255 - zz, 255));
            else
            {
                zz = (int)(zz - 255)/2;
                if (zz < 255)
                    bmp.SetPixel(xx, yy, Color.FromArgb(zz, zz, 255 - zz));
                else
                {
                    zz = (int)(zz - 255) / 2;
                    if (zz < 255)
                        bmp.SetPixel(xx, yy, Color.FromArgb(255, 255-zz , 0));
                    else
                    {
                        zz = (int)(zz - 255) / 2;
                        if (zz > 255)
                            zz = 255;
                        bmp.SetPixel(xx, yy, Color.FromArgb(255 -zz , zz, zz));
                    }                                
                }
            }
        }
    }
}
bmp.Save("TW.bmp");

以本案例來說,
成果的檔案約113MB,
差不多可以放大到這種程度。



後續我們也測試了很多相關的進階應用,
例如把所有資料存入資料庫,
處理成Mesh、計算挖填方體積等,
從工程需求的角度來處理地形資料,
但就比較不適合套用在這個資料集進行展示。

最後,感謝政府資料開放平台 及 內政部地政司 提供資料,
讓我們可以在合理的授權下對資料進行加值應用及展示。