xrspatial.curvature.curvature

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

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

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")
../../_images/xrspatial-curvature-curvature-1_00.png

(png, hires.png, pdf)

../../_images/xrspatial-curvature-curvature-1_01.png

(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