1

I am using the Python simple editor in GRASS GIS to write a script that can import multiple raster files into GRASS GIS mapset. The first for loop successfully pulls in the orthophoto and its relevant bands into the mapset. Then the second for loop executes without error, but when I print the computation region extent info, this run_command with g.region did not result in my computational extent being updated. I was referencing this coding tutorial by Brendan Harmon. The difference here is that the orthophotos I am working with are found in different locations in a relatively small region.

This for loop successfully brings in the imagery, also set_resolution is a var set to 0.0068

# set path
data = os.path.join(
    gisdbase,
    r"C:\myfilepath")

# set region
gscript.run_command('g.region', res=set_resolution)

# iterate through files in directory
for file in os.listdir(data):
    filename = os.path.splitext(file)[0]
    # iterate through geotiffs
    if file.endswith('.tif'):
        # check if raster is already in mapset
        gscript.run_command('r.in.gdal',
            input=os.path.join(data, file),
            output=filename,
            overwrite=overwrite)
    else:
        pass 

But then when I run this for loop on the list of rasters in the mapset

    # list rasters in mapset
raster_list = gscript.list_grouped('rast',
    pattern='*ORTHO*')[mapset]

for raster in raster_list:
    # set region
    gscript.run_command(
        'g.region',
        raster=raster,
        res=set_resolution)

This is my output for g.region -p

zone:       0
datum:      wgs84
ellipsoid:  wgs84
north:      34:24:33.787198N
south:      34:24:30.533559N
west:       119:52:30.888396W
east:       119:52:21.920483W
nsres:      0:00:03.253639
ewres:      0:00:08.967913
rows:       1
cols:       1
cells:      1

Note that this is an easy task manually, but I would like to automate this. Since I am new to GRASS GIS, I would print my environment settings to give you all more context, but I am not sure what command that is.

1 Answer 1

0

A couple of things seem strange here.

First:

# set path
data = os.path.join(
    gisdbase,
    r"C:\myfilepath")

The gisdbase variable should point to your GRASS database. Why are you adding that drive and path to the gisdbase? You should not have any other files in the GRASS database.

Next:

raster_list = gscript.list_grouped('rast',
    pattern='*ORTHO*')[mapset]

Can we assume that the string "ORTHO" appears in all raster names? Is that from the original tif files?

Finally:

for raster in raster_list:
    # set region
    gscript.run_command(
        'g.region',
        raster=raster,
        res=set_resolution)

Are the extents or resolution of the original tifs different? I suggest to add a 'p' flag to the g.region command within the loop to get the region output for each raster.

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.