r.sun.daily: add average method (#933)

src/raster/r.sun.daily/r.sun.daily.html

@@ -8,7 +8,7 @@ <h3>Output parameters explanation</h3>
 There are two basic options:
 <ul>
 <li>output series of maps (one for each day): options containing basename in their name</li>
-<li>output one map which is a sum of the intermediate maps</li>
+<li>output one map which is an aggregation of the intermediate maps</li>
 </ul>
 
 You can choose any combination of parameters: e.g. total map of diffuse radiance and

src/raster/r.sun.daily/r.sun.daily.py

@@ -154,7 +154,7 @@
 # %option
 # % key: step
 # % type: double
-# % description: Time step when computing all-day radiation sums [decimal hours]
+# % description: Time step when computing all-day radiation [decimal hours]
 # % answer: 0.5
 # %end
 
@@ -168,39 +168,39 @@
 # % key: beam_rad
 # % type: string
 # % gisprompt: new,cell,raster
-# % description: Output beam irradiation raster map cumulated for the whole period of time [Wh.m-2.day-1]
+# % description: Output beam irradiation raster map aggregated for the whole period of time [Wh.m-2.day-1]
 # % required: no
 # %end
 
 # %option
 # % key: diff_rad
 # % type: string
 # % gisprompt: new,cell,raster
-# % description: Output diffuse irradiation raster map cumulated for the whole period of time [Wh.m-2.day-1]
+# % description: Output diffuse irradiation raster map aggregated for the whole period of time [Wh.m-2.day-1]
 # % required: no
 # %end
 
 # %option
 # % key: refl_rad
 # % type: string
 # % gisprompt: new,cell,raster
-# % description: Output ground reflected irradiation raster map cumulated for the whole period of time [Wh.m-2.day-1]
+# % description: Output ground reflected irradiation raster map aggregated for the whole period of time [Wh.m-2.day-1]
 # % required: no
 # %end
 
 # %option
 # % key: glob_rad
 # % type: string
 # % gisprompt: new,cell,raster
-# % description: Output global (total) irradiance/irradiation raster map cumulated for the whole period of time [Wh.m-2.day-1]
+# % description: Output global (total) irradiance/irradiation raster map aggregated for the whole period of time [Wh.m-2.day-1]
 # % required: no
 # %end
 
 # %option
 # % key: insol_time
 # % type: string
 # % gisprompt: new,cell,raster
-# % description: Output insolation time raster map cumulated for the whole period of time [h]
+# % description: Output insolation time raster map aggregated for the whole period of time [h]
 # % required: no
 # %end
 
@@ -235,7 +235,7 @@
 # %option
 # % key: insol_time_basename
 # % type: string
-# % label: Base name for output insolation time raster map cumulated for the whole period of time [h]
+# % label: Base name for output insolation time raster map aggregated for the whole period of time [h]
 # % description: Underscore and day number are added to the base name for daily maps
 # %end
 
@@ -248,6 +248,16 @@
 # % description: If not specified, r.sun default will be used.
 # %end
 
+# %option
+# % key: method
+# % type: string
+# % required: no
+# % multiple: no
+# % label: Method for daily maps aggregation
+# % options: sum,average
+# % answer: sum
+# %end
+
 # %option
 # % key: nprocs
 # % type: integer
@@ -431,13 +441,23 @@ def check_daily_map_names(basename, mapset, start_day, end_day, day_step):
             )
 
 
-def sum_maps(sum_, basename, suffixes):
+def maps_sum(out_, basename, suffixes):
     """
     Sum up multiple raster maps
     """
     maps = "+".join([basename + suf for suf in suffixes])
+    grass.mapcalc(f"{out_} = {maps}", overwrite=True, quiet=True)
+
+
+def maps_avg(out_, basename, suffixes):
+    """
+    Get average value from multiple raster maps.
+    """
+    maps = "+".join([basename + suf for suf in suffixes])
     grass.mapcalc(
-        "{sum_} = {new}".format(sum_=sum_, new=maps), overwrite=True, quiet=True
+        f"{out_} = ({maps}) / {len(suffixes)}",
+        overwrite=True,
+        quiet=True,
     )
 
 
@@ -458,6 +478,7 @@ def main():
     albedo_value = options["albedo_value"]
     horizon_basename = options["horizon_basename"]
     horizon_step = options["horizon_step"]
+    method = options["method"]
 
     # outputs
     beam_rad = options["beam_rad"]
@@ -492,7 +513,11 @@ def main():
     end_day = int(options["end_day"])
     day_step = int(options["day_step"])
 
-    if day_step > 1 and (beam_rad or diff_rad or refl_rad or glob_rad or insol_time):
+    if (
+        day_step > 1
+        and method == "sum"
+        and (beam_rad or diff_rad or refl_rad or glob_rad or insol_time)
+    ):
         grass.fatal(
             _("Day step higher then 1 would produce" " meaningless cumulative maps.")
         )
@@ -649,16 +674,23 @@ def main():
             # Empty process list
             proc_list = []
 
+    # process multiple maps into the desired one
+    if method == "average":
+        stats_func = maps_avg
+    else:
+        # sum (original behavior) as a fallback
+        stats_func = maps_sum
+
     if beam_rad:
-        sum_maps(beam_rad, beam_rad_basename, suffixes_all)
+        stats_func(beam_rad, beam_rad_basename, suffixes_all)
     if diff_rad:
-        sum_maps(diff_rad, diff_rad_basename, suffixes_all)
+        stats_func(diff_rad, diff_rad_basename, suffixes_all)
     if refl_rad:
-        sum_maps(refl_rad, refl_rad_basename, suffixes_all)
+        stats_func(refl_rad, refl_rad_basename, suffixes_all)
     if glob_rad:
-        sum_maps(glob_rad, glob_rad_basename, suffixes_all)
+        stats_func(glob_rad, glob_rad_basename, suffixes_all)
     if insol_time:
-        sum_maps(insol_time, insol_time_basename, suffixes_all)
+        stats_func(insol_time, insol_time_basename, suffixes_all)
 
     # FIXME: how percent really works?
     # core.percent(1, 1, 1)