Estimation of Direct and Diffuse Radiation Load on Extruded Polygon Surfaces

This is a wrapper function for calculating total diffuse, direct and total radiation load per unit area on extruded polygon surfaces. The function operates on obstacle geometry and a set of sun positions with associated meteorological estimates for direct and diffuse radiation (see Details below).

radiation(
  grid,
  obstacles,
  obstacles_height_field,
  solar_pos = solarpos2(obstacles, time),
  time = NULL,
  solar_normal,
  solar_diffuse,
  radius = Inf,
  returnList = FALSE,
  parallel = getOption("mc.cores")
)

Arguments

grid	A 3D `SpatialPointsDataFrame` layer, such as returned by function `surfaceGrid`, specifying the locations where radiation is to be estimated. The layer must include an attribute named `type`, with possible values being `"roof"` or `"facade"`, expressing surface orientation per 3D point. The layer must also include an attribute named `facade_az`, specifying facade azimuth (only for `"facade"` points, for `"roof"` points the value should be `NA`). The `type` and `facade_az` attributes are automatically created when creating the grid with the `surfaceGrid` function
obstacles	A `SpatialPolygonsDataFrame` object specifying the obstacles outline, inducing self- and mutual-shading on the grid points
obstacles_height_field	Name of attribute in `obstacles` with extrusion height for each feature
solar_pos	A `matrix` with two columns representing sun position(s); first column is the solar azimuth (in decimal degrees from North), second column is sun elevation (in decimal degrees); rows represent different sun positions corresponding to the `solar_normal` and the `solar_diffuse` estimates. For example, if `solar_normal` and `solar_diffuse` refer to hourly measurements in a Typical Meteorological Year (TMY) dataset, then `solar_pos` needs to contain the corresponding hourly sun positions
time	When `solar_pos` is unspecified, `time` can be passed to automatically calculate `solar_pos` based on the time and the centroid of `obstacles`, using function `maptools::solarpos`. In such case `obstacles` must have a defined CRS (not `NA`). The `time` value must be a `POSIXct` or `POSIXlt` object
solar_normal	Direct Normal Irradiance (e.g. in Wh/m^2), at sun positions corresponding to `solar_pos`. Must be a vector with the same number of elements as the number of rows in `solar_pos`
solar_diffuse	Diffuse Horizontal Irradiance (e.g. in Wh/m^2), at sun positions corresponding to `solar_pos`. Must be a vector with the same number of elements as the number of rows in `solar_pos`
radius	Effective search radius (in CRS units) for considering obstacles when calculating shadow and SVF. The default is to use a global search, i.e. `radius=Inf`. Using a smaller radius can be used to speed up the computation when working on large areas. Note that the search radius is not specific per grid point; instead, a buffer is applied on all grid points combined, then "dissolving" the individual buffers, so that exactly the same obstacles apply to all grid points
returnList	Logical, determines whether to return summed radiation over the entire period per 3D point (default, `FALSE`), or to return a list with all radiation values per time step (`TRUE`)
parallel	Number of parallel processes or a predefined socket cluster. With `parallel=1` uses ordinary, non-parallel processing. Parallel processing is done with the `parallel` package

Value

If returnList=FALSE (the default), then returned object is a data.frame, with rows corresponding to grid points and four columns corresponding to the following estimates:

svf Computed Sky View Factor (see function SVF)
direct Total direct radiation for each grid point
diffuse Total diffuse radiation for each grid point
total Total radiation (direct + diffuse) for each grid point

Each row of the data.frame gives summed radiation values for the entire time period in solar_pos, solar_normal and solar_diffuse

If returnList=TRUE then returned object is a list with two elements:

direct Total direct radiation for each grid point
diffuse Total diffuse radiation for each grid point

Each of the elements is a matrix with rows corresponding to grid points and columns corresponding to time steps in solar_pos, solar_normal and solar_diffuse

Details

Input arguments for this function comprise the following:

An extruded polygon obstacles layer (obstacles and obstacles_height_field) inducing shading on the queried grid
A grid of 3D points (grid) where radiation is to be estimated. May be created from the 'obstacles' layer, or a subset of it, using function surfaceGrid. For instance, in the code example (see below) radiation is estimated on a grid covering just one of four buildings in the build layer (the first building), but all four buildings are taken into account for evaluating self- and mutual-shading by the buildings.
Solar positions matrix (solar_pos)
Direct and diffuse radiation meteorological estimate vectors (solar_normal and solar_diffuse)

Given these inputs, the function goes through the following steps:

