Slope, Aspect and Hillshade... No SA?
Like is says.
Additions: 2017-12-20
Added the hillshade calculation that was missing from the original.
Modified: 2016-12-05
Simplified aspect calculation and added the ability of how to classify 'flat' during aspect calculations. this should clean up those pesky random seemingly random aspect values in areas of small slope.
This is simply an implementation of the 3rd order finite difference method of determining slope and aspect (Horton 1981, Burrough et al. various) and used by the Spatial Analyst extension. It is also the base for a variety of terrain derivatives. I have improved the base algorithm to facilitate using larger arrays/rasters in a shorter period of time. This implementation doesn't account for nodata values... I will leave that for the Code Sharing site. I will leae the code documentation and discussion attached, but bear in mind it is much slower than the code below... but the fundamentals are the same.
This approach provides an estimate of slope which attempts to account for the cell values (ie elevations) in the 8 neighbours surrounding a core cell. A measure of the slope in both the X and Y directions is made and averaged.
Alternate measures of slope can be simply implemented using variants of the fundamental algorithms documented in the attachment. The author has also worked with implementations of simple hillshading and related neighbour functions which employ the same or similar methods of parsing 2D arrays into a 4D array.
The block or tiling process that is used to produce the subarrays has been documented in previous threads and here.
- 
| Sample surface... |
|---|
Slopes form a pyramid. Edges trimmed... |
| Slope calculations... | Aspect calculations |
|---|---|
| |
Now for some functions
The code below provides a check on the arrays then a few helper functions before the slope_a and aspect_a functions. A demo function is used to show the results for a 3xN window with know properties. I have used this easily on arrays in the 1000's in both x and y.
import sys
import numpy as np
from numpy.lib.stride_tricks import as_strided
from textwrap import dedent, indent
import matplotlib.pyplot as plt
ft = {'bool': lambda x: repr(x.astype('int32')),
'float': '{: 0.1f}'.format}
np.set_printoptions(edgeitems=10, linewidth=100, precision=2,
suppress=True, threshold=200,
formatter=ft)
np.ma.masked_print_option.set_display('-')
def _check(a, r_c, subok=False):
"""Performs the array checks necessary for stride and block.
: see arr_tools
"""
if isinstance(r_c, (int, float)):
r_c = (1, int(r_c))
r, c = r_c
if a.ndim == 1:
a = np.atleast_2d(a)
r, c = r_c = (min(r, a.shape[0]), min(c, a.shape[1]))
a = np.array(a, copy=False, subok=subok)
return a, r, c, tuple(r_c)
def stride(a, r_c=(3, 3)):
"""Provide a 2D sliding/moving view of an array.
: see arr_tools
"""
a, r, c, r_c = _check(a, r_c)
shape = (a.shape[0] - r + 1, a.shape[1] - c + 1) + r_c
strides = a.strides * 2
a_s = (as_strided(a, shape=shape, strides=strides)) #.squeeze()
return a_s
def filter_a(a_s, filter=None, cell_size=1):
"""Used by aspect, slope and hillshade to filter a raster/array
:Requires:
:--------
: a_s - a r*c*3x3 strided array
: filter - a 3x3 filter to apply to a_s
: cell_size - for slope (actual size*8, aspect (8 is required)
"""
f_dxyz = filter
cell_size= cell_size*8.0
if filter is None:
f_dxyz = np.array([[1, 2, 1], [2, 0, 2], [1, 2, 1]], dtype="float64")
a_s = a_s * f_dxyz
dz_dx = (a_s[...,:,2] - a_s[...,:,0]).sum(axis=-1)/cell_size
dz_dy = (a_s[...,2,:] - a_s[...,0,:]).sum(axis=-1)/cell_size
return dz_dx, dz_dy
def slope_a(a, cell_size=1, degrees=True, verbose=False, keepdims=False):
"""Return slope in degrees for an input array using 3rd order
: finite difference method for a 3x3 moing window view into the array.
:
:Requires:
:--------
: a - an input 2d array. X and Y represent coordinates of the Z values
: cell_size - must be in the same units as X and Y
: degrees - True, returns degrees otherwise radians
: verbose - True, to print results
: keepdims - False, to remove/squeeze extra dimensions
: filter - np.array([[1, 2, 1], [2, 0, 2], [1, 2, 1]]) **current default
:
:Notes: dzdx: sum(col2 - col0)/8*cellsize
:----- dzdy: sum(row2 - row0)/8*celsize
: Assert the array is ndim=4 even if (1,z,y,x)
: general dzdx + dzdy = dxyz
: [[a, b, c], [[1, 0, 1], [[1, 2, 1] [[1, 2, 1]
: [d, e, f] [2, 0, 2], + [0, 0, 0] = [2, 0, 2],
: [g, h, i] [1, 0, 1]] [1, 2, 1]] [1, 2, 1]]
:
"""
frmt = """\n :----------------------------------------:
:{}
:input array...\n {}
:slope values...\n {!r:}
:----------------------------------------:
"""
# ---- stride the data and calculate slope for 3x3 sliding windows
a_s = stride(a, r_c=(3,3))
if a_s.ndim < 4:
new_shape = (1,) * (4-len(a_s.shape)) + a_s.shape
a_s = a_s.reshape(new_shape)
r = a_s.shape[0]
c = a_s.shape[1]
f_dxyz = np.array([[1, 2, 1], [2, 0, 2], [1, 2, 1]], dtype="float64")
# ---- default filter, apply the filter to the array ----
#
dz_dx, dz_dy = filter_a(a_s, filter=f_dxyz, cell_size=cell_size)
#
s = np.sqrt(dz_dx**2 + dz_dy**2)
if degrees:
s = np.rad2deg(np.arctan(s))
if not keepdims:
s = np.squeeze(s)
if verbose:
from textwrap import indent, dedent # if not imported
p = " "
args = ["Results for slope_a... ",
indent(str(a), p), indent(str(s), p)]
print(dedent(frmt).format(*args))
return s
def aspect_a(a, cell_size=1, flat=0.1):
"""Return the aspect of a slope in degrees from North.
:Requires:
:--------
: a - an input 2d array. X and Y represent coordinates of the Z values
: flat - degree value, e.g. flat surface <= 0.05 deg
: 0.05 deg => 8.7e-04 0.10 deg => 1.7e-02
:
"""
if not isinstance(flat, (int, float)):
flat = 0.1
a_s = stride(a, r_c=(3,3))
if a_s.ndim < 4:
new_shape = (1,) * (4-len(a_s.shape)) + a_s.shape
a_s = a_s.reshape(new_shape)
r = a_s.shape[0]
c = a_s.shape[1]
f_dxyz = np.array([[1, 2, 1], [2, 0, 2], [1, 2, 1]], dtype="float64")
a_s = a_s * f_dxyz
# filter the array, using default filter f_dxyz above
#
dz_dx, dz_dy = filter_a(a_s, filter=f_dxyz, cell_size=1)
#
asp = np.arctan2(dz_dy, -dz_dx) # relative to East
s = np.sqrt((dz_dx*cell_size)**2 + (dz_dy*cell_size)**2) # get the slope
asp = np.rad2deg(asp)
asp = np.mod((450.0 - asp), 360.) # simplest way to get azimuth
asp = np.where(s <= flat, -1, asp)
return asp
def _demo():
"""Demonstration of calculations
:
"""
frmt = """
:------------------------------------------------------------------
:{}
:input array
{}
:slope values
{}
:zero?
{}
:aspect values
{}
:
:------------------------------------------------------------------
"""
p = " "
t = 1.0e-12
d = [[0, 0, 0, 0, t, t, t, 0, 0, 0, 0, 0, 0, 4, 4, 4, 2, 2, 2],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, t, t, t, 4, 4, 4, 2, 2, 2],
[0, 0, 0, 0, 0, 0, 0, t, t, t, 0, 0, 0, 4, 4, 4, 2, 2, 2]]
a = np.array(d, dtype='float64')
sl = slope_a(a, cell_size=5, degrees=True, verbose=False, keepdims=False)
sboo = np.where(sl==0.0, " T", " F")
asp = aspect_a(a)
z = np.ma.asarray(a)
z.mask = (z==t)
args = ["Surface properties",
z,
indent(str(sl), p),
indent(str(sboo), p),
indent(str(asp), p)]
for i in d: print(i)
print(dedent(frmt).format(*args))
return a, sl, asp, d
# ---------------------------------------------------------------------
if __name__ == "__main__":
"""Main section... """
#print("Script... {}".format(script))
a, sl, asp, d = _demo()
That's all for now.
Once slope and aspect are obtained, the hillshade can be derived. The code is fairly self-explanatory
def hillshade_a(a, cell_size=1, sun_azim=315, sun_elev=45, degrees=True):
"""Hillshade calculation as outlined in Burrough and implemented by
: esri in ArcMap and ArcGIS Pro. All measures in radians.
: z, az => sun's zenith angle an azimuth
: sl, asp => surface properties
: 255.0 * ((cos(z) * cos(sl)) + (sin(z) * sin(sl) * cos(az-asp)))
:
"""
s_azi = np.deg2rad(sun_azim)
s_elev = np.deg2rad(90.0 - sun_elev)
a_a = aspect_a(a, degrees=False)
a_s = slope_a(a, cell_size=cell_size, degrees=False)
out = 255*((np.cos(s_elev) * np.cos(a_s)) +
(np.sin(s_elev) * np.sin(a_s) * np.cos(s_azi - a_a)))
out = np.where(out < 0, 0, out)
return out.astype('int')Other terrain derivates will be added soon.