DAY 28
0
AI & Data

## Day28 20米DEM資料處理

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

• 網格檔案GRD操作
• 坡度坡向

## 網格檔案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)
``````

``````import numpy as np
height
``````

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)
``````

``````def hillshade(array, azimuth, angle_altitude):

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)
``````

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

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

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

## 坡度

How Slope works—Help | ArcGIS for Desktop

``````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:

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

``````aspect
``````

## 參考資料

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