利用python和DEM數(shù)據(jù)實現(xiàn)坡度計算

利用python和DEM數(shù)據(jù)實現(xiàn)坡度計算今天我們來了解下如何利用python實現(xiàn)DEM數(shù)據(jù)的坡度計算。坡度計算一般采用擬合曲面法。擬合曲面一般采用二次曲面，即3×3的窗口，每一個窗口的中心為一個高程點。圖中列舉了四種坡度的計算方法，其中算法1的精度最高，計算效率也最高，其次是算法2。ERDASImagine中采用的是算法4，ArcMap采用的是算法2。本次也采用算法2進行實現(xiàn)。下面就開始利用python實現(xiàn)DEM數(shù)據(jù)的坡度計算。Part1導(dǎo)入相應(yīng)的庫import

osimport

numpy

as

npfrom

osgeo

import

gdalPart2坡度計算要計算的DEM數(shù)據(jù)如下所示。先讀取DEM數(shù)據(jù)，并初始化輸出。dem_ds

=

gdal.Open(dem)cell_width

=

dem_ds.GetGeoTransform()[1]cell_height

=

dem_ds.GetGeoTransform()[5]band

=

dem_ds.GetRasterBand(1)in_data

=

band.ReadAsArray().astype(np.float)out_data

=

np.ones((band.YSize,

band.XSize))

*

-99

#初始化輸出利用make_slices函數(shù)(詳見附錄)進行數(shù)據(jù)切片，make_slices返回的是一個長度為9的列表，每一個列表中的值分別對應(yīng)3*3窗口中9個輸入像素的對應(yīng)集合。然后計算坡度，注意此時計算出來的結(jié)果不包括原DEM數(shù)據(jù)中最外圍的像素點，所以輸出的時候out_data最外圍的一圈仍是初始化的值。slices

=

make_slices(in_data,(3,3))rise=

((slices[6]

+

(2

*slices[7])

+

slices[8])-

(slices[0]+(2*slices[1])+slices[2]))/(8*cell_height)run

=

((slices[2]

+

(2*slices[5])+slices[8])-(slices[0]+(2*slices[3])+slices[6]))/

(8*cell_width)dist

=

np.sqrt(np.square(rise)+np.square(run))out_data[1:-1,1:-1]

=

np.arctan(dist)*180

/np.pi利用make_raster函數(shù)(詳見附錄)輸出坡度結(jié)果，保存為tif。輸出結(jié)果如圖所示。save_raster(dem_ds,out_tif,out_data,gdal.GDT_Float32,-99)del

dem_dsDEM數(shù)據(jù)坡度圖Part3總結(jié)本文實現(xiàn)了DEM數(shù)據(jù)的坡度計算，感興趣的同學(xué)可以實現(xiàn)下另外三種坡度計算算法。如需獲取練習(xí)數(shù)據(jù)，可以關(guān)注本公眾號后，后臺發(fā)送練習(xí)數(shù)據(jù)01，即可獲取下載鏈接。Part4附錄def

make_slices(data,win_size):#計算切片大小rows

=

data.shape[0]

-

win_size[0]

+1cols

=

data.shape[1]

-

win_size[1]+1slices

=

[]for

i

in

range(win_size[0]):for

j

in

range(win_size[1]):slices.append(data[i:rows+i,j:cols+j])return

slicesdef

save_raster(input_data,

out_file,output_data,data_type,nodata=None):#

創(chuàng)建文件driver

=

gdal.GetDriverByName('GTiff')

#

數(shù)據(jù)類型必須有，因為要計算需要多大內(nèi)存空間dataset

=

driver.Create(out_file,

input_data.RasterXSize,

input_data.RasterYSize,

1,

data_type)im_proj

=

input_data.GetProjection()im_geotrans

=

input_data.GetGeoTransform()dataset.SetGeoTransform(im_geotrans)

#

寫入仿射變換參數(shù)dataset.SetPro