Programmatically performing viewshed on points from a shapefile?

Question

I want to perform a viewshed analysis on multiple points in a shapefile. I'm working with a solution as seen in this post by Jakob

Viewshed / Line of Sight analysis- as a batch process?

I modified his code such that I am able to pass the extents into r.viewshed. this is given by the variable "extcoordStr" in the below code

import processing
from PyQt4.QtCore import QFileInfo

rasterLayer = qgis.utils.iface.activeLayer()
      pointLayer=QgsMapLayerRegistry.instance().mapLayersByName('Tasmania_select_two')[0]

i=0
for ft in pointLayer.getFeatures():
   point = ft.geometry().asPoint()
   coordStr = '%d,%d' % (point.x(),point.y())
   coordArray = coordStr.split(',')

   #extent generation, +- 50,000Metres from point origin(coordStr)
   xMin = int(coordArray[0]) - 50000
   xMax = int(coordArray[0]) + 50000
   yMin = int(coordArray[1]) - 50000
   yMax = int(coordArray[1]) + 50000

   #Type juggling to get xmin,xmax,ymin,ymax
   extcoordStr = str(xMin) + ',' + str(xMax) + ',' + str(yMin) + ',' + str(yMax)

   outputViewshed = 'C:\Users\nightswitch57\Desktop\Script\viewshed_%i.tif' % i

   #running viewshed with observer and target heights set to 20Metres at max distance of 50Km
   outputs_0=processing.runalg("grass7:r.viewshed", rasterLayer,coordStr,20,20,50000,False, extcoordStr,0,outputViewshed)
   i = i + 1

problem when I click run, my computer freezes up a bit task manager shows GRASS is been used but at the end there is nothing in the output folder. in Jakobs answer he uses, I feel this might be the reason why I do not have any output in the folder I specified

outputViewshed = 'd:/temp/viewshed_%i.tif' % i

while I use

outputViewshed = 'C:\Users\nightswitch57\Desktop\Script\viewshed_%i.tif' % i

Also in addition I would much rather have the viewshed named after the point features "Site ID" which is in a column within the point shapefile (Tasmania_select_two) as currently it is supposed to name each viewshed as the number of the point feature. to achieve this i'm guessing I would need to append the point features Site ID to the outputViewshed variable somehow but I don't know how to.

Answer 1

I hope it's not too late for answering.

Your solution won't work because you reproduced the syntax of the r.los algorithm, but it is now deprecated: you probably knew that because you called the r.viewshed algorithm, but you forgot looking for the right way of running it.

In addition to this, you didn't see any error because you disabled the Python Console. In fact, I ran your code and this message appeared in the Python Console every time there was a wrong call to the algorithm (i.e. equal to the number of points):

Error: Wrong number of parameters
ALGORITHM: r.viewshed - Computes the viewshed of a point on an elevation raster map.
    input <ParameterRaster>
    coordinates <ParameterString>
    observer_elevation <ParameterString>
    target_elevation <ParameterString>
    max_distance <ParameterString>
    refraction_coeff <ParameterNumber>
    memory <ParameterNumber>
    -c <ParameterBoolean>
    -r <ParameterBoolean>
    -b <ParameterBoolean>
    -e <ParameterBoolean>
    GRASS_REGION_PARAMETER <ParameterExtent>
    GRASS_REGION_CELLSIZE_PARAMETER <ParameterNumber>
    output <OutputRaster>

The above lines are self-explanatory since they take in account many parameters that you neglected (if you want to know what these new parameters mean, please refer to the documentation of the r.viewshed algorithm).

I edited your code and it works by now (remember to adapt the parameters to your specific case):

import processing
from PyQt4.QtCore import QFileInfo

rasterLayer = qgis.utils.iface.activeLayer()
pointLayer=QgsMapLayerRegistry.instance().mapLayersByName('Tasmania_select_two')[0]

i=0
for ft in pointLayer.getFeatures():
    point = ft.geometry().asPoint()
    coordStr = '%d,%d' % (point.x(),point.y())
    coordArray = coordStr.split(',')

    #extent generation, +- 50,000Metres from point origin(coordStr)
    xMin = int(coordArray[0]) - 50000
    xMax = int(coordArray[0]) + 50000
    yMin = int(coordArray[1]) - 50000
    yMax = int(coordArray[1]) + 50000

    #Type juggling to get xmin,xmax,ymin,ymax
    extcoordStr = str(xMin) + ',' + str(xMax) + ',' + str(yMin) + ',' + str(yMax)

    outputViewshed = 'C:/Users/nightswitch57/Desktop/Script/viewshed_%i.tif' % i

    #running viewshed with observer and target heights set to 20Metres at max distance of 50Km
    outputs_0=processing.runalg('grass7:r.viewshed', rasterLayer, coordStr, '0', '0', '800', 0.14286, 500, False, False, False, False, extcoordStr, 0, outputViewshed)
    i = i + 1

Two final notes:

1. You need to use slashes when setting the output path, otherwise the above code won't work;
2. I posted a working solution for the code which you proposed. If you want to rename your output using a value from a field instead of a number, you may use this line: outputViewshed = 'C:/Users/nightswitch57/Desktop/Script/viewshed_%s.tif' % (ft["ID_IN"]) (the changes are in the using of %s and the name of the field).