4

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) Viewing Position data
  • 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)

DTM, Points A(black dots) and Points B (blue dots)

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: r.los

  • Do you know python? The best approach is to write a little python script, which passes one point geometry at a time – Curlew Oct 22 '13 at 9:34
5
+50

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.

  • Thanks afalciano. I would ask, is the 'max_dist' set to metres as default? And do I need to specify the exact path for the input/output files? I haven't used GRASS before you see. – Rob Lodge Oct 24 '13 at 11:53
  • Don't mention it. The max_dist parameter is expressed in meters, as you can see from the r.los manual. And you don't have to specify the file paths, but their location and mapset. So before starting you need to create a location and a mapset (e.g. LOS), then import your DTM. Finally, after running script.sh you have to export the results if you want to use them outside GRASS. You can easily accomplish all these tasks in QGIS using the GRASS plugin. – Antonio Falciano Oct 24 '13 at 12:49
4

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.

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!

  • I will have to repeat the requested operation a number of times, so ideally the method to obtain the lines of sight should be reproducible and as streamlined as possible. – Rob Lodge Oct 22 '13 at 11:42
  • I appreciate, for the LOS described above I've used it for 145 sites. While the processing/running time is lengthy, the script writing time is not. – MikeRSpencer Oct 24 '13 at 8:47

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.