0

I'm new to GRASS and struggling to do a viewshed across multiple tiles. I have got the GRASS r.viewshed function working from Python as follows for a single tile:

gis_path = os.path.join(output_dir, 'gisdb')
output_filename = '{}-viewshed.tif'.format(filename)
transmitter_coords = str(x_transmitter) + ',' + str(y_transmitter)

with Session(gisdb=gis_path, location="location", create_opts="EPSG:27700"):

    print(gcore.parse_command("g.gisenv", flags="s"))
    gcore.run_command('r.external', input=tile_path, output=filename, overwrite=True)
    gcore.run_command('r.external.out', directory=output_dir, format="GTiff")
    gcore.run_command('g.region', raster=filename)
    gcore.run_command('r.viewshed',flags='b',input=filename,output=output_filename,coordinates=transmitter_coords,overwrite=True)
    gcore.run_command('r.external.out', flags="r")

What is the right way to do this for multiple tiles? I have the following attempt with a loop that correctly reads in a list of file paths:

with Session(gisdb=gis_path, location="location", create_opts="EPSG:27700"):

    print(gcore.parse_command("g.gisenv", flags="s"))

    for tile in list_of_tiles:

        output_filename = '{}-viewshed.tif'.format(tile[-10:-4])
        gcore.message('Importing %s -> %s@%s...' % (tile, tile[-10:], output_dir))
        gcore.run_command('r.external', input=tile, output=output_filename, overwrite=True)
        gcore.run_command('r.external.out', directory=output_dir, format="GTiff")

I'm unsure if I now need to merge together the individual files (not sure how to do this) or to feed r.viewshed a list (which I haven't been able to do).

2 Answers 2

1

Yes, you should patch together all tiles and then do r.viewshed on the patched raster. Here is what it might look like:

import os
with Session(gisdb=gis_path, location="location", create_opts="EPSG:27700"):
    print(gcore.parse_command("g.gisenv", flags="s"))
    for tile in list_of_tiles:
        base, ext = os.path.splitext(os.path.split(tile)[1])
        tile_name = "tile_%s" % base
        gcore.run_command('r.import', input=tile, output=tile_name, overwrite=True)

    rast_list = gcore.read_command('g.list', type='rast', pattern="tile_*", separator="comma").strip().split('\n')
    gcore.run_command('g.region', rast=rast_list, flags='p')
    gcore.run_command('r.patch', input=rast_list, output="all_tiles")
    gcore.run_command('r.viewshed',flags='b',input=all_tiles, output="viewshed", coordinates=transmitter_coords, overwrite=True)

Now export "viewshed" to Geotiff with r.out.gdal or r.external.out as you did above.

2
  • Dear Micha, this partially worked but do you know why my rast_list returns this output: b'tile_TL45NE,tile_TL45NW' ? I am unsure of by the 'b' is included? Dec 28, 2018 at 15:48
  • What version of python? The b'...' indicates a byte string. You should be able to just to decode('UTF-8') on the rast_list string.
    – Micha
    Dec 29, 2018 at 16:34
0

For an easy approach, you can use r.buildvrt in order to build a VRT (Virtual Raster) from the list of input raster maps. It works similar to gdalbuildvrt but straight with GRASS GIS raster maps. On the resulting virtual map just run r.viewshed.

Note that you need to use at least GRASS GIS 7.6 in order to have r.buildvrt. The official release shall be in the next weeks but it is already available as a weekly snapshot from https://grass.osgeo.org/grass76/

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.