xrspatial.curvature.curvature

xrspatial.curvature.curvature(agg: xarray.core.dataarray.DataArray, name: Optional[str] = 'curvature') → xarray.core.dataarray.DataArray

Calculates, for all cells in the array, the curvature (second derivative) of each cell based on the elevation of its neighbors in a 3x3 grid. A positive curvature indicates the surface is upwardly convex. A negative value indicates it is upwardly concave. A value of 0 indicates a flat surface.

Units of the curvature output raster are one hundredth (1/100) of a z-unit.

Parameters

• agg (xarray.DataArray) – 2D NumPy, CuPy, NumPy-backed Dask, or Cupy-backed Dask array of elevation values. Must contain res attribute.

• name (str, default='curvature') – Name of output DataArray.

Returns

curvature_agg – 2D aggregate array of curvature values. All other input attributes are preserved.

Return type

xarray.DataArray, of the same type as agg

References

• arcgis: https://pro.arcgis.com/en/pro-app/latest/tool-reference/spatial-analyst/how-curvature-works.htm # noqa

Examples

import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

from xrspatial import generate_terrain, curvature

# Generate Example Terrain
W = 500
H = 300

template_terrain = xr.DataArray(np.zeros((H, W)))
x_range=(-20e6, 20e6)
y_range=(-20e6, 20e6)

terrain_agg = generate_terrain(
    template_terrain, x_range=x_range, y_range=y_range
)

# Edit Attributes
terrain_agg = terrain_agg.assign_attrs(
    {
        'Description': 'Example Terrain',
        'units': 'km',
        'Max Elevation': '4000',
    }
)

terrain_agg = terrain_agg.rename({'x': 'lon', 'y': 'lat'})
terrain_agg = terrain_agg.rename('Elevation')

# Create Curvature Aggregate Array
curvature_agg = curvature(agg = terrain_agg, name = 'Curvature')

# Edit Attributes
curvature_agg = curvature_agg.assign_attrs(
    {
        'Description': 'Curvature',
        'units': 'rad',
    }
)

# Plot Terrain
terrain_agg.plot(cmap = 'terrain', aspect = 2, size = 4)
plt.title("Terrain")
plt.ylabel("latitude")
plt.xlabel("longitude")

# Plot Curvature
curvature_agg.plot(aspect = 2, size = 4)
plt.title("Curvature")
plt.ylabel("latitude")
plt.xlabel("longitude")

(png, hires.png, pdf)

(png, hires.png, pdf)

>>> print(terrain_agg[200:203, 200:202])
<xarray.DataArray 'Elevation' (lat: 3, lon: 2)>
array([[1264.02296597, 1261.947921  ],
       [1285.37105519, 1282.48079719],
       [1306.02339636, 1303.4069579 ]])
Coordinates:
* lon      (lon) float64 -3.96e+06 -3.88e+06
* lat      (lat) float64 6.733e+06 6.867e+06 7e+06
Attributes:
    res:            (80000.0, 133333.3333333333)
    Description:    Example Terrain
    units:          km
    Max Elevation:  4000

>>> print(curvature_agg[200:203, 200:202])
<xarray.DataArray 'Curvature' (lat: 3, lon: 2)>
array([[-8.14280993e-08, -5.46873377e-08],
       [ 2.85253552e-08,  1.24471207e-08],
       [ 8.91172878e-08,  1.58492666e-07]])
Coordinates:
* lon      (lon) float64 -3.96e+06 -3.88e+06
* lat      (lat) float64 6.733e+06 6.867e+06 7e+06
Attributes:
    res:            (80000.0, 133333.3333333333)
    Description:    Curvature
    units:          rad
    Max Elevation:  4000

previous

xrspatial.aspect.aspect

next

xrspatial.hillshade.hillshade