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.

  • Have you tried using forward slashes (or double quoting the backward slashes) in your path separators? – BradHards Aug 16 '16 at 4:08
  • you mean like this ? outputViewshed = 'C:/Users/nightswitch57/Desktop/Script/viewshed_%i.tif' % i – Ashan Lazarus Aug 16 '16 at 4:19
  • That is what I meant. – BradHards Aug 16 '16 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 – Ashan Lazarus Aug 22 '16 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 '16 at 5:09
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).
  • hey mate thank you for this, i'll give it a go in a weeks time when I renter the project. – Ashan Lazarus Feb 2 '17 at 2:31
  • @AshanLazarus any news? Does it work? – mgri Apr 26 '17 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. – Ashan Lazarus May 2 '17 at 0:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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