iT邦幫忙

2019 iT 邦幫忙鐵人賽

DAY 28
0

DEM(Digital Elevation Model)以二維網格的方式呈現地形,也是GIS資料常見得網格式資料,我們今天將會使用內政資料開放平臺的20米全台DEM資料,他的副檔名是grd格式(ascii),我們可以使用Rasterio讀取資料,並做一些處理。(ps.如果要使用全球的DEM,可以使用SRTM)

大綱:

  • 網格檔案GRD操作
  • 暈渲圖(Hillshading)
  • 坡度坡向

網格檔案GRD操作

grd格式可以使用GDAL中讀取,當然rasterio也可以直接讀取

import rasterio
src=rasterio.open('data/95191010dem.grd')

import matplotlib.pyplot as pltfrom 
from rasterio import plot
%matplotlib inline
plot.show(src)

https://ithelp.ithome.com.tw/upload/images/20181112/20107816K4qFZL49ZR.png

在資料處理方面,可以使用numpy以array方式處理

import numpy as np
height=src.read(1)  # 通常就是一個band
height

https://ithelp.ithome.com.tw/upload/images/20181112/20107816c3lSI4lY3A.png

如果要轉格式,也可以把array另存成tif等格式
tif,png等等都可以是GIS網格資料的格式
(配合wordfile請參考[Day 14] WebGIS中的網格資料)

with rasterio.open('output/dem.tif', 'w', driver='GTiff', height=height.shape[0], width=height.shape[1], count=1, dtype=height.dtype) as dst:
    dst.write(height, 1)

暈渲圖(Hillshading)

暈渲圖是地形資料常見的呈現方式,主要概念是將光源投射在地形上的時候,將光照區及陰影區一同呈現在圖上,常見於地質極大地測量等應用。

暈渲圖的計算,是假設源在某個方向和某個太陽高度下,賦予每個像元其周遭的像元的值計算而成。

def hillshade(array, azimuth, angle_altitude):

    # 取自: http://geoexamples.blogspot.com.br/2014/03/shaded-relief-images-using-gdal-python.html

    x, y = np.gradient(array)
    slope = np.pi/2. - np.arctan(np.sqrt(x*x + y*y))
    aspect = np.arctan2(-x, y)
    azimuth = azimuth*np.pi / 180.
    altitude = angle_altitude*np.pi / 180.


    shade = np.sin(altitude) * np.sin(slope)  + np.cos(altitude) * np.cos(slope)   * np.cos(azimuth - aspect)
    return 255*(shade + 1)/2

計算太陽方向及入射角皆為30的暈渲圖

plot.show(hillshade(height, 30, 30))

https://ithelp.ithome.com.tw/upload/images/20181112/20107816aSloCrDU9O.png

換個角度可以看到,暈渲圖的效果不同:

plot.show(hillshade(height, 280, 25))

https://ithelp.ithome.com.tw/upload/images/20181112/20107816WmmSZ89qxH.png

關於暈渲圖計算可以參考以下的詳細說明<https://hk.saowen.com/a/56273e57ea6fbfc59bffd030dd8a7e5741cf1d72ec4a2436dd027012ea8ff8ec>

坡度

坡度及坡向也是DEM常用的計算
坡度坡向的計算,可以使用numpy的計算
計算方式請參考
How Slope works—Help | ArcGIS for Desktop

而下面則是直接使用GDAL的模組計算

import gdal
import numpy as np
import rasterio

# 坡度
gdal.DEMProcessing('output/slope.tif', 'data/95191010dem.grd', 'slope')
with rasterio.open('output/slope.tif') as dataset:
    slope=dataset.read(1)

# 坡向
gdal.DEMProcessing('output/aspect.tif', 'data/95191010dem.grd', 'aspect')
with rasterio.open('output/aspect.tif') as dataset:
    aspect=dataset.read(1)
slope

https://ithelp.ithome.com.tw/upload/images/20181112/20107816YbmMxTl9nT.png

aspect

https://ithelp.ithome.com.tw/upload/images/20181112/201078162ufZTyAfy3.png

參考資料

http://pangea.stanford.edu/~samuelj/musings/dems-in-python-pt-1-intro.html


上一篇
Day27 計算NDVI
下一篇
Day29 Pysal空間自相關
系列文
30天精通GIS資料分析-使用Python30

尚未有邦友留言

立即登入留言