r.viewshed.cva.py

#!/usr/bin/python
#
############################################################################
#
# MODULE:       	r.viewshed.cva.py
# AUTHOR(S):	Isaac Ullah
# PURPOSE:		Undertakes a "cumulative viewshed analysis" using a vector points map as input "viewing" locations, and utilizing the faster r.viewshed instead of r.los.
# COPYRIGHT:	(C) 2013 by Isaac Ullah
# REFERENCES:    r.viewshed
#			This program is free software under the GNU General Public
#			License (>=v2). Read the file COPYING that comes with GRASS
#			for details.
#
#############################################################################


#%Module
#%  description: Undertakes a "cumulative viewshed analysis" using a vector points map as input "viewing" locations, and the faster r.viewshed instead of r.los. NOTE: this routine requires the grass addon module r.viewshed, which can be added with g.extention
#%End

#%option
#% key: elev
#% type: string
#% gisprompt: old,cell,raster
#% description: Input elevation map (DEM)
#% required : yes
#%END 

#%option
#% key: output
#% type: string
#% gisprompt: old,cell,raster
#% description: Output CVA raster
#% required : yes
#%END

#%option
#% key: vect
#% type: string
#% gisprompt: old,vector,vector
#% description: Name of input vector points map containg the set of sites for this analysis.
#% required : yes
#%END

#%option
#% key: x_column
#% type: string
#% description: Column containing x values for site coordinates
#% required : yes
#%END

#%option
#% key: y_column
#% type: string
#% description: Column containing y values for site coordinates
#% required : yes
#%END

#%option
#% key: name_column
#% type: string
#% description: Column with unique identifiers for each point in input vector map (e.g. "cat" column)
#% required : yes
#%END

#%option
#% key: obs_elev
#% type: string
#% description: Height of observation points off the ground
#%answer: 0.0
#% required : yes
#%END

#%option
#% key: tgt_elev
#% type: string
#% description: Height of target areas off the ground
#%answer: 1.75
#% required : yes
#%END

#%option
#% key: max_dist
#% type: string
#% description: Maximum viewing distance (-1 = infinity)
#%answer: -1
#% required : yes
#%END

#%option
#% key: mem
#% type: string
#% description: Amount of memory to use (in MB)
#%answer: 1500
#% required : yes
#%END

#%option
#% key: refraction_coef
#% type: string
#% description: Refraction coefficient (with flag -r)
#%answer: 0.14286
#% required : no
#%END


#%flag
#% key: k
#% description: -k Keep all interim viewshed maps produced by the routine (maps will be named "vshed_'name'", where 'name' is the value in "name_column" for each input point)
#%END

#%flag
#% key: c
#% description: -c Consider the curvature of the earth (current ellipsoid)
#%END

#%flag
#% key: r
#% description: -r Consider the effect of atmospheric refraction
#%END

#%flag
#% key: b
#% description: -b   Output format is {0 (invisible) 1 (visible)}
#%END

#%flag
#% key: e
#% description:  -e   Output format is invisible = NULL, else current elev - viewpoint_elev
#%END



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


#main block of code starts here
def main():
    #bring in input variables
    elev = os.getenv("GIS_OPT_elev")
    out = os.getenv("GIS_OPT_output")
    vect = os.getenv("GIS_OPT_vect")
    xcol = os.getenv("GIS_OPT_x_column")
    ycol = os.getenv("GIS_OPT_y_column")
    namecol = os.getenv("GIS_OPT_name_column")
    obs_elev = os.getenv("GIS_OPT_obs_elev")
    tgt_elev = os.getenv("GIS_OPT_tgt_elev")
    max_dist = os.getenv("GIS_OPT_max_dist")
    mem = os.getenv("GIS_OPT_mem")
    refraction_coef = os.getenv("GIS_OPT_refraction_coef")
    #assemble flag string
    if os.getenv('GIS_FLAG_r') == '1':
        f1 = "r"
    else:
        f1 = ""
    if os.getenv('GIS_FLAG_c') == '1':
        f2 = "c"
    else:
        f2 = ""
    if os.getenv('GIS_FLAG_b') == '1':
        f3 = "b"
    else:
        f3 = ""
    if os.getenv('GIS_FLAG_e') == '1':
        f4 = "e"
    else:
        f4 = ""
    flagstring = f1 + f2 + f3 +f4
    #read in info from the table of the vector sites map, and parse it into a list of lists of info for each site
    s1 = grass.read_command("v.db.select", quiet = "True", map = vect, columns = xcol + "," + ycol + "," + namecol, fs = ",", nv = "False").strip()
    masterlist = []
    for item in s1.split("\n"):
        masterlist.append(item.strip("\n").split(","))
    #the first row is the column names, so pop that out of our master list
    index = masterlist.pop(0)
    #now, loop through the master list and run r.viewshed for each of the sites, and append the viewsheds to a list (so we can work with them later)
    vshed_list = []
    for site in masterlist:
        grass.message('Calculating viewshed for location %s,%s (point name = %s)\n' % (site[0], site[1], site[2]))
        tempry = "vshed_%s" % site[2]
        vshed_list.append(tempry)
        grass.run_command("r.viewshed", quiet = "True",  flags = flagstring,  input = elev, output = tempry, coordinate = site[0] + "," + site[1], obs_elev = obs_elev, tgt_elev = tgt_elev, max_dist = max_dist, mem = mem,  refraction_coef = refraction_coef)
    #now make a mapcalc statement to add all the viewsheds together to make the outout cumulative viewsheds map
    grass.message("Calculating \"Cumulative Viewshed\" map")
    #grass.mapcalc("${output}=${command_string}", quiet = "True", output = out, command_string = ("+").join(vshed_list))
    grass.run_command("r.series", quiet = "True", input = (",").join(vshed_list), output = out, method = "count")
    #Clean up temporary maps, if requested
    if os.getenv('GIS_FLAG_k') == '1':
        grass.message("Temporary viewshed maps will not removed")
    else:
        grass.message("Removing temporary viewshed maps")
        grass.run_command("g.remove",  quiet = "True",  rast = (",").join(vshed_list))
    return
    
# here is where the code in "main" actually gets executed. This way of programming is neccessary for the way g.parser needs to run.
if __name__ == "__main__":
    if ( len(sys.argv) <= 1 or sys.argv[1] != "@ARGS_PARSED@" ):
        os.execvp("g.parser", [sys.argv[0]] + sys.argv)
    else:
        main()