GRASS r.horizon only works from 100–360 degrees

[twest820]

What is the bug or the crash?

What happens here is QGIS starts GRASS and then calls r.horizon on a sequence of sightline angles. If the angle is from 0 to 99 degrees the GRASS command lines fail because QGIS does not consistently zero pad the angle to three digits as GRASS does. Angles of 100–360 work as expected because they're large enough not to require zero padding.

For example, QGIS will send something like this to GRASS

g.proj -c wkt=".../crs.prj"
r.in.gdal input="...\DEM.tif" band=1 output="rast_632650235f80710" --overwrite -o
g.region n=225753.8777 s=190753.8777 e=135246.5282 w=102446.5282 res=100.0
r.horizon elevation=rast_632650235f80710 step=90 start=0 end=360 distance=1 output=output5fc5a6d0bd9f4e38a9632dfb72afa090 --overwrite
g.region raster=output5fc5a6d0bd9f4e38a9632dfb72afa090_0
r.out.gdal -t -m input="output5fc5a6d0bd9f4e38a9632dfb72afa090_0" output="...\000.tif" format="GTiff" createopt="TFW=YES,COMPRESS=LZW" --overwrite

and it will fail at g.region raster=output5fc5a6d0bd9f4e38a9632dfb72afa090_0

ERROR: Raster map or group <output5fc5a6d0bd9f4e38a9632dfb72afa090_0> not found

because the correct call is g.region raster=output5fc5a6d0bd9f4e38a9632dfb72afa090_000 (note the _0 versus _000 at the end of the raster name). The same layer name change is needed for the following r.out.gdal call's input parameter.

Steps to reproduce the issue

Select any DEM layer that's handy, launch r.horizon from the toolbox, set the angle step size to some positive number as required, and click run. Since the default is to output horizon angles on sightlines from 0 to 360 degrees any r.horizon execution will then hit this issue on at least the 0 degree step.

Versions

3.22.9

Supported QGIS version

• [X] I'm running a supported QGIS version according to the roadmap.

New profile

• [ ] I tried with a new QGIS profile

Additional context

I haven't tested with non-integer angles or non-integer steps. Probably those are also worth a look.

A basic workaround sequence is to

1. Use the automation QGIS provides to generate output rasters in the 100–360 degree range.
2. Copy QGIS's output trace and edit the 0–99 degree export commands QGIS generates to have the needed zero padding and suitable paths.
3. Set up a GRASS project and paste the setup commands, r.horizon call, and edited exports into GRASS's command window.

[agiudiceandrea]

@twest820, as a workaround, for "start angle" and "end angle" parameters values with or without decimals and only for "angle step size" parameter value without decimals, you can fix the issue changing the line 48 of r_horizon.py in your QGIS_install_dir\python\plugins\grassprovider\ext\ https://github.com/qgis/QGIS/blob/b19df398e8858c442706ffebeac16387d8fbf31b/python/plugins/grassprovider/ext/r_horizon.py#L48 to

        grassName = '{}_{}'.format('output{}'.format(alg.uniqueSuffix), '{0:0>3}'.format(int(num)))

[nicogodet]

With https://github.com/qgis/QGIS/blob/6f9da191b4d8edbea8e8cb3ab20b8acd199b4282/python/plugins/grassprovider/ext/r_horizon.py#L40-L51

def processOutputs(alg, parameters, context, feedback):
    # There will be as many outputs as the difference between start and end divided by steps
    start = alg.parameterAsDouble(parameters, 'start', context)
    end = alg.parameterAsDouble(parameters, 'end', context)
    step = alg.parameterAsDouble(parameters, 'step', context)
    num = start
    directory = alg.parameterAsString(parameters, 'output', context)
+   os.mkdir(directory)
    while num < end:
-       grassName = '{}_{}'.format('output{}'.format(alg.uniqueSuffix), int(num))
+       grassName = '{}_{}'.format('output{}'.format(alg.uniqueSuffix), '{0:0>3}'.format(int(num)))
        fileName = '{}.tif'.format(os.path.join(directory, '{0:0>3}'.format(int(num))))
        alg.exportRasterLayer(grassName, fileName)
        num += step

I'm not a big fan of os.mkdir and how it will be handle be temp processing folder cleaning by QGIS.

Another issue: it completely fails if direction param != 0

Calculating map 1 of 5 (angle 4.00, raster map ) 0..100 Calculating map 2 of 5 (angle 5.00, raster map ) 0..100 Calculating map 3 of 5 (angle 6.00, raster map ) 0..100 Calculating map 4 of 5 (angle 7.00, raster map ) 0..100 Calculating map 5 of 5 (angle 8.00, raster map ) 0..100 ERREUR : Carte raster non trouvée ERREUR : Raster map or group not found ERREUR : Carte raster non trouvée ERREUR : Raster map or group not found ERREUR : Carte raster non trouvée ERREUR : Raster map or group not found ERREUR : Carte raster non trouvée ERREUR : Raster map or group not found Checking GDAL data type and nodata value... 50..100 Using GDAL data type Exporting raster data to GTiff format...

[nicogodet]

I'm too slow...

[agiudiceandrea]

@nicogodet the "direction" parameter must be added to the angle value computed on each step. Something like:

    num = start + direction
    directory = alg.parameterAsString(parameters, 'output', context)
    while num < end + direction:
        grassName = '{}_{}'.format('output{}'.format(alg.uniqueSuffix), '{0:0>3}'.format(int(num)))
        fileName = '{}.tif'.format(os.path.join(directory, '{0:0>3}'.format(int(num))))
        alg.exportRasterLayer(grassName, fileName)
        num += step

[agiudiceandrea]

Anyway, those changes doesn't fix the algorithm if the "step" parameter value contains decimal digits, in which case the grass raster layer is named like output3348a4ab93d440929195546779451aca_AAA_BBBB with AAA always 3 digits, and BBBB whose number of digits depends on the "step" parameter number of decimal digits.

The exact logic is in: https://github.com/OSGeo/grass/blob/c726ae1a5f8e933c13746c1ad7fb236283559ca2/raster/r.horizon/main.c#L1160-L1176 https://github.com/OSGeo/grass/blob/c726ae1a5f8e933c13746c1ad7fb236283559ca2/lib/gis/basename.c#L51-L105 https://github.com/OSGeo/grass/blob/c726ae1a5f8e933c13746c1ad7fb236283559ca2/lib/gis/basename.c#L159-L184

[nicogodet]

@twest820 Are you able to test the linked PR ? Feedback on this would be nice

[twest820]

@nicogodet, did a quick test based on copying the diffs manually into my QGIS installation. Looks ok, thanks!

I think line 75 of the PR is normal numerical precision.