radiation: Estimation of Direct and Diffuse Radiation Load on Extruded...
In shadow: Geometric Shadow Calculations

Description Usage Arguments Details Value Examples

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")
)

`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

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

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

# 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
)
rad1

## Not run: 

# 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


## End(Not run)

## Not run: 

### 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)

## End(Not run)