Skip to content

[Bug] r.horizon creates strange offset if bufferzone is not multiple of raster resolution  #3266

@matzegoebel

Description

@matzegoebel

Describe the bug
When setting bufferzone to a value that is not a multiple of the resolution of the input raster, an offset between the input raster and the horizon angle may occur without an error or warning message. I used a raster with 3 m resolution. With bufferzone=5, the offset occurs, with bufferzone=6 it doesn't.
In my case, there is no offset in one part of the area (upper right corner). In the rest of the area there is an offset, but not everywhere in the same direction.

To Reproduce
Here is a minimal code example for Python:

from pathlib import Path
import os
from grass_session import Session
import grass.script as grass

tile = "test_3m"
angle = 90

base_dir = Path(__file__).parent
dsm_path = base_dir / f"{tile}.tif"
output_path = base_dir / f"horizon_{tile}_{angle:03d}.tif"
grass_db = base_dir / "grass_db"
grass_dir = grass_db / tile

# init GRASS
os.system(f"grass -c {dsm_path} -e {grass_dir}")
os.system(f"grass {grass_dir / 'PERMANENT/'} --exec r.external input={dsm_path} output=dsm")

with Session(
    gisdb=str(grass_db),
    location=tile,
    mapset=tile,
    create_opts="",
):
    # compute horizon angle
    grass.run_command(
        "r.horizon",
        elevation="dsm",
        start=angle,
        end=angle + 1,
        step=1,
        bufferzone=5,
        output=tile,
    )

    # write as geotiff
    grass.run_command(
        "r.out.gdal",
        input=f"{tile}_{angle:03d}",
        output=str(output_path),
        flags="c",
    )

Here is the input raster and the output that I got:
horizon_test_3m.zip

Expected behavior
There should not be an offset, no matter what the buffer size is, or an error should be raised stating that the buffer size must be a multiple of the raster resolution.

Screenshots
Here is a QGIS visualization of a part of the input raster and of the two horizon outputs with buffer size 5 and 6:
horizon_offset

System description:

  • Operating System: Ubuntu 22.04
  • GRASS GIS version: 7.8.7
  • Python: 3.11.5
  • details about further software components
    version=7.8.7
    date=2022
    revision=exported
    build_date=2022-02-23
    build_platform=x86_64-pc-linux-gnu
    build_off_t_size=8
    libgis_revision=2022-02-23T23:17:09+00:00
    libgis_date=2022-02-23T23:17:09+00:00
    proj=8.2.1
    gdal=3.4.1
    geos=3.10.2
    sqlite=3.37.2

I also tried with version 8.3.1, but couldn't get the code to run. It crashed with exit code -11. No idea why.

Metadata

Metadata

Assignees

No one assigned

    Labels

    bugSomething isn't working

    Type

    No type

    Projects

    No projects

    Milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions