This paper introduces the *shadow* package for R. The package provides functions for shadow-related calculations in the urban environment, namely shadow height, shadow footprint and Sky View Factor (SVF) calculations, as well as a wrapper function to estimate solar radiation while taking shadow effects into account. All functions operate on a layer of polygons with a height attribute, also known as “extruded polygons” or 2.5D vector data. Such data are associated with accuracy limitations in representing urban environments. However, unlike 3D models, polygonal layers of building outlines along with their height are abundantly available and their processing does not require specialized closed-source 3D software. The present package thus brings spatio-temporal shadow, SVF and solar radiation calculation capabilities to the open-source spatial analysis workflow in R. Package functionality is demonstrated using small reproducible examples for each function. Wider potential use cases include urban environment applications such as evaluation of micro-climatic influence for urban planning, studying urban climatic comfort and estimating photovoltaic energy production potential.

Spatial analysis of the urban environment ((Biljecki et al. 2015)) frequently requires estimating whether a given point is shaded or not, given a representation of spatial obstacles (e.g. buildings) and a time-stamp with its associated solar position. For example, we may be interested in -

- Calculating the amount of time a given roof or facade is shaded, to
determine the utility of installing photovoltaic cells for
**electricity production**(e.g. (Redweik et al. 2013)). - Calculating shadow footprint on vegetated areas, to determine the
expected influence of a tall new building on the surrounding
**microclimate**(e.g. (Bourbia and Boucheriba 2010)).

Such calculations are usually carried out using GIS-based models
((Freitas et al. 2015)), in either **vector-based 3D** or
**raster-based 2.5D** settings. Both approaches have their advantages
and limitations, as discussed in the following paragraphs.

Shadow calculations on vector-based 3D models of the urban environment
are mostly restricted to proprietary closed-source software such as
**ArcGIS** (ESRI 2017) or **SketchUp** (Google 2017), though
recently some open-source models such as **SURFSUN3D** have been
developed (Liang et al. 2015). One of the drawbacks of using closed-source
software in this context is the difficulty of adjusting the software for
specific needs and uncommon scenarios. This problem is especially acute
in research settings, where flexibility and extensibility are essential
for exploring new computational approaches. The other difficulty with
using 3D software in urban spatial analysis concerns interoperability of
file formats. Since ordinary vector spatial data formats, such as the
*ESRI Shapefile*, cannot represent three-dimensional surfaces, 3D
software is associated with specialized file formats. The latter cannot
be readily imported to a general-purpose geocomputational environment
such as R or Python (Van Rossum and Drake 2011), thus fragmenting the analysis
workflow. Moreover, most 3D software, such as those mentioned above, are
design-oriented, thus providing advanced visualization capabilities but
limited quantitative tools (calculating areas, angles, coordinates,
etc.). Finally, true-3D databases of large urban areas are difficult to
obtain, while vector-based 2.5D databases (building outline and height,
see below) are almost universal. The advantages of true-3D software are
“wasted” when the input data are 2.5D, while the disadvantages, such as
lack of quantitative procedures and data interoperability difficulties,
still remain.

Raster-based 2.5D solutions, operating on a Digital Elevation Model
(DEM) raster, are much simpler and have thus been more widely
implemented in various software for several decades
(Kumar et al. 1997; Ratti and Richens 2004). For example, raster-based
shadow calculations are available in open-source software such as the
**r.sun** command (Hofierka and Suri 2002) in *GRASS GIS*
(GRASS Development Team 2017), the **UMEP** plugin (Lindberg et al. 2018) for
*QGIS* (QGIS Development Team 2017) and package
*insol* (Corripio 2014) in R. In
the proprietary **ArcGIS** software, raster-based shadow calculations
are provided through the **Solar Analyst** extension (Fu and Rich 1999).
Thanks to this variety of tools, raster-based shadow modelling can be
easily incorporated within a general spatial analysis workflow. However,
raster-based models are more suitable for large-scale analysis of
natural terrain, rather than fine-scale urban environments, for the
following reasons -

- A raster representing surface elevation, known as a DEM, at sufficiently high resolution for the urban context, may not be available and is expensive to produce, e.g. using airborne Light Detection And Ranging (LiDAR) surveys (e.g. (Redweik et al. 2013)). Much more commonly, municipalities and other sources such as OpenStreetMap (Haklay and Weber 2008) offer 2.5D vector-based data on cities, i.e. polygonal layers of building outlines associated with height attributes.
- Rasters are composed of pixels, which have no natural association to specific urban elements, such as an individual building, thus making it more difficult to associate analysis results with the corresponding urban elements.
- Vertical surfaces, such as building facades, are rare in natural terrain yet very common in urban environments. Raster-based representation of facades is problematic since the latter correspond to (vertical) discontinuities in the 2.5D digital elevation model, requiring unintuitive workarounds (Redweik et al. 2013).

