0

I'm trying to use r.viewshed in PyQGIS to perform Viewshed analysis on a point placed on a raster layer. The algorithm seems to run correctly, as no errors come up, but the output .tif file is not created, and I've tried it using multiple different locations, using absolute and relative paths, both with no luck.

Here is a snippet of my python code below.

rasterLayer is the QgsRasterLayer object, the main raster of the project for the viewshed to be done upon.

pointPath =  r"C:/Users/me/Documents/ZVIs/zvi.tif"
maxRange = 500
pointElevationFromDTM = 100
coordStr = "230049.90500026295, 2841692.311715042"
result = processing.run("grass7:r.viewshed",
                           {'input' : rasterLayer,
                            'coordinates': coordStr,
                            'observer_elevation': pointElevationFromDTM,
                            'max_distance': maxRange,
                            'memory': 300,
                            'output': pointPath})
            
if result['output']:
    print("Viewshed analysis completed. Output saved to:", pointPath)
else:
    print("Viewshed analysis failed")
if os.path.isfile(pointPath):
    print(f"ZVI file was created.")
else:
    print(f"ZVI file was not created")

I added in some checks at the end so you can see what I mean. It seems that an output is being created, just not saved for whatever reason.


I've been able to print out the processing feedback after the answer by @Taras. The main errors I can see are below:

ami_single_temp_name: G_mkstemp() failed: : No such file or directory

Assertion failed: 0, file ami_stream.cpp, line 91 

ERROR: Raster map <outputb422948942aa4ef6a24c5ce880e5fe23> not found

ERROR: Raster map or group <outputb422948942aa4ef6a24c5ce880e5fe23> not found

But when using isValid() on the raster layer it is returning True.

Where do I go next?

I've googled the errors but haven't had any luck yet.

4
  • What does the log say? Commented Feb 22 at 10:40
  • How do I check the log? If you mean the console output, it simply says that the viewshed analysis completed, but the ZVI file was not created. This is just output from the checks at the bottom of my code however.
    – oliver
    Commented Feb 22 at 10:42
  • I can confirm coordStr lies within the rasterLayer. The CRS of the rasterLayer is EPSG:32615 - WGS 84 / UTM zone 15N. How would I know if the rasterLayer CRS matches the coordStr? I've just ran it with a maxRange of 10, which should definitely lie within the rasterLayer, and it gave the same issues.
    – oliver
    Commented Feb 22 at 11:19
  • I don't think I can as it's private data, but I ran the same code on the publicly available raster SR_50M_alaska_nad.tif found at https://download.qgis.org/downloads/data/, and it is behaving exactly the same, when tested with coordStr = "1100248.3159823148, 6360501.794258764", with pointElevationFromDTM =100 and maxRange=10. Hope this helps.
    – oliver
    Commented Feb 22 at 12:20

1 Answer 1

2

There are several suggestions regarding your code:

  • check firstly alg input and output by means of the processing.algorithmHelp("grass7:r.viewshed")

  • as was suggested by @bugmenot123, check the logs. One can access them via putting the QgsProcessingFeedback() class into processing.run(...) and after applying the textLog()-method:

    feedback = QgsProcessingFeedback()
    processing.run("grass7:r.viewshed", params, feedback=feedback)
    print(feedback.textLog())
    
  • use the Python Try Except blocks to handle errors and exceptions, e.g. QgsProcessingException.

    try:
       params = {
          'input': rlayer,
          'coordinates': coordStr,
          'observer_elevation': pntElevationFromDTM,
          'max_distance': maxRange,
          'output': outputPath
       }
       processing.run("grass7:r.viewshed", params, feedback=feedback)
    except QgsProcessingException as e:
       print(e)
    
  • apply some useful functions e.g. basename, splitext, normpath, isfile, dirname, and join from the os.path library

  • use the isValid()-method to check the status of the QgsRasterLayer

So, in the end, your code may look like this (Note: it was tested within a standalone application):

# imports
from sys import path
from qgis.core import (QgsApplication, QgsRasterLayer, QgsProcessingFeedback, QgsProcessingException)
from os.path import (basename, splitext, normpath, isfile, dirname, join)

# See https://gis.stackexchange.com/a/155852/4972 for details about the prefix
QgsApplication.setPrefixPath('C:\\OSGeo4W\\apps\\qgis-ltr\\plugins', True)
qgs = QgsApplication([], False)
qgs.initQgis()

# Append the path where processing plugin can be found
path.append('C:\\OSGeo4W\\apps\\qgis-ltr\\python\\plugins')

import processing
from processing.core.Processing import Processing

Processing.initialize()

raster_path = normpath("D:\\qgis_test\\qgis_sample_data\\raster\\SR_50M_alaska_nad.tif")

if isfile(raster_path):

    rlayer = QgsRasterLayer(raster_path, splitext(basename(raster_path))[0])
    coordStr = "3178941, 3504031"
    pntElevationFromDTM = 65  # draped from the raster
    maxRange = 50000
    outputPath = join(dirname(raster_path), "AlaskaViewshed.tif")

    if rlayer.isValid():
        feedback = QgsProcessingFeedback()

        try:
            params = {
                'input': rlayer,
                'coordinates': coordStr,
                'observer_elevation': pntElevationFromDTM,
                'max_distance': maxRange,
                'output': outputPath
            }
            processing.run("grass7:r.viewshed", params, feedback=feedback)
            print(feedback.textLog())  # read the logs

        except QgsProcessingException as e:
            print(e)

With the above code, I could get the output for the SR_50M_alaska_nad.tif raster on QGIS version 3.28.10-Firenze.

output


References:

3
  • I thought the issue was fixed but it seems to still be happening. I've printed the feedback to console, and the main errors seem to be: ami_single_temp_name: G_mkstemp() failed: : No such file or directory, Assertion failed: 0, file ami_stream.cpp, line 91, ERROR: Raster map <outputb422948942aa4ef6a24c5ce880e5fe23> not found, ERROR: Raster map or group <outputb422948942aa4ef6a24c5ce880e5fe23> not found I've tested and the raster returns True from isValid() though. Not sure what the next steps would be
    – oliver
    Commented Mar 4 at 15:57
  • It seems like my code works some times and other times it does not. Yesterday it was giving those errors, and I haven't changed it, but today it's working fine again. isfile(raster_path) returns true at the moment, but I'll let you know what it says if the code stops working again. Really not sure why it's working again today when it wasn't yesterday though.
    – oliver
    Commented Mar 5 at 8:37
  • 1
    That's what I'm thinking. Although with the errors above it seemed unable to find the raster? Only way to really know is to wait for it to happen again. I'll let you know if it does.
    – oliver
    Commented Mar 5 at 9:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.