shadow (version 0.7.1)

inShadow: Logical shadow calculation (is given point shaded?) for 3D points considering sun position and obstacles

Description

This function determines whether each given point in a set of 3D points (location), is shaded or not, taking into account:

• Obstacles outline (obstacles), given by a polygonal layer with a height attribute (obstacles_height_field), or alternatively a Raster* which is considered as a grid of ground locations

• Sun position (solar_pos), given by azimuth and elevation angles

Alternatively, the function determines whether each point is in shadow based on a raster representing shadow height shadowHeightRaster, in which case obstacles, obstacles_height_field and solar_pos are left unspecified.

Usage

# S4 method for SpatialPoints,Raster,missing,missing
inShadow(
  location,
  shadowHeightRaster,
  obstacles,
  obstacles_height_field,
  solar_pos
)
# S4 method for SpatialPoints,missing,ANY,ANY
inShadow(
  location,
  shadowHeightRaster,
  obstacles,
  obstacles_height_field,
  solar_pos = solarpos2(location, time),
  time = NULL,
  ...
)
# S4 method for Raster,missing,ANY,ANY
inShadow(
  location,
  shadowHeightRaster,
  obstacles,
  obstacles_height_field,
  solar_pos = solarpos2(pnt, time),
  time = NULL,
  ...
)

Arguments

location

A SpatialPoints* or Raster* object, specifying the location(s) for which to calculate logical shadow values. If location is SpatialPoints*, then it can have 2 or 3 dimensions. A 2D SpatialPoints* is considered as a point(s) on the ground, i.e. 3D point(s) where \(z=0\). In a 3D SpatialPoints* the 3rd dimension is assumed to be elevation above ground \(z\) (in CRS units). Raster* cells are considered as ground locations

shadowHeightRaster

Raster representing shadow height

obstacles

A SpatialPolygonsDataFrame object specifying the obstacles outline

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 degrees from North), second column is sun elevation (in degrees); rows represent different positions (e.g. at different times of day)

time

When both shadowHeightRaster and solar_pos are unspecified, time can be passed to automatically calculate solarpos based on the time and the centroid of location, using function maptools::solarpos. In such case location must have a defined CRS (not NA). The time value must be a POSIXct or POSIXlt object.

...

Other parameters passed to shadowHeight, such as parallel

Value

Returned object is either a logical matrix or a Raster* with logical values -

• If input location is a SpatialPoints*, then returned object is a matrix where rows represent spatial locations (location features), columns represent solar positions (solar_pos rows) and values represent shadow state

• If input location is a Raster*, then returned object is a RasterLayer or RasterStack, where raster layers represent solar positions (solar_pos rows) and pixel values represent shadow state

In both cases the logical values express shadow state:

• TRUE means the location is in shadow

• FALSE means the location is not in shadow

• NA means the location 3D-intersects an obstacle

Examples

# NOT RUN {
# Method for 3D points - Manually defined

opar = par(mfrow = c(1, 3))

# Ground level
location = sp::spsample(
  rgeos::gBuffer(rgeos::gEnvelope(build), width = 20),
  n = 80,
  type = "regular"
)
solar_pos = as.matrix(tmy[9, c("sun_az", "sun_elev")])
s = inShadow(
  location = location,
  obstacles = build,
  obstacles_height_field = "BLDG_HT",
  solar_pos = solar_pos
)
plot(location, col = ifelse(s[, 1], "grey", "yellow"), main = "h=0")
plot(build, add = TRUE)

# 15 meters above ground level
coords = coordinates(location)
coords = cbind(coords, z = 15)
location1 = SpatialPoints(coords, proj4string = CRS(proj4string(location)))
solar_pos = as.matrix(tmy[9, c("sun_az", "sun_elev")])
s = inShadow(
  location = location1,
  obstacles = build,
  obstacles_height_field = "BLDG_HT",
  solar_pos = solar_pos
)
plot(location, col = ifelse(s[, 1], "grey", "yellow"), main = "h=15")
plot(build, add = TRUE)