It should be noted that more specialized approaches have been recently developed to address some of the above-mentioned difficulties, but they are usually not available as software packages (e.g. (Hofierka and Zlocha 2012; Redweik et al. 2013)).

The *shadow* package
(Dorman 2019) aims at addressing these limitations by introducing a simple
2.5D vector-based algorithm for calculating shadows, Sky View Factor
(SVF) and solar radiation estimates in the urban environment. The
algorithms operate on a polygonal layer extruded to 2.5D, also known as
*Levels-of-Detail (LoD) 1* in the terminology of the CityGML standard
(Gröger and Plümer 2012). On the one hand, the advantages of individual
urban element representation (over raster-based approach) and input data
availability (over both raster-based and full 3D approaches) are
maintained. On the other hand, the drawbacks of closed-source software
and difficult interoperability (as opposed to full 3D environment) are
avoided.

As demonstrated below, functions in the *shadow* package operate on a
vector layer of obstacle outlines (e.g. buildings) along with their
heights, passed as a `"SpatialPolygonsDataFrame"`

object defined in
package *sp* (Pebesma and Bivand 2005; Bivand et al. 2013). The
latter makes incorporating shadow calculations in
*Spatial* analysis workflow
in R straightforward. Functions to calculate shadow height, shadow
ground footprint, Sky View Factor (SVF) and solar radiation are
implemented in the package.

All functions currently included in package *shadow* are based on
trigonometric relations in the triangle defined by the sun’s rays, the
ground - or a plane parallel to the ground - and an obstacle.

For example, **shadow height** at any given ground point can be
calculated based on (1) sun elevation, (2) the height of the building(s)
that stand in the way of sun rays and (3) the distance(s) between the
queried point and the building(s) along the sun rays projection on the
ground. Figure 1 depicts a scenario
where shadow is being cast by building A onto the facade of building B,
given the solar position defined by its elevation angle \(\alpha_{elev}\)
and azimuth angle \(\alpha_{az}\). Once the intersection point is
identified (marked with x in Figure
1), shadow height (\(h_{shadow}\)) at
the queried point (\(viewer\)) can be calculated based on (1) sun
elevation (\(\alpha_{elev}\)), (2) the height of building A (\(h_{build}\))
and (3) the distance (\(dist_{1}\)) between the \(viewer\) and intersection
point x (Equation (1)).

\[\label{eq:eq_shadow_h} h_{shadow} = h_{build} - dist_{1} \cdot tan(\alpha_{elev}) \tag{1}\]

The latter approach can be extended to the general case of shadow height calculation at any ground location and given any configuration of obstacles. For example, if there is more than one obstacle potentially casting shadow on the queried location, we can calculate \(h_{shadow}\) for each obstacle and then take the maximum value.

Once the shadow height is determined, we may evaluate whether any given 3D point is in shadow or not. This is done simply by comparing the Z-coordinate (i.e. height) of the queried point with the calculated shadow height at the same X-Y (i.e. ground) location.

Instead of calculating shadow height at a pre-specified point (e.g. the \(viewer\) in Figure 1), we can set \(h_{shadow}\) to zero and calculate the distance (\(dist_{2}\)) where the shadow intersects ground level (Equation (2)).

\[\label{eq:eq_shadow_footprint} dist_{2} = \frac{h_{build}}{tan(\alpha_{elev})} \tag{2}\]

Shifting the obstacle outline by the resulting distance (\(dist_{2}\)) in
a direction opposite to sun azimuth (\(\alpha_{az}\)) yields a **shadow
footprint** outline (Weisthal 2014). Shadow footprints are useful to
calculate the exact ground area that is shaded at specific time. For
example, Figure 2 shows the shadow
footprints produced by a single building at different times of a given
day.

The **Sky View Factor**
(Grimmond et al. 2001; Erell et al. 2011; Beckers 2013) is the extent
of sky observed from a point as a proportion of the entire sky
hemisphere. The SVF can be calculated based on the maximal angles
(\(\beta\)) formed in triangles defined by the queried location and the
obstacles (Figure 3), evaluated in multiple
circular cross-sections surrounding the queried location. Once the
maximal angle \(\beta_{i}\) is determined for a given angular section \(i\),
\(SVF_{i}\) for that particular section is defined (Gál and Unger 2014) in
Equation (3).

\[\label{eq:eq_svf} SVF_{i} = 1 - sin^2(\beta_{i}) \tag{3}\]

For example, in case (\(\beta_{i}=45^\circ\)), as depicted in Figure 3, \(SVF_{i}\) is equal to -

