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.

Your output needs to be flipped, not rotated =). Anyway, have you checked the DEM after you convert the data (your `convertedXYZ.tif'`)? do they look ok? Also, have you read the documentation for richdem? the `no_data` option says: "no_data A value indicating which cell values should be treated as special NoData cells. (See Concepts, TODO)." Is it this what you are trying to achieve by setting it to -9999? — CAPSLOCK, Apr 24 '19 at 11:38
@Gio - Yes you are right, it should be flipped rather than rotated. I have not checked the `convertedXYZ.tif` as I do not know how to perform the checking process. What should I be looking for exactly? And regarding the `no_data` value I have selected, the `LoadGDAL` class requires to assign a value to `no_data` in order to run. As I have understood from their documentation before, it is the value that represents cells that are not part of the data (so i choose -9999 as per their examples) — Tabbakhh, Apr 24 '19 at 12:41
@Gio - I opened the `convertedXYZ.tif` and added the plot to the question (FYI). — Tabbakhh, Apr 24 '19 at 12:52
Is the input fully gridded, as it appears from the first few lines? If so, couldn’t you just apply a finite difference? — Davis Herring, Apr 24 '19 at 13:23
@DavisHerring - What do you mean by fully gridded? The input file goes on for 300k+ lines (same structure) — Tabbakhh, Apr 24 '19 at 13:32
@Tabbakhh: I mean that each *x* value exists for every *y* value and vice versa. It’s even easier if they’re uniformly spaced (as they appear to be). — Davis Herring, Apr 24 '19 at 14:06
@DavisHerring - Yes they are indeed. 38 x 38 m grid to be exact. May you lead me to how to perform the finite difference that you have suggested? — Tabbakhh, Apr 24 '19 at 14:29

Tabbakhh · Accepted Answer · 2019-07-04T07:37:37.480

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:

Davis Herring · Answer 2 · 2019-06-20T23:47:25.957

0

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])))

edited Jun 20 '19 at 23:47

answered Apr 24 '19 at 15:13

Davis Herring

36,443
4
48
76

Thank you for your answer. I think you meant to type m=data[:,2].reshape(ny,nx) rather than m=data[:,3].reshape(ny,nx) in the first line. I ran your code and what I understand is that `mag` needs to be converted to angles afterwards correct?. Also, it is important for me to have the final array with the same length as the input data. What other solution can I look into in regards to this? Are there libraries that can take care of the math for me rather than hard coding it? – Tabbakhh Apr 28 '19 at 11:10
@Tabbakhh: Yes, 2: edited. `mag` is slopes; they could be converted to pitch angles (with `atan`, as always) if desired. For having the same number of samples, you just have to use centered and, on the edges, one-sided differences; I can edit those in later, but you can just use [`numpy.gradient`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.gradient.html) with some of the corrections above for units instead of what I think you meant to be “hand coding”. – Davis Herring Apr 28 '19 at 14:20
Is it still possible to help me out with this by modifying the code. I couldn't get it to work. – Tabbakhh Jun 20 '19 at 05:57

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

2 Answers2