Clear sky Radiation on a digital elevation model¶

The radiation module can be used to create clear sky radiation estimates using a (high resolution) digital elevation model. In turn, these results can be used by the evaporation module to better estimate local evaporation.

Usage¶

Below the command-line usage of the e2o_radiation.py script is described:

e2o_radiation -D DEM [-d DEM][-O outputdir][-S start day][-E end day]
          [-M][-x lon][-y lat][-h][-l loglevel][-T minutes][-L directory]
          [-t transmissivity]

-D DEM Filename of the digital elevation model
-d DEM option low resolution dem to determine height correction for
-O outputdir (default is . )
-S Startday - Start day of the simulation (1 Jan is 1)
-E EndDay - End day of the simulation
-T minutes - timeresolution in minutes (60 default is 1 hour)
-M The DEM xy units are in metres (instead of lat/lon)
-x longitute of the map left (if map xy in metres)
-y lattitude of the map bottom (if map xy in metres)
-l loglevel Set loglevel (DEBUG, INFO, WARNING,ERROR)
-s start hour (per day) of the calculations (default =1)
-e end hour (per day) of the calculations (default = 23)
-h This information
-f output format as a gdal type string (http://www.gdal.org/formats_list.html) Default is PCRaster
-p add the file format string as a postfix to all filename (default is False)
-t specify thau transmissivity (default value = 0.6) (usually between 0.4 and 0.8)
-L determine optical correction via linke turbidity (in stead of thau transmissivity). The
   argument should point to a directory in which monthly linke maps are present in pcraster format

The following command takes the local digital elevation model (wflow_dem.map) and calculates the radiation components for day 1 to 366 and saves these as a pcraster mapstack. The internal timesteps is set to 60 minutes and the calculations are done for 5 hr in the morning to 22 hour at night.

e2o_radiation.py -S 1 -E 366 -D wflow_dem.map -O output_rad -s 5 -e 22 -T 60 -f PCRaster

Methods¶

Adjusted after van Dam 2000

This section gives a short description. Another description can be found at http://re.jrc.cec.eu.int/pvgis/pv/solres/solres.htm.

Potential solar radiation is the radiation of an unobstructed or cloudless sky. The magnitude of the potential solar radiation depends on the position of the sun the solar altitude or solar angle during the day, the inclination of the solar rays with the earth’s surface, the amount of radiation at the outer layer of the earth’s atmosphere, the transmissivity of the sky and the altitude of the earth’s surface.

Solar declination is the annual fluctuation of the sun between the two tropics and varies between –23 and +23 degrees latitude. Solar declination is calculated per Day (Julian day number):

$\delta = -23.4 cos(360 (Day + 10 )/365)$

The hour angle describes the movement of the sun around the earth in 24 hours, which equals 15 degrees longitude per hour (360 $deg$ /24h). The hour angle n is calculated for each Hour (whole hour of the day):

$n = ( 15 \dot Hour - 12)$

The position or height of the sun above the horizon is called the solar altitude or solar angle. Solar altitude $\alpha$ (deg) is calculated for each location, determined by the location’s latitude $\psi$ (deg), declination and hour angle:

$sin(\alpha) = sin(\psi) sin(\delta) + cos(\psi) cos(\delta) cos(n)$

Solar azimuth is the angle between the solar rays and the North-South axis of the earth. Solar azimuth ${\beta}_s$ (deg) is calculated by:

$&cos({\beta}_s) = (sin(\delta) cos(\psi) - cos(\delta) sin(\psi) cos(n))/cos(\alpha) \\ &for Hour \le 12: {\beta}_s = {\beta}_s \\ &for Hour > 12: {\beta}_s = 360 - {\beta}_s$

Surface azimuth or aspect ${\beta}_1$ (deg) is the orientation of the land surface or slope to the North-South axis of the sun. Slope $\varphi$ (deg) is the maximum rate of change in elevation.

The angle of incidence is the angle between the perpendicular plane of the incoming solar rays and the surface on which they are projected, defined by the aspect and slope of that surface. The angle of incidence $\vartheta$ (deg) is calculated with the solar angle $\alpha$ (deg), the slope of the land surface $\varphi$ (deg), the azimuth of the sun ${\beta}_s$ (deg) and azimuth of the land surface ${\beta}_1$ (deg):

$cos(\vartheta) = cos(\alpha) sin(\varphi) cos({\beta}_s - {\beta}_1) + sin(\alpha) cos(\varphi)$

The second section of the radiation module calculates the potential solar energy. The amount of solar radiation that reaches the outer atmosphere is decreased by the travelling distance of the solar rays through the sky to the surface, the transmissivity of the sky and the cloud factor.

Solar energy at the outer layer of the atmosphere $Sout (Wm^2)$ is calculated by (Kreider & Kreith 1975):

$S_{out} = S_c (1 + 0.034 cos(360 Day/365))$

where $S_c (Wm^2)$ is the solar constant of 1367 $Wm^2$ (Duffie & Beckman 1991). The solar ‘constant’ is subject to much discussion. Gates (1980) gives a value of 1360 $Wm^2$ . The NASA reports a value of 1353 $Wm^2$ (Jansen 1985), while Duncan et al. (1982) give a value of 1367 $Wm^2$ . Monteith and Unsworth (1990) measured the highest value of 1373 W.m-2. The World Radiation Centre uses a value of 1367 $Wm^2$ (Duffie & Beckman 1991) and this value is also used in this study.

The solar radiation energy that reaches the earth’s surface is decreased due to the length of the air mass it has to pass through and the transmissivity $\tau$ (% or fraction) of the sky. The radiation flux through a hypothetical plane normal to the beam ( $S_{nor} Wm^2$ ) is given by (Gates 1980):

$S_{nor} = S_{out} \tau^{Mh}$

in which Mh (% or fraction) is the relative path length of the optical air mass at altitude h (m). Transmissivity ( $\tau$ ) is usually between 0.5 and 0.8, but can be as low as 0.4 in the tropics (Whitmore et al. 1993), but mostly a value of 0.6 is used (Gates 1980). To calculate the relative path length of an optical air mass at altitude h (m), the relative path length of an optical air mass at sea level M0 (% or fraction) is corrected for the atmospheric pressure at altitude h. Mh (% or fraction) is calculated using (Kreider & Kreith 1975):

$Mh = M_0 P_h/P_0$

in which $P_h / P_0$ (mbar.mbar-1) is an atmospheric pressure correction. The relative path length of the optical air mass at sea level M0 is obtained by (Kreider & Kreith 1975):

$M_0 = \sqrt(1299 + (614 sin(\alpha))^2) - 614 sin(\alpha)$

The atmospheric pressure correction $P_h / P_0$ is written as (List 1984):

$P_h / P_0 = ((288 0.0065h) / 288)^5.256$

The incoming radiation normal to the beam Snor must be corrected by the orientation and slope of the surface, defined by the angle of incidence $\vartheta$ , to calculate the incoming radiation Sdir ( $Wm^2$ ) on the earth’s surface:

$S_dir = S_nor cos(\vartheta)$

Direct light is scattered in the atmosphere. This daylight scattering or diffuse radiation is approximately 15% of direct radiation (Gates 1980). A more accurate empirical estimation for diffuse radiation Sdif ( $Wm^2$ ) in a clear not dust-free sky reads as (Liu and Jordan in Gates 1980):

$S_dif = S_out (0.271 - 0.294 \vartheta^{Mh} sin(\alpha)$

During daylight when the sun is above the horizon, it is assumed that all cells receive the same amount of diffuse radiation. Total incoming radiation Sin ( $Wm^2$ ) is the sum of direct and diffuse radiation:

$S_in = S_dir + S_dif$

Total incoming radiation Sin as calculated with the above is actually a radiation flux for that moment. In the procedure given above, radiation is calculated per time step. If this amount of radiation is used in a water balance model, the amount of radiation and therewith the amount of evapotranspiration will be overestimated or under estimated, depending on the time of the day and the position of the sun.

Most of the work done for the shading is implemented in the pcraster horizontan function.

Example¶

The image below show the difference between the average daily clear sky radiation on a flat surface compared to the actual radiation received by each grid cell for julian day 180. IN the steep terrain of the Snowy Mountains of the mUrumbidgee catchment differences can be up to 30 $W/m^2$ on this 1x1km DEM.

Difference in $W/m^2$ between horizontal surface radiation and inclined surface radiation

Literature¶

https://www.google.nl/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&cad=rja&uact=8&ved=0CCYQFjAA&url=http%3A%2F%2Fwww.physics.arizona.edu%2F~cronin%2FSolar%2FReferences%2FIrradiance%2520Models%2520and%2520Data%2FWOC01.pdf&ei=CHNcVfKXLMf-UojagfAN&usg=AFQjCNG4m1HBoQMPeZYcPxIQRYtC4xcp6A&sig2=L02qLYq-2vNXwgwW_vUwFA

Implementation¶

Determine clear sky radiation over a Digital Elevation Model.

Usage:

e2o_radiation -D DEM [-d DEM][-O outputdir][-S start day][-E end day]
          [-M][-x lon][-y lat][-h][-l loglevel][-T minutes][-L directory]
          [-t transmissivity]

-D DEM Filename of the digital elevation model
-d DEM option low resolution dem to determine height correction for
    use -d WRR1 or -d WRR2 for default WRR1 and WRR2 DEM's
-O outputdir (default is . )
-S Startday - Start day of the simulation (1 Jan is 1)
-E EndDay - End day of the simulation
-T minutes - timeresolution in minutes (60 default is 1 hour)
-M The DEM xy units are in metres (instead of lat/lon)
-x longitute of the map left (if map xy in metres)
-y lattitude of the map bottom (if map xy in metres)
-l loglevel Set loglevel (DEBUG, INFO, WARNING,ERROR)
-s start hour (per day) of the calculations (default =1)
-e end hour (per day) of the calculations (default = 23)
-h This information
-f output format as a gdal type string (http://www.gdal.org/formats_list.html) Default is PCRaster
-p add the file format string as a postfix to all filename (default is False)
-t specify thau transmissivity (default value = 0.6) (usually between 0.4 and 0.8)
-L determine optical correction via linke turbidity (in stead of thau transmissivity).

The program produces the following map stacks, one for each day of the year:

COR00000.??? - Total clear sky radiation on DEM
SUN00000.??? - Nr of time intervals a pixel was in the sun
FLAT0000.??? - Total clear sky radiation on a flat surface
CORDIR00.??? - Direct clear sky radiation on DEM
FLATDIR0.??? - Direct clear sky radiation on a flat surface

e2o_dstools.e2o_radiation.GenRadMaps(SaveDir, Lat, Lon, Slope, Aspect, Altitude, DegreeDem, AltDem, logje, Trans=0.6, start=1, end=2, interval=60, shour=1, ehour=23, outformat='PCRaster', Addpostfix=False, linkemapstack=None)¶

Generate radiation masp for a number of days

Parameters:

SaveDir – When to save the maps
Lat – LAttitude for each pixel
Lon – Longitude for each pixel
Slope – Slope for each pixel
Aspect – Aspect
Altitude – Altitude
DegreeDem – Altitude with elevation rescaled to degree lat/lon
AltDem – DEM to convert optcor to
logje – Logger
start – Start day of the year
end – End day of the year
interval – Interval for calculations in minutes
shour – Start hour of the calculations
outformat – File format of results as a gdal string

Returns:

Nothing

e2o_dstools.e2o_radiation.correctrad(Day, Hour, Lat, Slope, Aspect, Altitude, Altitude_UnitLatLon, AltAltitude, Trans=0.6)¶

Determines radiation over a DEM assuming clear sky for a specified hour of a day

Return StotFlat:
Variables:	Day – Day of the year (1-366) Hour – Hour of the day (0-23) Lat – map with latitudes for each grid cell Slope – Slope in degrees Aspect – Aspect in degrees relative to north for each cell Altitude – Elevation in metres Altitude_Degree – Elevation in degrees. If the actual pcraster maps are in lat lon this maps should hold the Altitude converted to degrees. If the maps are in metres this maps should also be in metres AltAltitude – DEM to transfer optcorr to (CorFac) Trans – Transmissivity tau (Gates, 1980), default = 0.6
Return Stot:	Total radiation on the dem, shadows not taken into account
Return StotCor:	Total radiation on the dem taking shadows into acount
	Total radiation on the dem assuming a flat surface
Return SUN:	Map with shade (0) or no shade (1) pixels
Return OptCOr:	Correction factor for atmosphere of Altitude DEM
Return AltOptcor:
	Correction factor for atmosphere of AltAltitude DEM

e2o_dstools.e2o_radiation.correctrad_alt(Day, Hour, Lat, Slope, Aspect, Altitude, Altitude_UnitLatLon, AltAltitude, Linke, LinkeDEM)¶

Determines radiation over a DEM assuming clear sky for a specified hour of a day. This version used a fomulation based on linke turbidity to estimate clear sky radiation, similarly to r.sun

Return StotFlat:
Variables:	Day – Day of the year (1-366) Hour – Hour of the day (0-23) Lat – map with latitudes for each grid cell Slope – Slope in degrees Aspect – Aspect in degrees relative to north for each cell Altitude – Elevation in metres Altitude_Degree – Elevation in degrees. If the actual pcraster maps are in lat lon this maps should hold the Altitude converted to degrees. If the maps are in metres this maps should also be in metres AltAltitude – DEM to transfer optcorr to (CorFac) Linke – map with linke turbidity LinkeDEM – elevation for each cell for the linke map (fro correction to actual elevation)
Return Stot:	Total radiation on the dem, shadows not taken into account
Return StotCor:	Total radiation on the dem taking shadows into acount
	Total radiation on the dem assuming a flat surface
Return SUN:	Map with shade (0) or no shade (1) pixels
Return OptCOr:	Correction factor for atmosphere of Altitude DEM
Return AltOptcor:
	Correction factor for atmosphere of AltAltitude DEM

e2o_dstools.e2o_radiation.detRealCellLength(ZeroMap, sizeinmetres)¶

Determine cellength. Always returns the length in meters.

Parameters:	ZeroMap – pcraster map object used as clone sizeinmetres – is set to one the currect cell size is assumed to be in metres, otherwise in lat,log
Returns:

e2o_dstools.e2o_radiation.lattometres(lat)¶

Determines the length of one degree lat/long at a given latitude (in meter). Code taken from http:www.nga.mil/MSISiteContent/StaticFiles/Calculators/degree.html

Parameters:	lat – map with lattitude values for each cell
Returns:	length of a cell lat, length of a cell long

e2o_dstools.e2o_radiation.main(argv=None)¶

Parameters:	argv – See usage
Returns:

e2o_dstools.e2o_radiation.usage(*args)¶

Parameters:	args –
Returns:

Table Of Contents

Previous topic

Next topic

This Page