For latitude-longitude coordinates it requires that the elevation map is in meters. The rules are:
The model considers a shadowing effect of the local topography unless switched off with the -p flag. r.sun works in two modes: In the first mode it calculates for the set local time a solar incidence angle [degrees] and solar irradiance values [W.m-2]. In the second mode daily sums of solar radiation [Wh.m-2.day-1] are computed within a set day. By a scripting the two modes can be used separately or in a combination to provide estimates for any desired time interval. The model accounts for sky obstruction by local relief features. Several solar parameters are saved in the resultant maps' history files, which may be viewed with the r.info command.
The solar incidence angle raster map incidout is computed specifying
elevation raster map elevation, aspect raster map aspect, slope
steepness raster map slope, given the day day and local time
time. There is no need to define latitude for locations with known
and defined projection/coordinate system (check it with the
command). If you have undefined projection, (x,y) system, etc. then the latitude
can be defined explicitly for large areas by input raster map lat_in
with interpolated latitude values. All input raster maps must
be floating point (FCELL) raster maps. Null data in maps are excluded from
the computation (and also speeding-up the computation), so each output raster
map will contain null data in cells according to all input raster maps. The
user can use r.null
command to create/reset null file for your input raster maps.
The specified day day is the number of the day of the general year where January 1 is day no.1 and December 31 is 365. Time time must be a local (solar) time (i.e. NOT a zone time, e.g. GMT, CET) in decimal system, e.g. 7.5 (= 7h 30m A.M.), 16.1 = 4h 6m P.M..
The solar declination parameter is an option to override the value computed by the internal routine for the day of the year. The value of geographical latitude can be set as a constant for the whole computed region or, as an option, a grid representing spatially distributed values over a large region. The geographical latitude must be also in decimal system with positive values for northern hemisphere and negative for southern one. In similar principle the Linke turbidity factor (linke, lin ) and ground albedo (albedo, alb) can be set.
Besides clear-sky radiations, the user can compute a real-sky radiation (beam, diffuse) using coeff_bh and coeff_dh input raster maps defining the fraction of the respective clear-sky radiations reduced by atmospheric factors (e.g. cloudiness). The value is between 0-1. Usually these coefficients can be obtained from a long-terms meteorological measurements provided as raster maps with spatial distribution of these coefficients separately for beam and diffuse radiation (see Suri and Hofierka, 2004, section 3.2).
The solar irradiation or irradiance raster maps beam_rad, diff_rad, refl_rad are computed for a given day day, latitude lat_in, elevation elevation, slope slope and aspect aspect raster maps. For convenience, the output raster given as glob_rad will output the sum of the three radiation components. The program uses the Linke atmosphere turbidity factor and ground albedo coefficient. A default, single value of Linke factor is lin=3.0 and is near the annual average for rural-city areas. The Linke factor for an absolutely clear atmosphere is lin=1.0. See notes below to learn more about this factor. The incidence solar angle is the angle between horizon and solar beam vector.
The solar radiation maps for a given day are computed by integrating the relevant irradiance between sunrise and sunset times for that day. The user can set a finer or coarser time step used for all-day radiation calculations with the step option. The default value of step is 0.5 hour. Larger steps (e.g. 1.0-2.0) can speed-up calculations but produce less reliable (and more jagged) results. As the sun moves through approx. 15° of the sky in an hour, the default step of half an hour will produce 7.5° steps in the data. For relatively smooth output with the sun placed for every degree of movement in the sky you should set the step to 4 minutes or less. step=0.05 is equivalent to every 3 minutes. Of course setting the time step to be very fine proportionally increases the module's running time.
The output units are in Wh per squared meter per given day [Wh/(m*m)/day]. The incidence angle and irradiance/irradiation maps are computed with the shadowing influence of relief by default. It is also possible for them to be computed without this influence using the planar flag (-p). In mountainous areas this can lead to very different results! The user should be aware that taking into account the shadowing effect of relief can slow down the speed of computation, especially when the sun altitude is low.
When considering the shadowing effect, speed and precision of computation can be controlled by the distance_step parameter, which defines the sampling density at which the visibility of a grid cell is computed in the direction of a path of the solar flow. It also defines the method by which the obstacle's altitude is computed. When choosing a distance_step less than 1.0 (i.e. sampling points will be computed at distance_step * cellsize distance), r.sun takes the altitude from the nearest grid point. Values above 1.0 will use the maximum altitude value found in the nearest 4 surrounding grid points. The default value distance_step=1.0 should give reasonable results for most cases (e.g. on DEM). The distance_step value defines a multiplying coefficient for sampling distance. This basic sampling distance equals to the arithmetic average of both cell sizes. The reasonable values are in the range 0.5-1.5. The values below 0.5 will decrease and values above 1.0 will increase the computing speed. Values greater than 2.0 may produce estimates with lower accuracy in highly dissected relief. The fully shadowed areas are written to the output maps as zero values. Areas with NULL data are considered as no barrier with shadowing effect.
The maps' history files are generated containing the following listed
parameters used in the computation:
- Solar constant used W.m-2
- Extraterrestrial irradiance on a plane perpendicular to the solar beam [W.m-2]
- Day of the year
- Declination [radians]
- Decimal hour (Alternative 1 only)
- Sunrise and sunset (min-max) over a horizontal plane
- Daylight lengths
- Geographical latitude (min-max)
- Linke turbidity factor (min-max)
- Ground albedo (min-max)
The user can use a nice shellcript with variable day to compute radiation for some time interval within the year (e.g. vegetation or winter period). Elevation, aspect and slope input values should not be reclassified into coarser categories. This could lead to incorrect results.
Currently, there are two modes of r.sun. In the first mode it calculates solar incidence angle and solar irradiance raster maps using the set local time. In the second mode daily sums of solar irradiation [Wh.m-2.day-1] are computed for a specified day.
The clear-sky solar radiation model applied in the r.sun is based on the work undertaken for development of European Solar Radiation Atlas (Scharmer and Greif 2000, Page et al. 2001, Rigollier 2001). The clear sky model estimates the global radiation from the sum of its beam, diffuse and reflected components. The main difference between solar radiation models for inclined surfaces in Europe is the treatment of the diffuse component. In the European climate this component is often the largest source of estimation error. Taking into consideration the existing models and their limitation the European Solar Radiation Atlas team selected the Muneer (1990) model as it has a sound theoretical basis and thus more potential for later improvement.
Details of underlying equations used in this program can be found in the reference literature cited below or book published by Neteler and Mitasova: Open Source GIS: A GRASS GIS Approach (published in Kluwer Academic Publishers in 2002).
Average monthly values of the Linke turbidity coefficient for a mild climate in the northern hemisphere (see reference literature for your study area):
Planned improvements include the use of the SOLPOS algorithm for solar geometry calculations and internal computation of aspect and slope.
To overcome this problem, the user can use the option civil_time=<timezone_offset> in r.sun to make it use real-world (wall clock) time. For example, for Central Europe the timezone offset is +1, +2 when daylight saving time is in effect.
# without input raster map partitioning: # memory requirements: 4 bytes per raster cell # rows,cols: rows and columns of current region (find out with g.region) # IR: number of input raster maps without horizon maps # OR: number of output raster maps memory_bytes = rows*cols*(IR*4 + horizon_steps + OR*4) # with input raster map partitioning: memory_bytes = rows*cols*((IR*4+horizon_steps)/npartitions + OR*4)
g.region raster=elevation -p # calculate horizon angles (to speed up the subsequent r.sun calculation) r.horizon elevation=elevation step=30 bufferzone=200 output=horangle \ maxdistance=5000 # slope + aspect r.slope.aspect elevation=elevation aspect=aspect.dem slope=slope.dem # calculate global radiation for day 180 at 2p.m., using r.horizon output r.sun elevation=elevation horizon_basename=horangle horizon_step=30 \ aspect=aspect.dem slope=slope.dem glob_rad=global_rad day=180 time=14 # result: output global (total) irradiance/irradiation [W.m-2] for given day/time r.univar global_rad
Calculation of the integrated daily irradiation for a region in North-Carolina for a given day of the year at 30m resolution. Here day 172 (i.e., 21 June in non-leap years):
g.region raster=elev_ned_30m -p # considering cast shadows r.sun elevation=elev_ned_30m linke_value=2.5 albedo_value=0.2 day=172 \ beam_rad=b172 diff_rad=d172 \ refl_rad=r172 insol_time=it172 d.mon wx0 # show irradiation raster map [Wh.m-2.day-1] d.rast.leg b172 # show insolation time raster map [h] d.rast.leg it172
>>> import datetime >>> datetime.datetime(2014, 6, 21).timetuple().tm_yday 172
Latest change: Thu Feb 3 11:10:06 2022 in commit: 73413160a81ed43e7a5ca0dc16f0b56e450e9fef
© 2003-2022 GRASS Development Team, GRASS GIS 8.0.3dev Reference Manual