Elevation (XYZ) data to slope/gradient map using python

Question

I have a text file with Easting (x), Northing (y), and Elevation data (z) as shown below:

   x            y         z
241736.69   3841916.11  132.05
241736.69   3841877.89  138.76
241736.69   3841839.67  142.89
241736.69   3841801.45  148.24
241736.69   3841763.23  157.92
241736.69   3841725.02  165.01
241736.69   3841686.80  171.86
241736.69   3841648.58  178.80
241736.69   3841610.36  185.26
241736.69   3841572.14  189.06
241736.69   3841533.92  191.28
241736.69   3841495.71  193.27
241736.69   3841457.49  193.15
241736.69   3841419.27  194.85
241736.69   3841381.05  192.31
241736.69   3841342.83  188.73
241736.69   3841304.61  183.68
241736.69   3841266.39  176.97
241736.69   3841228.18  160.83
241736.69   3841189.96  145.69
241736.69   3841151.74  129.09
241736.69   3841113.52  120.03
241736.69   3841075.30  111.84
241736.69   3841037.08  104.82
241736.69   3840998.86  101.63
241736.69   3840960.65  97.66
241736.69   3840922.43  93.38
241736.69   3840884.21  88.84
...

I can get an elevation map from the above data easily with plt.contour and plt.contourf as shown below:

However,I am trying to get a slope map of the data I have, something like this:

What I tried to do is to convert my XYZ data to DEM using GDAL as explained here, and loading the DEM with richdem, as explained here, but I am getting wrong slope values.

The results I get from converting to .tif:

This is the code I have tried with richdem:

import richdem as rd

dem_path = 'convertedXYZ.tif'
dem = rd.LoadGDAL(dem_path, no_data=-9999)
slope = rd.TerrainAttribute(dem, attrib='slope_riserun')
rd.rdShow(slope, axes=True, cmap='gist_yarg', figsize=(16, 9))

the plot I am getting:

The values on the colorbar are too high to be correct and the plot must be inverted to match the above plots (not my main issue right now).

I am not an expert when using python for GIS purposes (I mainly use python for data analysis) and I am hoping this is not as complicated as I think it is.

Answer 1

eI was able to write a function that does the job correctly but first I need to give credit to this answer for saving me some time with writing my own moving windows function (works perfectly!):

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from tqdm import trange


def window3x3(arr, shape=(3, 3)):
    r_win = np.floor(shape[0] / 2).astype(int)
    c_win = np.floor(shape[1] / 2).astype(int)
    x, y = arr.shape
    for i in range(x):
        xmin = max(0, i - r_win)
        xmax = min(x, i + r_win + 1)
        for j in range(y):
            ymin = max(0, j - c_win)
            ymax = min(y, j + c_win + 1)
            yield arr[xmin:xmax, ymin:ymax]