# 30 meters above ground level
coords = coordinates(location)
coords = cbind(coords, z = 30)
location2 = SpatialPoints(coords, proj4string = CRS(proj4string(location)))
solar_pos = as.matrix(tmy[9, c("sun_az", "sun_elev")])
s = inShadow(
  location = location2,
  obstacles = build,
  obstacles_height_field = "BLDG_HT",
  solar_pos = solar_pos
)
plot(location, col = ifelse(s[, 1], "grey", "yellow"), main = "h=30")
plot(build, add = TRUE)

par(opar)

# Shadow on a grid covering obstacles surface
# }
# NOT RUN {
# Method for 3D points - Covering building surface

obstacles = build[c(2, 4), ]
location = surfaceGrid(
  obstacles = obstacles,
  obstacles_height_field = "BLDG_HT",
  res = 2,
  offset = 0.01
)
solar_pos = tmy[c(9, 16), c("sun_az", "sun_elev")]
solar_pos = as.matrix(solar_pos)
s = inShadow(
  location = location,
  obstacles = obstacles,
  obstacles_height_field = "BLDG_HT",
  solar_pos = solar_pos
)
location$shadow = s[, 1]
plotGrid(location, color = c("yellow", "grey")[as.factor(location$shadow)], size = 0.5)
location$shadow = s[, 2]
plotGrid(location, color = c("yellow", "grey")[as.factor(location$shadow)], size = 0.5)

# Method for ground locations raster

ext = as(raster::extent(build) + 20, "SpatialPolygons")
location = raster::raster(ext, res = 2)
proj4string(location) = proj4string(build)
obstacles = build[c(2, 4), ]
solar_pos = tmy[c(9, 16), c("sun_az", "sun_elev")]
solar_pos = as.matrix(solar_pos)
s = inShadow(                ## Using 'solar_pos'
  location = location,
  obstacles = obstacles,
  obstacles_height_field = "BLDG_HT",
  solar_pos = solar_pos,
  parallel = 3
)
time = as.POSIXct(tmy$time[c(9, 16)], tz = "Asia/Jerusalem")
s = inShadow(               ## Using 'time'
  location = location,
  obstacles = obstacles,
  obstacles_height_field = "BLDG_HT",
  time = time,
  parallel = 3
)
plot(s)

# Method for pre-calculated shadow height raster

ext = as(raster::extent(build), "SpatialPolygons")
r = raster::raster(ext, res = 1)
proj4string(r) = proj4string(build)
r[] = rep(seq(30, 0, length.out = ncol(r)), times = nrow(r))
location = surfaceGrid(
  obstacles = build[c(2, 4), ],
  obstacles_height_field = "BLDG_HT",
  res = 2,
  offset = 0.01
)
s = inShadow(
  location = location,
  shadowHeightRaster = r
)
location$shadow = s[, 1]
r_pnt = raster::as.data.frame(r, xy = TRUE)
coordinates(r_pnt) = names(r_pnt)
proj4string(r_pnt) = proj4string(r)
r_pnt = SpatialPointsDataFrame(
  r_pnt,
  data.frame(
    shadow = rep(TRUE, length(r_pnt)),
    stringsAsFactors = FALSE
    )
 )
pnt = rbind(location[, "shadow"], r_pnt)
plotGrid(pnt, color = c("yellow", "grey")[as.factor(pnt$shadow)], size = 0.5)

# Automatically calculating 'solar_pos' using 'time' - Points
location = sp::spsample(
  rgeos::gBuffer(rgeos::gEnvelope(build), width = 20),
  n = 500,
  type = "regular"
)
time = as.POSIXct("2004-12-24 13:30:00", tz = "Asia/Jerusalem")
s = inShadow(
  location = location,
  obstacles = build,
  obstacles_height_field = "BLDG_HT",
  time = time
)
plot(location, col = ifelse(s[, 1], "grey", "yellow"), main = time)
plot(build, add = TRUE)

# }
# NOT RUN {
# }