* Vaclav Petras <wenzesl...@gmail.com> [2015-03-26 23:12:50 -0400]: > On Thu, Mar 26, 2015 at 5:04 AM, Nikos Alexandris <n...@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
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)
_______________________________________________ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev