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 grass.run_command("r.slope.aspect", 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
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).