import os
import subprocess
from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt
import math

# syntax: command (cmd) = gdaltool mode input file (+path if necessary) output file
cmd = "gdaldem slope data/dem.tif slope2.tif -compute_edges"
#os.system(cmd)
subprocess.call(cmd.split())

# slope
slp1 = gdal.Open("slope2.tif")
slp1Array = slp1.GetRasterBand(1).ReadAsArray()
plt.figure()
plt.imshow(slp1Array)
plt.colorbar()
plt.show()

# aspect
dem = gdal.Open('data/dem.tif')
aspect = gdal.DEMProcessing("aspect.tif", dem, "aspect", computeEdges=True)
aspectArray = aspect.GetRasterBand(1).ReadAsArray()
plt.figure()
plt.imshow(aspectArray, cmap = "magma")
plt.colorbar()
plt.show()

# slope3
dem = gdal.Open('data/dem.tif')
slp3 = gdal.DEMProcessing("slope3.tif", dem, "slope", computeEdges=True)
slp3Array = slp3.GetRasterBand(1).ReadAsArray()
plt.figure()
plt.imshow(slp3Array)
plt.colorbar()
plt.show()


# hillshade
dem = gdal.Open('data/dem.tif')
hillshade = gdal.DEMProcessing("hillshade.tif", dem, "hillshade", computeEdges=True)
hsArray = hillshade.GetRasterBand(1).ReadAsArray()
plt.figure()
plt.imshow(hsArray, cmap = "gray")
plt.colorbar()
plt.show()


# calculate difference
diff = slp3Array - aspectArray
plt.figure()
plt.imshow(diff)
plt.colorbar()
plt.show()

# mean slope angle
print(np.mean(slp3Array))

### all with richdem ###

## somehow does not work yet ##

# import richdem as rd


# dem = rd.LoadGDAL("data/dem.tif")
# slp4 = rd.TerrainAttribute(dem, attrib="slope_degrees")
# plt.figure()
# plt.imshow(slp4)
# plt.colorbar()
# plt.show()
# rd.SaveGDAL("slope4.tif", slp4)


slp1 = slp3 = None







import os
import subprocess
from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt
import math

# syntax: command (cmd) = gdaltool mode input file (+path if necessary) output file
cmd = "gdaldem slope data/dem.tif slope2.tif -compute_edges"
#os.system(cmd)
subprocess.call(cmd.split())

# slope
slp1 = gdal.Open("slope2.tif")
slp1Array = slp1.GetRasterBand(1).ReadAsArray()
plt.figure()
plt.imshow(slp1Array)
plt.colorbar()
plt.show()

# aspect
dem = gdal.Open('data/dem.tif')
aspect = gdal.DEMProcessing("aspect.tif", dem, "aspect", computeEdges=True)
aspectArray = aspect.GetRasterBand(1).ReadAsArray()
plt.figure()
plt.imshow(aspectArray, cmap = "magma")
plt.colorbar()
plt.show()

# slope3
dem = gdal.Open('data/dem.tif')
slp3 = gdal.DEMProcessing("slope3.tif", dem, "slope", computeEdges=True)
slp3Array = slp3.GetRasterBand(1).ReadAsArray()
plt.figure()
plt.imshow(slp3Array)
plt.colorbar()
plt.show()


# hillshade
dem = gdal.Open('data/dem.tif')
hillshade = gdal.DEMProcessing("hillshade.tif", dem, "hillshade", computeEdges=True)
hsArray = hillshade.GetRasterBand(1).ReadAsArray()
plt.figure()
plt.imshow(hsArray, cmap = "gray")
plt.colorbar()
plt.show()


# calculate difference
diff = slp3Array - aspectArray
plt.figure()
plt.imshow(diff)
plt.colorbar()
plt.show()

# mean slope angle
print(np.mean(slp3Array))

### all with richdem ###

## somehow does not work yet ##

# import richdem as rd


# dem = rd.LoadGDAL("data/dem.tif")
# slp4 = rd.TerrainAttribute(dem, attrib="slope_degrees")
# plt.figure()
# plt.imshow(slp4)
# plt.colorbar()
# plt.show()
# rd.SaveGDAL("slope4.tif", slp4)


slp1 = slp3 = None







# syntax: command (cmd) = gdaltool mode input file (+path if necessary) output file #os.system(cmd) # slope # aspect # slope3 # hillshade # calculate difference # mean slope angle ### all with richdem ### ## somehow does not work yet ## # import richdem as rd # dem = rd.LoadGDAL("data/dem.tif") # slp4 = rd.TerrainAttribute(dem, attrib="slope_degrees") # plt.figure() # plt.imshow(slp4) # plt.colorbar() # plt.show() # rd.SaveGDAL("slope4.tif", slp4)