def gradient(XYZ_file, min=0, max=15, figsize=(15, 10), **kwargs):

    """

    :param XYZ_file: XYZ file in the following format: x,y,z (inlcuding headers)
    :param min: color bar minimum range.
    :param max: color bar maximum range.
    :param figsize: figure size.
    :param kwargs:
           plot: to plot a gradient map. Default is True.
    :return: returns an array with the shape of the grid with the computed slopes


    The algorithm calculates the gradient using a first-order forward or backward difference on the corner points, first
    order central differences at the boarder points, and a 3x3 moving window for every cell with 8 surrounding cells (in
    the middle of the grid) using a third-order finite difference weighted by reciprocal of squared distance

    Assumed 3x3 window:

                        -------------------------
                        |   a   |   b   |   c   |
                        -------------------------
                        |   d   |   e   |   f   |
                        -------------------------
                        |   g   |   h   |   i   |
                        -------------------------


    """

    kwargs.setdefault('plot', True)

    grid = XYZ_file.to_numpy()

    nx = XYZ_file['x'].unique().size
    ny = XYZ_file['y'].unique().size

    xs = grid[:, 0].reshape(ny, nx, order='F')
    ys = grid[:, 1].reshape(ny, nx, order='F')
    zs = grid[:, 2].reshape(ny, nx, order='F')
    dx = abs((xs[:, 1:] - xs[:, :-1]).mean())
    dy = abs((ys[1:, :] - ys[:-1, :]).mean())

    gen = window3x3(zs)
    windows_3x3 = np.asarray(list(gen))
    windows_3x3 = windows_3x3.reshape(ny, nx)

    dzdx = np.empty((ny, nx))
    dzdy = np.empty((ny, nx))
    loc_string = np.empty((ny, nx), dtype="S25")

    for ax_y in trange(ny):
        for ax_x in range(nx):

            # corner points
            if ax_x == 0 and ax_y == 0:  # top left corner
                dzdx[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][0][1] - windows_3x3[ax_y, ax_x][0][0]) / dx
                dzdy[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][1][0] - windows_3x3[ax_y, ax_x][0][0]) / dy
                loc_string[ax_y, ax_x] = 'top left corner'

            elif ax_x == nx - 1 and ax_y == 0:  # top right corner
                dzdx[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][0][1] - windows_3x3[ax_y, ax_x][0][0]) / dx
                dzdy[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][1][1] - windows_3x3[ax_y, ax_x][0][1]) / dy
                loc_string[ax_y, ax_x] = 'top right corner'

            elif ax_x == 0 and ax_y == ny - 1:  # bottom left corner
                dzdx[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][1][1] - windows_3x3[ax_y, ax_x][1][0]) / dx
                dzdy[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][1][0] - windows_3x3[ax_y, ax_x][0][0]) / dy
                loc_string[ax_y, ax_x] = 'bottom left corner'

            elif ax_x == nx - 1 and ax_y == ny - 1:  # bottom right corner
                dzdx[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][1][1] - windows_3x3[ax_y, ax_x][1][0]) / dx
                dzdy[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][1][1] - windows_3x3[ax_y, ax_x][0][1]) / dy
                loc_string[ax_y, ax_x] = 'bottom right corner'

            # top boarder
            elif (ax_y == 0) and (ax_x != 0 and ax_x != nx - 1):
                dzdx[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][0][-1] - windows_3x3[ax_y, ax_x][0][0]) / (2 * dx)
                dzdy[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][1][1] - windows_3x3[ax_y, ax_x][0][1]) / dy
                loc_string[ax_y, ax_x] = 'top boarder'

            # bottom boarder
            elif ax_y == ny - 1 and (ax_x != 0 and ax_x != nx - 1):
                dzdx[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][1][-1] - windows_3x3[ax_y, ax_x][1][0]) / (2 * dx)
                dzdy[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][1][1] - windows_3x3[ax_y, ax_x][0][1]) / dy
                loc_string[ax_y, ax_x] = 'bottom boarder'

            # left boarder
            elif ax_x == 0 and (ax_y != 0 and ax_y != ny - 1):
                dzdx[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][1][1] - windows_3x3[ax_y, ax_x][1][0]) / dx
                dzdy[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][-1][0] - windows_3x3[ax_y, ax_x][0][0]) / (2 * dy)
                loc_string[ax_y, ax_x] = 'left boarder'

            # right boarder
            elif ax_x == nx - 1 and (ax_y != 0 and ax_y != ny - 1):
                dzdx[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][1][1] - windows_3x3[ax_y, ax_x][1][0]) / dx
                dzdy[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][-1][-1] - windows_3x3[ax_y, ax_x][0][-1]) / (2 * dy)
                loc_string[ax_y, ax_x] = 'right boarder'

            # middle grid
            else:
                a = windows_3x3[ax_y, ax_x][0][0]
                b = windows_3x3[ax_y, ax_x][0][1]
                c = windows_3x3[ax_y, ax_x][0][-1]
                d = windows_3x3[ax_y, ax_x][1][0]
                f = windows_3x3[ax_y, ax_x][1][-1]
                g = windows_3x3[ax_y, ax_x][-1][0]
                h = windows_3x3[ax_y, ax_x][-1][1]
                i = windows_3x3[ax_y, ax_x][-1][-1]

                dzdx[ax_y, ax_x] = ((c + 2 * f + i) - (a + 2 * d + g)) / (8 * dx)
                dzdy[ax_y, ax_x] = ((g + 2 * h + i) - (a + 2 * b + c)) / (8 * dy)
                loc_string[ax_y, ax_x] = 'middle grid'

    hpot = np.hypot(abs(dzdy), abs(dzdx))
    slopes_angle = np.degrees(np.arctan(hpot))
    if kwargs['plot']:
        slopes_angle[(slopes_angle < min) | (slopes_angle > max)]

        plt.figure(figsize=figsize)
        plt.pcolormesh(xs, ys, slopes_angle, cmap=plt.cm.gist_yarg, vmax=max, vmin=min)

        plt.colorbar()
        plt.tight_layout()
        plt.show()

    return slopes_angle


if __name__ == '__main__':

    XYZ = pd.read_csv('xyz_file')
    slopes = gradient(XYZ)

and the final plot:

Answer 2

Assuming your data is in an n x 3 Numpy array, first reinterpret just the elevation column as a matrix (representing the uniform grid):

m=data[:,2].reshape(ny,nx)

Then perform several slices and subtractions to get the derivatives at the cell centers:

dx=m[:,1:]-m[:,:-1]
dy=m[1:,:]-m[:-1,:]
mag=numpy.hypot(dx[1:,:]+dx[:-1,:],
                dy[:,1:]+dy[:,:-1])
mag*=abs(data[1][1]-data[1][0])/2

The coefficient corrects the units (which would otherwise be meters per cell, not per meter) and converts the sums to averages. (If the spacing in each dimension differed, you would scale the arguments to hypot separately.) Note that the resulting array is one smaller in each dimension than the input; other, more complicated differencing schemes can be used if the sizes need to be the same. numpy.gradient implements some of them, allowing a simple

mag=numpy.hypot(*numpy.gradient(m,abs(data[1][1]-data[1][0])))