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

Question

I wish to carry out a line of site analysis similar to the subject in here (preferably within QGIS): Visual Impact Assessment in QGIS or GRASS

Scenario:

I have two points, A and B. Point A is a signal provider and Point B is a signal receiver. Point A is at X altitude and Point B is at Y altitude. There may be hills between Points A and B (represented by either vector contour lines or a raster DTM) that block Point B from receiving signal.

Data:

• A vector shapefile containing hundreds of points (the viewing position / Point A)
• A vector shapefile containing several thousands of points (each point in this shapefile represents a position that needs to be viewed / Point B)
• A 50 m resolution DTM in raster format (currently looking for higher resolution free data for the UK)

Question:

How can I establish that Point As are able to view all of Point Bs (via say, a mask?). And how do I run this as a batch process? At the moment I can only run the line of sight (LOS) analysis for one coordinate at a time:

Answer 1

You can make the preliminary tasks of your visibility analysis semi-automatic in few steps. Firstly, convert the shapefile of the observation points (Final_Branches.shp) in the CSV format, exporting only the coordinates:

ogr2ogr -f CSV coords.csv Final_Branches.shp -sql "SELECT Easting, Northing FROM Final_Branches"

The content of coords.csv should be:

Easting,Northing
255745,131455
245819,130162
247348,130341
...and so on

Then, define a batch file (e.g. create_script.bat) to loop through each line of the CSV file (skipping the header) echoing r.los (copying and pasting directly from GRASS, as Mike stated before, and setting the coordinate parameter equal to the variable %%i, i.e. the X,Y to be read from the CSV file):

@ECHO OFF    
for /F "skip=1 tokens=*" %%i in (coords.csv) DO (
ECHO r.los input=DTM@LOS output=visible_area coordinate=%%i max_dist=10000
)

Redirect the output to another file, e.g. script.sh:

create_script.bat > script.sh 

Opening script.sh in a text editor you should have something like:

r.los input=DTM@LOS output=visible_area coordinate=255745,131455 max_dist=10000
r.los input=DTM@LOS output=visible_area coordinate=245819,130162 max_dist=10000
r.los input=DTM@LOS output=visible_area coordinate=247348,130341 max_dist=10000
...and so on

Finally, execute sh script.sh from the GRASS command line.

Answer 2

For the record, a QGIS Python console solution for generating multiple viewsheds:

import processing
from PyQt4.QtCore import QFileInfo

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

i=0
for ft in pointLayer.getFeatures():
    point = ft.geometry().asPoint()
    coordStr = '%d,%d' % (point.x(),point.y())
    outputViewshed = 'd:/temp/viewshed_%i.tif' % i
    outputs_0=processing.runalg("grass:r.los", rasterLayer, coordStr, None, 10, 10000, True, '700000.0,710000.0,6200000.0,6220000.0', 30, outputViewshed)
    i = i + 1

You can then on the signal receiver layer (layer B) use the point sampling tool plugin, to create a new point table with compiled viewshed data. The r.los does not take a parameter for looking at height above DEM, only looking from height. So your receiver is at ground level i.e. the DEM value. Converting the receiver point layer with a height column, to raster and adding receiver raster to the DEM could fix that problem.

Answer 3

Avoiding python...

You can create a very crude script in GRASS using a spreadsheet. In GRASS itself (not GRASS via QGIS) when you set up r.los (the line of sight analysis) it will generate code for the process at the bottom of the window, it will look like this:

r.los input=dtm@LOS output=visible_area coordinate=XXXX,YYYY max_dist=15000

Copy this into a spreadsheet and split it into separate columns, like:

col1:col2:col3:col4:col5:col6:col7

r.los input=1@LOS output=:test: coordinate=:XXXX:,:YYYY: max_dist=15000

You can then paste in your columns of x and y values, use fill down on the other columns and finally concatenate() to join it all together. Paste the concatenated column into a text file and then run it in GRASS from the GRASS terminal window. This should be something like:

Linux

sh ~/filepath/file.sh

Windows

sh C:\filepath\file.sh

As mentioned initially, this method is crude and a very long way from pretty - but it does work!