forked from uva-hydroinformatics/wetland_id
-
Notifications
You must be signed in to change notification settings - Fork 0
/
pygeonet_slope_curvature.py
119 lines (108 loc) · 5.2 KB
/
pygeonet_slope_curvature.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
import numpy as np
from scipy import stats
from time import clock
import statsmodels.api as sm
import matplotlib.pyplot as plt
from pygeonet_rasterio import *
from pygeonet_plot import *
def compute_dem_slope(filteredDemArray, pixelDemScale):
slopeXArray,slopeYArray = np.gradient(filteredDemArray, pixelDemScale)
slopeDemArray = np.sqrt(slopeXArray**2 + slopeYArray**2)
slopeMagnitudeDemArrayQ = slopeDemArray
slopeMagnitudeDemArrayQ = np.reshape(slopeMagnitudeDemArrayQ,
np.size(slopeMagnitudeDemArrayQ))
slopeMagnitudeDemArrayQ = slopeMagnitudeDemArrayQ[~np.isnan(slopeMagnitudeDemArrayQ)]
# Computation of statistics of slope
print ' slope statistics'
print ' angle min:', np.arctan(np.percentile(slopeMagnitudeDemArrayQ,0.1))*180/np.pi
print ' angle max:', np.arctan(np.percentile(slopeMagnitudeDemArrayQ,99.9))*180/np.pi
print ' mean slope:', np.nanmean(slopeDemArray[:])
print ' stdev slope:', np.nanstd(slopeDemArray[:])
return slopeDemArray
def compute_dem_curvature(demArray, pixelDemScale, curvatureCalcMethod):
gradXArray, gradYArray = np.gradient(demArray, pixelDemScale)
slopeArrayT = np.sqrt(gradXArray**2 + gradYArray**2)
if curvatureCalcMethod == 'geometric':
#Geometric curvature
print ' using geometric curvature'
gradXArrayT = np.divide(gradXArray,slopeArrayT)
gradYArrayT = np.divide(gradYArray,slopeArrayT)
elif curvatureCalcMethod=='laplacian':
# do nothing..
print ' using laplacian curvature'
gradXArrayT = gradXArray
gradYArrayT = gradYArray
gradGradXArray,tmpy = np.gradient(gradXArrayT,pixelDemScale)
tmpx,gradGradYArray = np.gradient(gradYArrayT,pixelDemScale)
curvatureDemArray = gradGradXArray + gradGradYArray
curvatureDemArray[np.isnan(curvatureDemArray)] = 0
del tmpy, tmpx
# Computation of statistics of curvature
print ' curvature statistics'
tt = curvatureDemArray[~np.isnan(curvatureDemArray[:])]
print ' non-nan curvature cell number:', tt.shape[0]
finiteCurvatureDemList = curvatureDemArray[np.isfinite(curvatureDemArray[:])]
print ' non-nan finite curvature cell number:', finiteCurvatureDemList.shape[0]
curvatureDemMean = np.nanmean(finiteCurvatureDemList)
curvatureDemStdDevn = np.nanstd(finiteCurvatureDemList)
print ' mean: ', curvatureDemMean
print ' standard deviation: ', curvatureDemStdDevn
return curvatureDemArray, curvatureDemMean, curvatureDemStdDevn
def compute_quantile_quantile_curve(x):
print 'getting qqplot estimate'
if not hasattr(defaults, 'figureNumber'):
defaults.figureNumber = 0
defaults.figureNumber = defaults.figureNumber + 1
plt.figure(defaults.figureNumber)
res = stats.probplot(x, plot=plt)
res1 = sm.ProbPlot(x, stats.t, fit=True)
print res1
return res
def main():
## plt.switch_backend('agg')
filteredDemArray = read_geotif_filteredDEM()
# filteredDemArray = read_dem_from_geotiff("NF_breachdep50.tif", r"D:\2ndStudy_ONeil\Terrain_Processing\whitebox\Site3_rt29")
# Computing slope
print 'computing slope'
slopeDemArray = compute_dem_slope(filteredDemArray, Parameters.demPixelScale)
slopeDemArray[np.isnan(filteredDemArray)] = np.nan
# Writing the curvature array
outfilepath = Parameters.geonetResultsDir
# outfilepath = r"D:\Wetlands_2\Rt29\Study2_Indices\A"
demName = Parameters.demFileName.split('.')[0]
# outfilename = demName +'_slope.tif'
outfilename = "NF_breachdep50_slope.tif"
write_geotif_generic(slopeDemArray, outfilepath, outfilename)
# Computing curvature
print 'computing curvature'
## curvatureDemArrayIn = filteredDemArray
curvatureDemArray, curvatureDemMean, \
curvatureDemStdDevn = compute_dem_curvature(filteredDemArray,
Parameters.demPixelScale,
defaults.curvatureCalcMethod)
curvatureDemArray[np.isnan(filteredDemArray)] = np.nan
# Writing the curvature array
outfilename = demName +'_curvature.tif'
# outfilename = "NF_breachdep50_curvature.tif"
write_geotif_generic(curvatureDemArray, outfilepath, outfilename)
# plotting the curvature image
if defaults.doPlot == 1:
raster_plot(curvatureDemArray, 'Curvature DEM')
finiteCurvatureDemList = curvatureDemArray[np.isfinite(curvatureDemArray[:])]
#*************************************************
#Compute curvature quantile-quantile curve
# This seems to take a long time ... is commented for now
#print 'computing curvature quantile-quantile curve'
#osm,osr = compute_quantile_quantile_curve(finiteCurvatureDemList)
#print osm[0]
#print osr[0]
thresholdCurvatureQQxx = 1
# have to add method to automatically compute the thresold
# .....
# .....
#*************************************************
if __name__ == '__main__':
t0 = clock()
main()
t1 = clock()
print "time taken to complete slope and curvature computation:", t1-t0, " seconds"