r.soildepth2

#!/usr/bin/env python
#
############################################################################
#
# MODULE:       r.soildepth
# AUTHOR:       Isaac Ullah, Arizona State University
# PURPOSE:      Create soil depth map based on slope, accumulated flow, and hillslope curvature.
# ACKN:         National Science Foundation Grant #BCS0410269
#               Based on information from Heimsath et al, 1997
# COPYRIGHT:    (C) 2016 by Isaac Ullah
#            This program is free software under the GNU General Public
#            License (>=v2). Read the file COPYING that comes with GRASS
#            for details.
#
#############################################################################


#%Module
#%  description: Create soil depth map based on slope, accumulated flow, and hillslope curvature.
#%End

#%option G_OPT_R_INPUT
#% key: elevation
#% description: Input elevation map (DEM)
#% required : yes
#%END
#%option G_OPT_R_OUTPUT
#% key: bedrock
#% description: Output bedrock elevation map
#% required : yes
#%END
#%option G_OPT_R_OUTPUT
#% key: soildepth
#% description: Output soil depth map
#% required : yes
#%END


#%option
#% key: min
#% type: double
#% description: What is the actual (empirical) minimum soil depth in this area (in meters)? (NOTE: The smoothing operation will trim extreme values, so resulting soildepth map may not have minimum depths exactly equal to this number)
#% answer: 0.0001
#% required : yes
#%END
#%option
#% key: max
#% type: double
#% description: What is the actual (empirical) maximum soil depth in this area (in meters)? (NOTE: The smoothing operation will trim extreme values, so resulting soildepth map may not have maximum depths exactly equal to this number)
#% answer: 20.0
#% required : yes
#%END
#%option
#% key: slopebreaks
#% type: string
#% description: Slope remapping pairs in the form: "x1,y1;x2,y2" where x is a slope value (in degrees), and y is a number between 0 and 1 to rescale to. (Two pairs are needed, in order, lower slope then higher slope.)
#% answer: 5,0.25;45,0.05
#% required : yes
#%END
#%option
#% key: flaccbreaks
#% type: string
#% description: Flacc remapping pairs in the form: "x1,y1;x2,y2" where x is a Flow Accumulation value, and y is a number between 0 and 1 to rescale to. (Two pairs are needed, in order, lower Flacc then higher Flacc.)
#% answer: 50,0.05;250,0.25
#% required : yes
#%END
#%option
#% key: curvebreaks
#% type: string
#% description: Breakpoints in mean curvature values for maximum and minimum soil depths in the form of: "x1,y1;x2,y2" where x1 and x2 are the concavity and convexity value breakpoints respectively, y1 is a number between 0 and -1 and y2 is a number between 0 and 1. (Note: if x1 and/or x2 are too large, the real max and/or min curvature will be used.)
#% answer: -0.0025,-0.5;0.0025,0.5
#% required : yes
#%END
#%option
#% key: weight
#% type: integer
#% description: Weight for the combination of flow accumulation and slope (number weights flow accumulation)
#% answer: 2
#% required : yes
#%END
#%option
#% key: smoothingtype
#% type: string
#% description: Type of smoothing to perform ("median" is recommended)
#% answer: median
#% options: average,median,mode,minimum,maximum
#% required : yes
#%END
#%option
#% key: smoothingsize
#% type: integer
#% description: Size (in cells) of the smoother's moving window
#% answer: 3
#% required : yes
#%END
#%flag
#% key: k
#% description: -k Keep the soil depth "potential" map (map name will be the same as specified in input option 'bedrock" with suffix "_depth_potential" appended to it)
#%END
#%flag
#% key: s
#% description: -s Print some useful statistics about the output soil depths (written to stdout)
#%END

import sys
import os
import tempfile
grass_install_tree = os.getenv('GISBASE')
sys.path.append(grass_install_tree + os.sep + 'etc' + os.sep + 'python')
import grass.script as grass


