shadow: R Package for Geometric Shadow Calculations in an Urban Environment

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.

Michael Dorman (BGU) , Evyatar Erell (BGU) , Adi Vulkan (BGU) , Itai Kloog (BGU)
2019-08-17

1 Introduction

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 -

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 -

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.

2 Theory

Shadow height

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

graphic without alt text
Figure 1: Shadow height calculation

\[\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.

Logical shadow flag

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.

Shadow footprint

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.

graphic without alt text
Figure 2: Shadow footprints cast by a building on a horizontal ground surface at hourly intervals on 2004-06-24. The building, indicated by the gray shaded area, is located at 31.97°N 34.78°E, and is 21.38 meters tall

Sky View Factor (SVF)

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

graphic without alt text
Figure 3: Sky View factor calculation

\[\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).

graphic without alt text
Figure 4: Angular cross sections for calculating the Sky View Factor (SVF)

Solar radiation

Components

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 -

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.

Direct Normal Irradiance

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 azimuth1. 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).