3

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.

7
  • Have you tried using forward slashes (or double quoting the backward slashes) in your path separators?
    – BradHards
    Aug 16, 2016 at 4:08
  • you mean like this ? outputViewshed = 'C:/Users/nightswitch57/Desktop/Script/viewshed_%i.tif' % i Aug 16, 2016 at 4:19
  • That is what I meant.
    – BradHards
    Aug 16, 2016 at 8:06
  • this produces the following error Traceback (most recent call last): File "<input>", line 1, in <module> File "c:/users/nights~1/appdata/local/temp/2/tmpwf9kd3.py", line 23 SyntaxError: Non-ASCII character '\xe2' in file c:/users/nights~1/appdata/local/temp/2/tmpwf9kd3.py on line 23, but no encoding declared; see python.org/peps/pep-0263.html for details Aug 22, 2016 at 4:22
  • That looks like a different problem - perhaps you want to close this question off and try posting a new question with the full command you are running, and the exact error.
    – BradHards
    Aug 22, 2016 at 5:09

1 Answer 1

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).
3
  • hey mate thank you for this, i'll give it a go in a weeks time when I renter the project. Feb 2, 2017 at 2:31
  • @AshanLazarus any news? Does it work?
    – mgri
    Apr 26, 2017 at 18:24
  • 1
    thanks for following up mate sorry this project has hit the back burner. i'll ry set up some time for it this week, I really do need to give it a shot. I sincerely thank you for helping me out with this. i'll keep you posted of the results when I get back into it. May 2, 2017 at 0:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.