\[SVF_{i} = 1 - sin^2(45^\circ) = 0.5\]

Averaging \(SVF_{i}\) values for all \(i=1,2,...,n\) circular cross-sections gives the final \(SVF\) estimate for the queried location (Equation (4)).

\[\label{eq:eq_svf_final} SVF = \frac{\sum_{i=1}^{n}SVF_{i}}{n} \tag{4}\]

The number of evaluated cross sections depends on the chosen angular resolution. For example, an angular resolution of \(5^\circ\) means the number of cross sections is \(n=360^\circ/5^\circ=72\) (Figure 4).

Frequently, evaluating whether a given location is shaded, and when, is just a first step towards evaluating the amount of solar radiation for a given period of time. The annual insolation at a given point is naturally affected by the degree of shading throughout the year, but shading is not the only factor.

The three components of the solar radiation are the **direct**,
**diffuse** and **reflected** radiation -

**Direct**radiation refers to solar radiation traveling on a straight line from the sun to the surface of the earth. Direct radiation can be estimated by taking into account: (1) shading, (2) surface orientation relatively to the sun, and (3) meteorological measurements of direct radiation on a horizontal plane or on a plane normal to the beam of sunlight.**Diffuse**radiation refers to solar radiation reaching the Earth’s surface after having been scattered from the direct solar beam by molecules or particulates in the atmosphere. Diffuse radiation can be estimated by taking into account: (1) SVF, and (2) meteorological measurements of diffuse radiation at an exposed location.**Reflected**radiation refers to the sunlight that has been reflected off non-atmospheric obstacles such as ground surface cover or buildings. Most urban surfaces have a low albedo: asphalt reflects only 5-10 percent of incident solar radiation, brick and masonry 20-30 percent, and vegetation about 20 percent. Because a dense urban neighborhood will typically experience multiple reflections, an iterative process is required for a complete analysis. Calculating reflected radiation requires taking into account reflective properties of the various surfaces, their geometrical arrangement (Givoni 1998) and their view factors from the receiving surface, which is beyond the scope of the*shadow*package.

The diffuse radiation component is the dominant one on overcast days, when most radiation is scattered, while the direct radiation component is dominant under clear sky conditions when direct radiation reaches the earth’s surface.

Equation (5) specifies the Coefficient of Direct Normal
Irradiance for a vertical **facade** surface, as function of solar
position given by the difference between facade azimuth and sun azimuth
angles, and sun elevation angle, at time \(t\).

\[\label{eq:coef} {\theta_{facade,t}} = cos({\alpha_{az,t}} - {\alpha'_{az}}) \cdot cos({\alpha_{elev,t}}) \tag{5}\]

In Equation (5), \(\theta_{facade,t}\) is the Coefficient of Direct Normal Irradiance on a facade at time \(t\), \(\alpha_{az,t}\) is the sun azimuth angle at time \(t\) (see Figure 1), \(\alpha'_{az}\) is the facade azimuth angle, i.e. the direction where the facade is facing, and \(\alpha_{elev,t}\) is sun elevation angle at time \(t\) (see Figure 1). Note that all of latter variables, with the exception of facade azimuth angle \(\alpha'_{az}\), are specific for the time interval \(t\) due to the variation in solar position.

Horizontal **roof** surfaces, unlike facades, are not tilted towards any
particular azimuth^{1}. Equation (5) thus simplifies to
Equation (6) when referring to a roof, rather than a
facade, surface.

\[\label{eq:coef_roof} {\theta_{roof,t}} = cos(90^\circ - {\alpha_{elev,t}}) \tag{6}\]

Figure 5 demonstrates the relation given in Equations
(5) and (6) for the entire relevant range of
solar positions relative to facade or roof orientation. Again, note that
for roof surfaces, the \({\theta_{roof,t}}\) coefficient is only dependent
on sun elevation angle \({\alpha_{elev,t}}\) (Equation
(6)) as illustrated on the **right** panel of Figure
5. (The code for producing Figure 5
can be found in the help page of function `coefDirect`

from *shadow*).

For example, the left panel in Figure 5 shows that maximal proportion of incoming solar radiation (i.e. \(\theta_{facade,t} = 1\)) on a facade surface is attained when facade azimuth is equal to sun azimuth and sun elevation is 0 (\({\alpha_{elev,t}}=0^\circ\), i.e. facade directly facing the sun). Similarly, the right panel shows that maximal proportion of solar radiation on a roof surface (i.e. \(\theta_{roof,t} = 1\)) is attained when the sun is at the zenith (\({\alpha_{elev,t}}=90^\circ\), i.e. sun directly above the roof).