4

gradientm is a matlab function that Calculate gradient, slope, and aspect of data grid

Syntax

[ASPECT, SLOPE, gradN, gradE] = gradientm(Z, R)

Description

[ASPECT, SLOPE, gradN, gradE] = gradientm(Z, R) computes the slope, aspect, and north and east components of the gradient for a regular data grid Z with respect to reference R. If the grid contains elevations in meters, the resulting aspect and slope are in units of degrees clockwise from north and up from the horizontal. The north and east gradient components are the change in the map variable per meter of distance in the north and east directions. The computation uses finite differences for the map variable on the default earth ellipsoid.

Problem

I want something equivalent in python. I have founded this python library PyDEM but it is for python 2.x but I am working with python 3.x

this is The way of Use the DEMProcessor to compute slope and aspect for the DEM files by PyDEM:

# needs to match above command
filename = 'Shasta-30m-DEM.tif'

# instantiate a processor object
processor = DEMProcessor(filename)

# get magnitude of slope and aspect
mag, aspect = processor.calc_slopes_directions()
pd shah
  • 1,346
  • 2
  • 14
  • 26

1 Answers1

15

step1: install GDAL

step2: install rasterio

step3: easy as py:

from osgeo import gdal
import numpy as np
import rasterio


def calculate_slope(DEM):
    gdal.DEMProcessing('slope.tif', DEM, 'slope')
    with rasterio.open('slope.tif') as dataset:
        slope=dataset.read(1)
    return slope

def calculate_aspect(DEM):
    gdal.DEMProcessing('aspect.tif', DEM, 'aspect')
    with rasterio.open('aspect.tif') as dataset:
        aspect=dataset.read(1)
    return aspect

slope=calculate_slope('DEM.tif')
aspect=calculate_aspect('DEM.tif')

print(type(slope))
print(slope.dtype)
print(slope.shape)
pd shah
  • 1,346
  • 2
  • 14
  • 26
  • 1
    how can I get data about which slope value related to which coordinate? Do you have any idea about that? Currently it's returning numpy.ndarray with bunch of values inside – Sabit Rakhim May 21 '21 at 09:11
  • @SabitRakhim you could use xArray instead of numpy for that. You can modify the code provided by @pd_shah adding the xArray import first: `import xarray as xr` and then modifying the returned data into: `return xr.open_rasterio('slope.tif')` – clausmichele Apr 20 '22 at 09:33