def main():
    soildepth = "temp.soildepth"
    tempsdepth = "temp.soildepth2"
    if options['soildepth'] is None:
        soildepth2 = "temp.soildepth3"
    else:
        soildepth2 = options['soildepth']
    elev = options['elevation']
    bedrock = options['bedrock']
    smoothingtype = options['smoothingtype']
    smoothingsize = options['smoothingsize']
    slopebreaks = options['slopebreaks'].split(";")
    flaccbreaks = options['flaccbreaks'].split(";")
    curvebreaks = options['curvebreaks'].split(";")
    weight = options['weight']
    smin = options['min']
    smax = options['max']
    slope = "tmp_slope_deletable"
    flacc = "tmp_flacc_deletable"
    pc = "tmp_pc_deletable"
    tc = "tmp_tc_deletable"
    mc = "tmp_mc_deletable"
    temprateflacc = "tmp_flacc_rate_deletable"
    temprateslope = "tmp_slop_rate_deletable"
    # let's grab the current resolution
    res = grass.region()['nsres']
    # make color rules for soil depth maps
    sdcolors = tempfile.NamedTemporaryFile()
    sdcolors.write(
        '100% 0:249:47\n20% 78:151:211\n6% 194:84:171\n0% 227:174:217')
    sdcolors.flush()
    grass.message('STEP 1: Calculating slope, flow accumulation, and curvatures\n')
    grass.run_command(
        'r.slope.aspect',
        quiet="True",
        overwrite="True",
        elevation=elev,
        slope=slope,
        pcurv=pc,
        tcurv=tc)
    grass.run_command(
        'r.watershed',
        quiet="True",
        overwrite="True",
        flags='a',
        elevation=elev,
        accumulation=flacc)
    # creat meancurvature map (in a manner compatible with older versions of
    # grass6), and then grab some stats from it for the later rescale
    # operation
    grass.mapcalc("${mc}=(${tc}+${pc})/2", quiet="True", mc=mc, tc=tc, pc=pc)
    mcdict = grass.parse_command('r.univar', flags="g", map=mc)
    # figuring out if user-supplied curvature breakpoints exceed actual limits
    # of curvature in the map, and adjusting if necessary
    if mcdict['min'] < curvebreaks[0].split(',')[0]:
        y1 = mcdict['min'] + ',' + curvebreaks[0].split(',')[1]
    else:
        y1 = curvebreaks[0]
    if mcdict['max'] > curvebreaks[1].split(',')[0]:
        y2 = mcdict['max'] + ',' + curvebreaks[1].split(',')[1]
    else:
        y2 = curvebreaks[1]

    grass.message(
        'STEP 2: Calculating "depth potential" across the landscape\n')
    # nested rescale of first flacc (to percentage of maximum soil depth
    # potential), and then curvature (to percentage offset from flac
    # function), and then combining the two measures to make final estimation
    # of soil depth potential. Final output depth potential map will be scaled
    # between 0 and 1, which maps to lowest depth potential to highest depth
    # potential.
    grass.mapcalc("${temprate}=eval(x=graph( ${flacc}, 0,0, 50,0.05,  250,0.25, 1000,1), y=graph( ${mc} , -0.0025,-0.5, 0,0, 0.0025,0.5), z=if(y < 0, x+(x*y), x+((1-x)*y)), if(z < 0, 0, if(z > 1, 1, z)))",
        quiet=True,
        temprate=temprateflacc,
        flacc=flacc, 
        mc=mc)
    # nested rescale of first slope (to percentage of maximum soil depth
    # potential), and then curvature (to percentage offset from slope
    # function), and then combining the two measures to make final estimation
    # of soil depth potential. Final output depth potential map will be scaled
    # between 0 and 1, which maps to lowest depth potential to highest depth
    # potential.
    grass.mapcalc(
        "${temprate}=eval(x=graph( ${slope}, 0,1 , ${x1}, ${x2}, 90,0), y=graph(${mc}, ${y1}, 0,0, ${y2}), z=if(y < 0, x+(x*y), x+((1-x)*y)), if(z < 0, 0, if(z > 1, 1, z)))",
        quiet="True",
        temprate=temprateslope,
        slope=slope,
        x1=slopebreaks[0],
        x2=slopebreaks[1],
        mc=mc,
        y1=y1,
        y2=y2)

    grass.message(
        'STEP 3: Calculating actual soil depths across the landscape (based on user input min and max soil depth values)\n')
    # create dictionary to record max and min rate so we can rescale it
    # according the user supplied max and min desired soil depths
    flaccratedict = grass.parse_command('r.univar', flags="g", map=temprateflacc)
    sloperatedict = grass.parse_command('r.univar', flags="g", map=temprateslope)
    # creating and running a linear regression to scale the calculated
    # landform soil depth potential into real soil depths using user specified
    # min and max soil depth values, and then average the two methods accroding
    # to the specified weighting of the Flacc method.
    grass.mapcalc(
        '${soildepth}=eval(x=graph(${temprateslope}, ${rmins},${smin}, ${rmaxs},${smax}), y=graph(${temprateflacc}, ${rminf},${smin}, ${rmaxf},${smax}), (x+(${weight}*y)/(1+${weight}) ))',
        quiet="True",
        soildepth=soildepth,
        temprateflacc=temprateflacc,
        temprateslope=temprateslope,
        rminf=flaccratedict['min'],
        rmaxf=flaccratedict['max'],
        rmins=sloperatedict['min'],
        rmaxs=sloperatedict['max'],
        smin=smin,
        smax=smax,
        weight=weight)
    unsmodict = grass.parse_command('r.univar', flags="g", map=soildepth)
    grass.run_command(
        'r.neighbors',
        quiet="True",
        input=soildepth,
        output=tempsdepth,
        size=smoothingsize,
        method=smoothingtype)
    # fix shrinking edge caused by neighborhood operation (and r.slope.aspect
    # above) by filling in the null areas. We do this by 0.98 * smax, since
    # the null cells are all the cells of the actual drainage divide with
    # slope basically = 0, and very mildly convex curvatures. This blends them
    # in nicely with the neighboring cells.
    grass.mapcalc(
        '${soildepth_real}=if(isnull(${input_sdepth}) && isnull(${elev}), null(), if(isnull(${input_sdepth}), 0.98*${smax},${input_sdepth}))',
        quiet='True',
        input_sdepth=tempsdepth,
        soildepth_real=soildepth2,
        elev=elev,
        smax=smax)
    # grab some stats if asked to
    if flags['s']:
        depthdict = grass.parse_command(
            'r.univar', flags="ge", map=soildepth2, percentile=90)

    grass.message('STEP 4: Calculating bedrock elevation map\n')
    grass.mapcalc(
        "${bedrock}=eval(x=(${elev} - ${soildepth}), if(isnull(x), ${elev}, x))",
        quiet="True",
        bedrock=bedrock,
        elev=elev,
        soildepth=soildepth)
    grass.message('Cleaning up...')
    grass.run_command(
        'g.remove', quiet=True, flags='f', type='raster',
        name=[
            pc,
            tc,
            mc,
            slope,
            flacc,
            soildepth,
            tempsdepth])
    if flags['k']:
        grass.run_command(
            'g.rename', quiet=True, raster="%s,%s_flacc_depth_potential" %
            (temprateflacc, bedrock))
        grass.run_command(
            'g.rename', quiet=True, raster="%s,%s_slope_depth_potential" %
            (temprateslope, bedrock))
    else:
        grass.run_command(
        'g.remove', quiet=True, flags='f', type='raster', name=[
            temprateflacc,
            temprateslope])
    if options['soildepth'] is None:
        grass.run_command(
        'g.remove', quiet=True, flags='f', type='raster', name=soildepth2)
    else:
        grass.run_command(
            'r.colors',
            quiet=True,
            map=soildepth2,
            rules=sdcolors.name)
    grass.message('DONE!\n')
    if flags['s']:
        grass.message(
            "min, max, and mean before smoothing: " +
            unsmodict['min'] +
            ", " +
            unsmodict['max'] +
            ", " +
            unsmodict['mean'])
        for key in depthdict.keys():
            grass.message('%s=%s' % (key, depthdict[key]))
        grass.message('Total volume of soil is %s cubic meters' %
                      (float(depthdict['sum']) * res * res))
    return

# Run above code.
if __name__ == "__main__":
    options, flags = grass.parser()
    main()
    exit(0)