我還有一個問題。在DEM上計算坡度
我必須計算數字高程模型中每個單元的斜率。應該爲形狀爲3×3的移動窗口中的每個中心單元計算斜率,並且根據公式:
斜率= max | x9-xi |/A;其中值i是從1到8的值,並且值x9是窗口的中心。 A是到相鄰小區中點的距離。因此,對於與中心對角的單元,距離(A)是sqrt(2)乘以分辨率,對於其他單元則等於分辨率。
所以我不知道的是如何編寫一個代碼,將不同於對角線,而不是?我創建了一個空的numpy數組,其中'沒有值',我想要具有相同分辨率的斜率值。我知道我必須遍歷行和列,並且必須忽略第一行和最後一行和列。 到目前爲止我的代碼:
import numpy as np
import matplotlib.pyplot as plt
dem=np.loadtxt('dem.txt',dtype='int',delimiter=',')
(rows,cols)=np.shape(dem)
slope2=np.zeros((rows,cols))
for r in xrange(1,rows-1):
for c in xrange(1,cols-1):
temp_win=dem[r-1:r+2,c-1:c+2]
mid= temp_win[1,1]
max_d=np.max([temp_win[0,0]],[temp_win[0,2]],[temp_win[2,0]],[temp_win[2,2]])
max_1=np.max([temp_win[1,1]],[temp_win[1,0]],[temp_win[1,2]],[temp_win[2,1]])
slope = (np.max([np.abs(mid-max_d)/np.sqrt(2)]),np.max([np.abs(mid-s1/np.sqrt(2))])
slope_2 = slope[r,c]
沒有人有任何想法,我會很感激一些幫助?
使用GIS包而不是重新發明輪子http://desktop.arcgis.com/en/desktop/latest/tools/spatial-analyst-toolbox/how-slope-works.htm但如果您需要模擬其中的代碼有變體或諮詢Burrough和McDonnel 1998,地理信息系統原理,牛津大學。按...或等效文本 – 2015-10-05 20:58:41
您通常使用'numpy.gradient'來代替迭代每個單元格。它將產生兩個具有y和x梯度的2D陣列,然後可以使用它們來計算斜率/方面。然而它使用了一個非常不同的算法。 –
我現在編輯我的代碼,但它說最後2行中的語法無效,所以我仍然需要一些幫助 – anjisa