[GRASS-dev] Linke turbidity in r.sun.daily
Nikos Alexandris
nik at nikosalexandris.net
Fri Mar 27 14:14:12 PDT 2015
* Vaclav Petras <wenzeslaus at gmail.com> [2015-03-26 23:12:50 -0400]:
> On Thu, Mar 26, 2015 at 5:04 AM, Nikos Alexandris <nik at nikosalexandris.net>
> wrote:
> > Disregard the previous one. This one without old 'rast' entries (instead,
> > new 'raster').
> >
> > Looks OK (I'm only looking at the source code). Just few things:
^^^--- please check your mail clients settings -- first line of your
replies get's added to the original sender. I'd like to know if it's
something messy in my side!
)
> It seems you removed the option to add timestamps without registering to
> temporal database. I would just leave the else there. (You may want to
> create timestamps but not register or register later.)
Misjudgement, lead by the comment "...or something in GRASS6.x". Added
back.
> It seems that tgis.open_new_stds is actually not according to PEP8. At
> least one contains whitespace at the end, rerun pep8.
Testing another editor these days, and got absorbed by it. I work
with Spyder normally, but I am considering to move into pure vim +
plugins.
> I'm not sure if this is part of PEP8 but "comment graphics" is often
> discouraged; I personally don't see a reason for:
>
> # add timestamps either via temporal framework ----------------------
> # ------------------------------------------ outputs timestamped ---
>
> It would be much better to include there actual sentence(s) which will be
> explanatory even for everybody. (I know what is behind "outputs
> timestamped" because I remember Soeren's commit but I don't think that this
> is well documented and generally known.)
>
> Update Copyright year (I think it is OK to use 2013-2015 for simplicity -
> precise info is in Subversion).
>
> When you are doing stylistic changes, you can also put spaces around
> assignment in r.mapcalc expressions.
You mean 'output = expression' instead of 'output=expression'?
> For your todo: picture to the manual would be nice (perhaps two or three
> from a series).
Thnx :-), Nikos
-------------- next part --------------
Index: r.sun.daily.py
===================================================================
--- r.sun.daily.py (revision 64893)
+++ r.sun.daily.py (working copy)
@@ -2,11 +2,12 @@
############################################################################
#
-# MODULE: r.sun.daily for GRASS 6 and 7; based on rsun_crop.sh from GRASS book
+# MODULE: r.sun.daily; based on rsun_crop.sh from GRASS book
# AUTHOR(S): Vaclav Petras, Anna Petrasova
+# Nikos Alexandris (updates for linke, albedo, latitude, horizon)
#
# PURPOSE:
-# COPYRIGHT: (C) 2013 by the GRASS Development Team
+# COPYRIGHT: (C) 2013 - 2015 by the GRASS Development Team
#
# This program is free software under the GNU General Public
# License (>=v2). Read the file COPYING that comes with GRASS
@@ -16,29 +17,115 @@
#%module
#% description: Runs r.sun for multiple days in loop (mode 2)
-#% keyword: raster
-#% keyword: sun
+#% keywords: raster
+#% keywords: sun
#%end
-#%option
+
+#%option G_OPT_R_ELEV
+#% key: elevation
#% type: string
+#% description: Name of the input elevation raster map [meters]
#% gisprompt: old,cell,raster
-#% key: elevation
-#% description: Name of the input elevation raster map [meters]
#% required : yes
#%end
+
#%option
+#% key: aspect
#% type: string
#% gisprompt: old,cell,raster
-#% key: aspect
#% description: Name of the input aspect map (terrain aspect or azimuth of the solar panel) [decimal degrees]
#%end
+
#%option
+#% key: slope
#% type: string
#% gisprompt: old,cell,raster
-#% key: slope
#% description: Name of the input slope raster map (terrain slope or solar panel inclination) [decimal degrees]
#%end
+
+#%option G_OPT_R_INPUT
+#% key: linke
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Name of the Linke atmospheric turbidity coefficient input raster map [-]
+#% required : no
+#%end
+
#%option
+#% key: linke_value
+#% key_desc: float
+#% type: double
+#% description: A single value of the Linke atmospheric turbidity coefficient [-]
+#% options: 0.0-7.0
+#% answer: 3.0
+#% required : no
+#%end
+
+#% rules
+#% exclusive: linke, linke_value
+#% end
+
+#%option G_OPT_R_INPUT
+#% key: albedo
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Name of the Linke atmospheric turbidity coefficient input raster map [-]
+#% required : no
+#%end
+
+#%option
+#% key: albedo_value
+#% key_desc: float
+#% type: double
+#% description: A single value of the ground albedo coefficient [-]
+#% options: 0.0-1.0
+#% answer: 0.2
+#% required : no
+#%end
+
+#% rules
+#% exclusive: albedo, albedo_value
+#% end
+
+#%option G_OPT_R_INPUT
+#% key: lat
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Name of input raster map containing latitudes (if projection undefined) [decimal degrees]
+#% required : no
+#%end
+
+#%option G_OPT_R_INPUT
+#% key: long
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Name of input raster map containing longitudes (if projection undefined) [decimal degrees]
+#% required : no
+#%end
+
+#%option G_OPT_R_BASENAME_INPUT
+#% key: horizon_basename
+#% key_desc: basename
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: The horizon information input map basename
+#% required : no
+#%end
+
+#%option
+#% key: horizon_step
+#% key_desc: stepsize
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Angle step size for multidirectional horizon [degrees]
+#% required : no
+#%end
+
+#% rules
+#% requires_all: horizon_basename, horizon_step
+#% end
+
+#%option
#% key: start_day
#% type: integer
#% description: Start day (of year) of interval
@@ -45,6 +132,7 @@
#% options: 1-365
#% required : yes
#%end
+
#%option
#% key: end_day
#% type: integer
@@ -52,6 +140,7 @@
#% options: 1-365
#% required : yes
#%end
+
#%option
#% key: day_step
#% type: integer
@@ -59,6 +148,7 @@
#% options: 1-365
#% answer: 1
#%end
+
#%option
#% key: step
#% type: double
@@ -65,39 +155,45 @@
#% description: Time step when computing all-day radiation sums [decimal hours]
#% answer: 0.5
#%end
+
#%option
#% key: civil_time
#% type: double
#% description: Civil time zone value, if none, the time will be local solar time
#%end
+
#%option
+#% key: beam_rad
#% type: string
#% gisprompt: new,cell,raster
-#% key: beam_rad
#% description: Output beam irradiation raster map cumulated for the whole period of time [Wh.m-2.day-1]
#% required: no
#%end
+
#%option
+#% key: diff_rad
#% type: string
#% gisprompt: new,cell,raster
-#% key: diff_rad
#% description: Output diffuse irradiation raster map cumulated for the whole period of time [Wh.m-2.day-1]
#% required: no
#%end
+
#%option
+#% key: refl_rad
#% type: string
#% gisprompt: new,cell,raster
-#% key: refl_rad
#% description: Output ground reflected irradiation raster map cumulated for the whole period of time [Wh.m-2.day-1]
#% required: no
#%end
+
#%option
+#% key: glob_rad
#% type: string
#% gisprompt: new,cell,raster
-#% key: glob_rad
#% description: Output global (total) irradiance/irradiation raster map cumulated for the whole period of time [Wh.m-2.day-1]
#% required: no
#%end
+
#%option
#% key: beam_rad_basename
#% type: string
@@ -104,6 +200,7 @@
#% label: Base name for output beam irradiation raster maps [Wh.m-2.day-1]
#% description: Underscore and day number are added to the base name for daily maps
#%end
+
#%option
#% key: diff_rad_basename
#% type: string
@@ -110,6 +207,7 @@
#% label: Base name for output diffuse irradiation raster maps [Wh.m-2.day-1]
#% description: Underscore and day number are added to the base name for daily maps
#%end
+
#%option
#% key: refl_rad_basename
#% type: string
@@ -116,6 +214,7 @@
#% label: Base name for output ground reflected irradiation raster maps [Wh.m-2.day-1]
#% description: Underscore and day number are added to the base name for daily maps
#%end
+
#%option
#% key: glob_rad_basename
#% type: string
@@ -122,6 +221,7 @@
#% label: Base name for output global (total) irradiance/irradiation raster maps [Wh.m-2.day-1]
#% description: Underscore and day number are added to the base name for daily maps
#%end
+
#%option
#% key: nprocs
#% type: integer
@@ -129,6 +229,7 @@
#% options: 1-
#% answer: 1
#%end
+
#%flag
#% key: t
#% description: Dataset name is the same as the base name for the output series of maps
@@ -149,6 +250,9 @@
def cleanup():
+ """
+ Clean up temporary maps
+ """
if REMOVE or MREMOVE:
core.info(_("Cleaning temporary maps..."))
for rast in REMOVE:
@@ -155,25 +259,36 @@
grass.run_command('g.remove', type='raster', name=rast, flags='f',
quiet=True)
for pattern in MREMOVE:
- grass.run_command('g.remove', type='raster', pattern='%s*' % pattern,
+ grass.run_command('g.remove', type='raster',
+ pattern='{pattern}*'.format(pattern=pattern),
flags='f', quiet=True)
-def is_grass_7():
- if core.version()['version'].split('.')[0] == '7':
- return True
- return False
-
-
def create_tmp_map_name(name):
+ """
+ Create temporary map names
+ """
return '{mod}{pid}_{map_}_tmp'.format(mod='r_sun_crop',
pid=os.getpid(),
map_=name)
-# add latitude map
-def run_r_sun(elevation, aspect, slope, day, step,
- beam_rad, diff_rad, refl_rad, glob_rad, suffix):
+def run_r_sun(elevation, aspect, slope, latitude, longitude,
+ linke, linke_value, albedo, albedo_value,
+ horizon_basename, horizon_step,
+ day, step, beam_rad, diff_rad, refl_rad, glob_rad, suffix):
+ '''
+ Execute r.sun using the provided input options. Except for the required
+ parameters, the function updates the list of optional/selected parameters
+ to instruct for user requested inputs and outputs.
+ Optional inputs:
+
+ - latitude
+ - longitude
+ - linke OR single linke value
+ - albedo OR single albedo value
+ - horizon maps
+ '''
params = {}
if beam_rad:
params.update({'beam_rad': beam_rad + suffix})
@@ -183,50 +298,66 @@
params.update({'refl_rad': refl_rad + suffix})
if glob_rad:
params.update({'glob_rad': glob_rad + suffix})
+ if linke:
+ params.update({'linke': linke})
+ if linke_value:
+ params.update({'linke_value': linke_value})
+ if albedo:
+ params.update({'albedo': albedo})
+ if albedo_value:
+ params.update({'albedo_value': albedo_value})
+ if horizon_basename and horizon_step:
+ params.update({'horizon_basename': horizon_basename})
+ params.update({'horizon_step': horizon_step})
- if is_grass_7():
- grass.run_command('r.sun', elevation=elevation, aspect=aspect,
- slope=slope, day=day, step=step,
- overwrite=core.overwrite(), quiet=True,
- **params)
- else:
- grass.run_command('r.sun', elevin=elevation, aspin=aspect,
- slopein=slope,
- day=day, step=step,
- overwrite=core.overwrite(), quiet=True,
- **params)
+ grass.run_command('r.sun', elevation=elevation, aspect=aspect,
+ slope=slope, day=day, step=step,
+ overwrite=core.overwrite(), quiet=True,
+ **params)
def set_color_table(rasters):
- if is_grass_7():
- grass.run_command('r.colors', map=rasters, col='gyr', quiet=True)
- else:
- for rast in rasters:
- grass.run_command('r.colors', map=rast, col='gyr', quiet=True)
+ """
+ Set 'gyr' color tables for raster maps
+ """
+ grass.run_command('r.colors', map=rasters, col='gyr', quiet=True)
def set_time_stamp(raster, day):
+ """
+ Timestamp script's daily output map
+ """
grass.run_command('r.timestamp', map=raster, date='%d days' % day,
quiet=True)
def format_order(number, zeros=3):
+ """
+ Add leading zeros to string
+ """
return str(number).zfill(zeros)
def check_daily_map_names(basename, mapset, start_day, end_day, day_step):
+ """
+ Check if maps exist with name(s) identical to the scripts intented outputs
+ """
if not basename:
return
for day in range(start_day, end_day + 1, day_step):
- map_ = '%s%s%s' % (basename, '_', format_order(day))
+ map_ = '{base}_{day}'.format(base=basename, day=format_order(day))
if grass.find_file(map_, element='cell', mapset=mapset)['file']:
- grass.fatal(_("Raster map <%s> already exists. Change the base "
- "name or allow overwrite.") % map_)
+ grass.fatal(_('Raster map <{name}> already exists. '
+ 'Change the base name or allow overwrite.'
+ ''.format(name=map_)))
def sum_maps(sum_, basename, suffixes):
+ """
+ Sum up multiple raster maps
+ """
maps = '+'.join([basename + suf for suf in suffixes])
- grass.mapcalc('{sum_}={sum_} + {new}'.format(sum_=sum_, new=maps),
+ grass.mapcalc('{sum_} = {sum_} + {new}'.format(sum_=sum_, new=maps),
overwrite=True, quiet=True)
@@ -233,20 +364,32 @@
def main():
options, flags = grass.parser()
+ # required
elevation_input = options['elevation']
aspect_input = options['aspect']
slope_input = options['slope']
+ # optional
+ latitude = options['lat']
+ longitude = options['long']
+ linke_input = options['linke']
+ linke_value = options['linke_value']
+ albedo_input = options['albedo']
+ albedo_value = options['albedo_value']
+ horizon_basename = options['horizon_basename']
+ horizon_step = options['horizon_step']
+
+ # outputs
beam_rad = options['beam_rad']
diff_rad = options['diff_rad']
refl_rad = options['refl_rad']
glob_rad = options['glob_rad']
-
beam_rad_basename = beam_rad_basename_user = options['beam_rad_basename']
diff_rad_basename = diff_rad_basename_user = options['diff_rad_basename']
refl_rad_basename = refl_rad_basename_user = options['refl_rad_basename']
glob_rad_basename = glob_rad_basename_user = options['glob_rad_basename']
+ # missing output?
if not any([beam_rad, diff_rad, refl_rad, glob_rad,
beam_rad_basename, diff_rad_basename,
refl_rad_basename, glob_rad_basename]):
@@ -270,10 +413,6 @@
nprocs = int(options['nprocs'])
- temporal = flags['t']
- if not is_grass_7() and temporal:
- grass.warning(_("Flag t has effect only in GRASS 7"))
-
if beam_rad and not beam_rad_basename:
beam_rad_basename = create_tmp_map_name('beam_rad')
MREMOVE.append(beam_rad_basename)
@@ -287,7 +426,7 @@
glob_rad_basename = create_tmp_map_name('glob_rad')
MREMOVE.append(glob_rad_basename)
- # here we check all the days
+ # check for existing identical map names
if not grass.overwrite():
check_daily_map_names(beam_rad_basename, grass.gisenv()['MAPSET'],
start_day, end_day, day_step)
@@ -315,13 +454,13 @@
quiet=True, **params)
if beam_rad:
- grass.mapcalc('%s=0' % beam_rad, quiet=True)
+ grass.mapcalc('{beam} = 0'.format(beam=beam_rad), quiet=True)
if diff_rad:
- grass.mapcalc('%s=0' % diff_rad, quiet=True)
+ grass.mapcalc('{diff} = 0'.format(diff=diff_rad), quiet=True)
if refl_rad:
- grass.mapcalc('%s=0' % refl_rad, quiet=True)
+ grass.mapcalc('{refl} = 0'.format(refl=refl_rad), quiet=True)
if glob_rad:
- grass.mapcalc('%s=0' % glob_rad, quiet=True)
+ grass.mapcalc('{glob} = 0'.format(glob=glob_rad), quiet=True)
grass.info(_("Running r.sun in a loop..."))
count = 0
@@ -339,8 +478,13 @@
suffix = '_' + format_order(day)
proc_list.append(Process(target=run_r_sun,
- args=(elevation_input, aspect_input,
- slope_input, day, step,
+ args=(elevation_input,
+ aspect_input, slope_input,
+ latitude, longitude,
+ linke_input, linke_value,
+ albedo_input, albedo_value,
+ horizon_basename, horizon_step,
+ day, step,
beam_rad_basename,
diff_rad_basename,
refl_rad_basename,
@@ -374,6 +518,7 @@
# Empty process list
proc_list = []
suffixes = []
+
# FIXME: how percent really works?
# core.percent(1, 1, 1)
@@ -391,22 +536,29 @@
refl_rad_basename_user, glob_rad_basename_user]):
return 0
- # add timestamps either via temporal framework in 7 or r.timestamp in 6.x
- if is_grass_7() and temporal:
+ # add timestamps and register to spatio-temporal raster data set
+ temporal = flags['t']
+ if temporal:
core.info(_("Registering created maps into temporal dataset..."))
import grass.temporal as tgis
def registerToTemporal(basename, suffixes, mapset, start_day, day_step,
title, desc):
+ """
+ Register daily output maps in spatio-temporal raster data set
+ """
maps = ','.join([basename + suf + '@' + mapset for suf in suffixes])
- tgis.open_new_stds(basename, type='strds',
- temporaltype='relative',
- title=title, descr=desc,
- semantic='sum', dbif=None,
- overwrite=grass.overwrite())
- tgis.register_maps_in_space_time_dataset(
- type='raster', name=basename, maps=maps, start=start_day, end=None,
- unit='days', increment=day_step, dbif=None, interval=False)
+ tgis.open_new_stds(basename, type='strds', temporaltype='relative',
+ title=title, descr=desc, semantic='sum',
+ dbif=None, overwrite=grass.overwrite())
+
+ tgis.register_maps_in_space_time_dataset(type='rast',
+ name=basename, maps=maps,
+ start=start_day, end=None,
+ unit='days',
+ increment=day_step,
+ dbif=None, interval=False)
+
# Make sure the temporal database exists
tgis.init()
@@ -430,6 +582,7 @@
start_day, day_step, title="Total irradiation",
desc="Output total irradiation raster maps [Wh.m-2.day-1]")
+ # just add timestamps, don't register
else:
for i, day in enumerate(days):
if beam_rad_basename_user:
@@ -441,6 +594,7 @@
if glob_rad_basename_user:
set_time_stamp(glob_rad_basename + suffixes_all[i], day=day)
+ # set color table for daily maps
if beam_rad_basename_user:
maps = [beam_rad_basename + suf for suf in suffixes_all]
set_color_table(maps)
More information about the grass-dev
mailing list