I need to do integrated viewshed analysis of a target with several observer locations. I am going to use r.viewshed analysis function in GRASS. With this function, it is possible to perform analysis based on single observer location. My aim is to add output of each single observer location into integrated viewshed output. Each output, i am going to flag with -b where o denotes non-visible while 1 denotes visible. Through a loop of observer locations, output of viewshed of each location can be integrated with r.mapcalc function. After integration, values = 0 in the rater map dentoes non-visible areas and value >= 1 denotes visible areas. Python code in GRASS for r.viewshed is ok. Under the loop of observer locations how to use r.mapcalc function with each viewshed adding? After getting the integrated viewshed raster, i need to extract its boundary depicting visible area in vector format. How to do that in GRASS? Any idea for this integrated analysis is highly appreciated.

If the integrated output is intoutput, can i use a code like this. in the loop if i=0 -> r.mapcalc 'intoutput = viewshedoutput' if i > 1 -> r.mapcalc 'intoutput = viewshedoutput + intoutput'


You could do that... or you could do it a less labour-intensive way in one step (instead of a step for every observer).

I assume you want to use free tools. However, in ArcGIS for instance you can use Spatial Analyst which sophisticated and allows setting of azimuths and even torus shaped viewing zones. You can also use the Observers module which will code your output raster according not only to how many observers can see the location, but specifically which ones can see it. However ArcGIS with Spatial Analyst is expensive.

I do a lot of this work and use a variety of tools but a nice free one is called OpenWind. OpenWind is designed for visibility calculations for windfarms, but if that is not what you are doing, don't worry, just fake a turbine to match your requirements (ie with a tip height of a human). When you use this tool (and others like it) you simply create a single shapefile of all your observer locations and OpenWind's ZTV (Zones of Theoretical Visibility) tool will tell you how many turbines can be seen at any one location (or in your case this equates to how many observers can see a given location).

Alternatively you could use SAGA GIS' viewshed tools. Again it is free and it can do integrated viewshed analysis. I'm less familiar with SAGA as other GIS I use have everything I need, but it has a good reputation.

  • Ok fine, i will look into OpenWind. However, if it is done with GRASS, how to extract vector boundary in the visible area of the analysis is somethong i need to know – Sanjeewa Jan 5 '13 at 14:47
  • 1
    You can use Raster->Conversion->Polygonize – MappaGnosis Jan 6 '13 at 15:36

This would be a lot easier to answer if you had posted an example of your code because I could set my edits within your existing loops using your existing object names. However, hopefully this should answer your specific question;

import grass as grass

# the below will convert 0s cells on your binary raster to null cells
# 'rastnulls' should be the name for this modified raster
# 'binaryrast' should be the name of the raster produced by r.viewshed

>> formula = '"%s" = if("%s" = 0, null(),1)' % (rastnulls, binaryrast)
>> grass.raster.mapcalc(formula, overwrite=True)

# then, run r.to.vect to convert these rasters to vector shapes
>> grass.run_command("r.to.vect", overwrite=True, input=rastnulls, output="losMap", feature="area")

Additionally, the GRASS 7.x plug-in r.viewshed.cva with the -k flag will produce multiple binary viewshed rasters for a set of points. You need to install it using the 'Install add-ons' option https://grass.osgeo.org/grass70/manuals/addons/r.viewshed.cva.html.

If you have trouble running the mapcalc script then insert this additional bit of code before running it.

formula = formula.replace('"', '')

Hope that helps.

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.