Finding point with largest viewshed using QGIS?

Question

I have a large number of points (2600) and I need to rank them by the size of their viewshed.

I have been using the Viewshed Analysis plugin in QGIS, but I can only select the option for a cumulative raster of all the viewsheds together. In the result raster, its not possible to determine the size of the viewedsheds for the individual points.

Answer 1

there are two parts to this.. creating individual viewshed rasters, and doing the report. I'll assume all your points are in one layer.

---

One raster per observer

You can at least create one binary raster per feature (observer) using the Advanced Viewshed Analysis plugin. The settings to use are

• output file: specify a pattern like /tmp/raster
• output: binary viewshed
• make sure cumulative is NOT checked

That should create one raster per feature (/tmp/raster_1_Binary, /tmp/raster_2_Binary and so on). Each raster will be the same size as your source DEM, and should be 1=visible, 0=not visible, Nodata = outside radius.

I found this by accident, I forgot to set 'cumulative', left it running, and had to tidy up thousands of files afterwards ;-)

The next part is to count the number of cells in each raster which are 1.

---

Counting visible cells

There are a couple of GRASS tools, r.stats and r.report. I think for multiple rasters these create cross-tabulations, which I don't think are suitable here... you want to consider each raster in isolation. (Someone correct me if I'm wrong)

Using gdalinfo from command line

The following bash shell code should work on Mac or Linux - I'm sure something equivalent will be possible on Windows, maybe using Powershell?

The units will be in cells, you'd need to convert areas yourself if you want it in square meters.

It makes use of the fact that gdalinfo can dump the histogram as a single line of 256 numbers separated by spaces. This scripts finds the histogram line, and picks out the first and last numbers (the others are all 0). This is a really dodgy way to do this, and I don't know how to change the number of "bins". But it worked on the output from the plugin:-

#!/bin/bash
for i in $( ls /tmp/raster_*_Binary.tiff); do
    echo $i `gdalinfo -mm -hist $i | grep -E '^(\s*[0-9]+\s+)*$' | cut -d ' ' -f 3,258`
done

That outputs three items per line; filename, number of 0 cells, number of 1 cells. If you don't have any spaces in your file paths, you can feed that into Excel. Example output:-

/Users/steven/data/Untitled_1_Binary.tiff 115503 10106
/Users/steven/data/Untitled_2_Binary.tiff 116432 9177
/Users/steven/data/Untitled_3_Binary.tiff 114720 10889

These two numbers in each file add up to 125609, corresponding to a 200m radius circle with my test rasters.