# Slope, Aspect and Hillshade... No SA?

Blog Post created by Dan_Patterson on Mar 12, 2016

Like is says.

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...

``>>> Dem...[[0 1 2 3 2 1 0] [1 2 3 4 3 2 1] [2 3 4 5 4 3 2] [3 4 5 5 5 4 3] [2 3 4 5 4 3 2] [1 2 3 4 3 2 1] [0 1 2 3 2 1 0]]Slope...[[ 15.79  15.79  11.31  15.79  15.79] [ 15.79  13.9    8.53  13.9   15.79] [ 11.31   8.53   0.     8.53  11.31] [ 15.79  13.9    8.53  13.9   15.79] [ 15.79  15.79  11.31  15.79  15.79]]Aspect...[[ 315.  315.    0.   45.   45.] [ 315.  315.    0.   45.   45.] [ 270.  270.  270.   90.   90.] [ 225.  225.  180.  135.  135.] [ 225.  225.  180.  135.  135.]]``

Slope calculations... Aspect calculations

``>>> # ---- Slope calculation using>>> # finite difference ---->>>>>>Slope face...[[1 1 1][3 3 3][5 5 5]]N     dx    deg.    %(1 )   0.5   71.6  400.0(2 )   1.0   56.3  200.0(3 )   2.0   36.9  100.0(4 )   4.0   20.6   50.0(5 )   6.0   14.0   33.3(6 )   8.0   10.6   25.0(7 )  10.0    8.5   20.0(8 )  12.5    6.8   16.0(9 )  15.0    5.7   13.3(10)  20.0    4.3   10.0(11)  25.0    3.4    8.0(12)  40.0    2.1    5.0(13)  80.0    1.1    2.5(14) 100.0    0.9    2.0>>> # ---- dx is the cell width ---->>> # It is a north facing slope,>>> # with dx=4 and dz=2, you get a>>> # 20.6 deg. or 50% slope>>>``

``>>> # ---- Some windows with>>> #   aspect determination ---->>> (0) Array aspect...0.0[[1 0 1] [2 2 2] [3 4 3]](1) Array aspect...315.0[[0 1 2] [1 2 3] [2 3 4]](2) Array aspect...270.0... snip ...(4) Array aspect...180.0[[3 4 3] [2 2 2] [1 0 1]](5) Array aspect...135.0[[4 3 2] [3 2 1] [2 1 0]](6) Array aspect...90.0[[3 2 1] [4 2 0] [3 2 1]](7) Array aspect...45.0[[2 1 0] etc... etc ... ``

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 sysimport numpy as npfrom numpy.lib.stride_tricks import as_stridedfrom textwrap import dedent, indentimport matplotlib.pyplot as pltft = {'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), min(c, a.shape))    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 - r + 1, a.shape - c + 1) + r_c    strides = a.strides * 2    a_s = (as_strided(a, shape=shape, strides=strides)) #.squeeze()    return a_sdef 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_dydef 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    c = a_s.shape    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    c = a_s.shape    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 aspdef _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.