I am trying to loop the r.sun function within a Python script.

Anyway, I seem not to understand how I can stack the raster maps each loop would produce into a single multi-band raster (preferably in geotiff format).

I think one way would be using the list produced with g.list or (even better) g.mlist as input parameter for the i.group function and then export the created group as multi-band file (as shown in the example here). Here is my simple code, modified from this page.

# import the needed modules
from osgeo import gdal
import numpy as np
import grass.script as grass

# prepare the needed input
dem = "dem4comparison@globrad_test"
doy = np.arange(1,366) # day of the year

# calculate slope and aspect maps
              elevation = dem,
              slope = "dem.slope@globrad_test",
              aspect = "dem.aspect@globrad_test")

# i is the day of year (DOY) to use for running r.sun
for d in doy:
    # want to get global radiation
    radname = "globrad" + str(d)
    # run the solar model
    grass.run_command("r.sun", flags = "s",
              elevin = dem,
              aspin = "dem.aspect@globrad_test",
              slopein = "dem.slope@globrad_test",
              insol_time = insolname, glob_rad = radname, day = d, step = 1)
    # export the outputs to GTiff files
    grass.run_command("r.out.gdal", input = radname + "@globrad_test", output = radname)

The last row would export each map to a geotiff file after each loop. Rather, I would stack them all into the same multi-band file with n*m*p dimensions, where m and n are rows and columns of my images (each of the same size of course) and p the number of output maps after the loops (that is, 365, the number of the days in a year).


Only a small modification is needed:

You need to generate the group prior to exporting (i.e., use i.group), then export to a multilayer file by assigning the group name as input in the r.out.gdal call.

  • fine thanks, but how do I go if I want to print a list of my maps created after the loop to be used inside i.group? As I said, I would like to do something like using the output from g.mlist as input for i.group. I red something about standard input and output, but I don't know if that's the way to follow. – umbe1987 Nov 3 '14 at 9:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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