Determining whether each grid point is shaded, at each solar position, using inShadow
Calculating the coefficient of Direct Normal Irradiance reduction, using coefDirect
Summing direct radiation considering (1) mutual shading, (2) direst radiation coefficient and (3) direct radiation estimates
Calculating the Sky View Factor (SVF) for each point, using SVF
Summing diffuse radiation load considering (1) SVF and (2) diffuse radiation estimates
Summing total (direct + diffuse) radiation load

Examples


# Create surface grid
grid = surfaceGrid(
  obstacles = build[1, ],
  obstacles_height_field = "BLDG_HT",
  res = 2
)

solar_pos = tmy[, c("sun_az", "sun_elev")]
solar_pos = as.matrix(solar_pos)

# Summed 10-hour radiation estimates for two 3D points
rad1 = radiation(
  grid = grid[1:2, ],
  obstacles = build,
  obstacles_height_field = "BLDG_HT",
  solar_pos = solar_pos[8:17, , drop = FALSE],
  solar_normal = tmy$solar_normal[8:17],
  solar_diffuse = tmy$solar_diffuse[8:17],
  returnList = TRUE
)
#> Determining Shadow
#> 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=======                                                               |  10%
  |                                                                            
  |==============                                                        |  20%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |============================                                          |  40%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |==========================================                            |  60%
  |                                                                            
  |=================================================                     |  70%
  |                                                                            
  |========================================================              |  80%
  |                                                                            
  |===============================================================       |  90%
  |                                                                            
  |======================================================================| 100%
#> Done
#> Calculating SVF
#> Done
rad1
#> $direct
#>          [,1]     [,2]     [,3]   [,4]     [,5]     [,6]     [,7]     [,8]
#> [1,] 131.5884 293.0234 399.9672 475.75 490.5326 437.2193 319.5239 175.1429
#> [2,] 131.5884 293.0234 399.9672 475.75 490.5326 437.2193 319.5239 175.1429
#>          [,9] [,10]
#> [1,] 32.61533     0
#> [2,] 32.61533     0
#> 
#> $diffuse
#>          [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]
#> [1,] 36.29564 45.97448 53.23361 53.23361 54.04018 55.65332 54.84675 45.16791
#> [2,] 44.38150 56.21657 65.09287 65.09287 66.07913 68.05164 67.06538 55.23032
#>          [,9]     [,10]
#> [1,] 17.74454 0.8065699
#> [2,] 21.69762 0.9862556
#> 

if (FALSE) {

# Same, using 'time' instead of 'solar_pos'

rad2 = radiation(
  grid = grid[1:2, ],
  obstacles = build,
  obstacles_height_field = "BLDG_HT",
  time = as.POSIXct(tmy$time[8:17], tz = "Asia/Jerusalem"),
  solar_normal = tmy$solar_normal[8:17],
  solar_diffuse = tmy$solar_diffuse[8:17],
  returnList = TRUE
)
rad2

# Differences due to the fact that 'tmy' data come with their own
# solar positions, not exactly matching those calulated using 'maptools::solarpos'
rad1$direct - rad2$direct
rad1$diffuse - rad2$diffuse

}

if (FALSE) {

### Warning! The calculation below takes some time.

# Annual radiation estimates for entire surface of one building
rad = radiation(
  grid = grid,
  obstacles = build,
  obstacles_height_field = "BLDG_HT",
  solar_pos = solar_pos,
  solar_normal = tmy$solar_normal,
  solar_diffuse = tmy$solar_diffuse,
  parallel = 3
)

# 3D plot of the results
library(plot3D)
opar = par(mfrow=c(1, 3))

scatter3D(
  x = coordinates(grid)[, 1],
  y = coordinates(grid)[, 2],
  z = coordinates(grid)[, 3],
  colvar = rad$direct / 1000,
  scale = FALSE,
  theta = 55,
  pch = 20,
  cex = 1.35,
  clab = expression(paste("kWh / ", m^2)),
  main = "Direct radiation"
)
scatter3D(
  x = coordinates(grid)[, 1],
  y = coordinates(grid)[, 2],
  z = coordinates(grid)[, 3],
  colvar = rad$diffuse / 1000,
  scale = FALSE,
  theta = 55,
  pch = 20,
  cex = 1.35,
  clab = expression(paste("kWh / ", m^2)),
  main = "Diffuse radiation"
)
scatter3D(
  x = coordinates(grid)[, 1],
  y = coordinates(grid)[, 2],
  z = coordinates(grid)[, 3],
  colvar = rad$total / 1000,
  scale = FALSE,
  theta = 55,
  pch = 20,
  cex = 1.35,
  clab = expression(paste("kWh / ", m^2)),
  main = "Total radiation"
)

par(opar)
}