Sunday, 13 May 2018

qgis - Viewshed Plugin gives incorrect results


QGIS 2.18.5


The radius and positioning of my viewsheds seem to be doing something peculiar. I have set my CRS to British National Grid for both the terrain raster (downloaded from USGS) and my observation points. On the fly projection is unchecked. From the image below you see I have 2 observation points, one by the coast to the south, another to the north further inland.


enter image description here


My settings in the viewshed analysis plugin:



enter image description here


The results show both viewsheds out of place. The plugin thinks the observation points are some distance to the south of where they actually are. Not only that, but why is the radius of the calculated area egg shaped? Both of these questions I thought related to a CRS issue. I tried saving the layers to World Geodetic System 1984 and WGS 84. Both attempts led to a python error (the same CRS applied to each layer at the time). Looking at the uppermost viewshed, why isn't the entire valley shaded in? If I was standing at that point, I'm pretty sure I could see the ground around me...the plugin seems to be underestimating the results, or perhaps I'm using the plugin incorrectly.


enter image description here



Answer



I was able to reproduce this. I've never had this problem before, as I tend to use the OS OpenTerrain dataset for mapping in the UK.


I used SRTM and the wonderful Python Elevation module to download the same area (and stitch it into a geotiff) at the command line...


eio clip -o /tmp/cumbria.tif --bounds -3.85 53.94 -2.56 54.91

This closed issue on the plugin github suggests the plugin can't handle non-square pixels, and that's what is throwing out the locations, rather than a projection issue.


In my case, SRTM pixels, when projected into OSGB 27700 come out as non-square (18.5542m x -30.5976m)



Your best bet is to probably use the freely available OS Open Terrain 50 dataset.


Here's an example... you can see my original SRTM-projected-into-OSGB 15km large ellipse. Using the TomBio plugin I was able to locate Sellafield in grid square NY00. I opened that one tile, and did a 5km viewshed from 100m up on one of the building corners. Looks reasonable...


enter image description here


That data comes in small 10km x 10km (50m resolution) tiles in ASCII format, so you'll need to stitch the tiles together yourself (using Merge/VRT) for larger areas. But the cells are square.


EDIT You might also try gdalwarp with the -tr option to set the cell size to something more square, but not tried this. Might be worth trying before downloading a whole set of tiles :)


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...