shadow (version 0.7.1)

shadowHeight: Shadow height calculation considering sun position and obstacles

Description

This function calculates shadow height at given points or complete grid (location), taking into account:

• Obstacles outline (obstacles), given by a polygonal layer with a height attribute (obstacles_height_field)

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

Usage

# S4 method for SpatialPoints
shadowHeight(
  location,
  obstacles,
  obstacles_height_field,
  solar_pos = solarpos2(location, time),
  time = NULL,
  b = 0.01,
  parallel = getOption("mc.cores"),
  filter_footprint = FALSE
)
# S4 method for Raster
shadowHeight(
  location,
  obstacles,
  obstacles_height_field,
  solar_pos = solarpos2(pnt, time),
  time = NULL,
  b = 0.01,
  parallel = getOption("mc.cores"),
  filter_footprint = FALSE
)

Arguments

location

A SpatialPoints* or Raster* object, specifying the location(s) for which to calculate 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 decimal degrees from North), second column is sun elevation (in decimal degrees); rows represent different positions (e.g. at different times of day)

time

When solar_pos is unspecified, time can be passed to automatically calculate solar_pos 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

b

Buffer size when joining intersection points with building outlines, to determine intersection height

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

filter_footprint

Should the points be filtered using shadowFootprint before calculating shadow height? This can make the calculation faster when there are many point which are not shaded

Value

Returned object is either a numeric matrix or a Raster* -

• 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 height

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

In both cases the numeric values express shadow height -

• NA value means no shadow

• A valid number expresses shadow height, in CRS units (e.g. meters)

• Inf means complete shadow (i.e. sun below horizon)

Examples

# NOT RUN {
# Single location
location = rgeos::gCentroid(build)
location_geo = matrix(c(34.7767978098526, 31.9665936050395), ncol = 2)
time = as.POSIXct("2004-12-24 13:30:00", tz = "Asia/Jerusalem")
solar_pos = maptools::solarpos(location_geo, time)
plot(build, main = time)
plot(location, add = TRUE)
sun = shadow:::.sunLocation(location = location, sun_az = solar_pos[1,1], sun_elev = solar_pos[1,2])
sun_ray = ray(from = location, to = sun)
build_outline = as(build, "SpatialLinesDataFrame")
inter = rgeos::gIntersection(build_outline, sun_ray)
plot(sun_ray, add = TRUE, col = "yellow")
plot(inter, add = TRUE, col = "red")
shadowHeight(
  location = location,
  obstacles = build,
  obstacles_height_field = "BLDG_HT",
  solar_pos = solar_pos
)

# Automatically calculating 'solar_pos' using 'time'
proj4string(build) = CRS("+init=epsg:32636")
proj4string(location) = CRS("+init=epsg:32636")
shadowHeight(
  location = location,
  obstacles = build,
  obstacles_height_field = "BLDG_HT",
  time = time
)

# }
# NOT RUN {
# Two points - three times
location0 = rgeos::gCentroid(build)
location1 = raster::shift(location0, 0, -15)
location2 = raster::shift(location0, -10, 20)
locations = rbind(location1, location2)
time = as.POSIXct("2004-12-24 13:30:00", tz = "Asia/Jerusalem")
times = seq(from = time, by = "1 hour", length.out = 3)
shadowHeight(                            ## Using 'solar_pos'
  location = locations,
  obstacles = build,
  obstacles_height_field = "BLDG_HT",
  solar_pos = maptools::solarpos(location_geo, times)
)
shadowHeight(                            ## Using 'time'
  location = locations,
  obstacles = build,
  obstacles_height_field = "BLDG_HT",
  time = times
)

# Grid - three times
time = as.POSIXct("2004-12-24 13:30:00", tz = "Asia/Jerusalem")
times = seq(from = time, by = "1 hour", length.out = 3)
ext = as(raster::extent(build), "SpatialPolygons")
r = raster::raster(ext, res = 2)
proj4string(r) = proj4string(build)
x = Sys.time()
shadow1 = shadowHeight(
  location = r,
  obstacles = build,
  obstacles_height_field = "BLDG_HT",
  time = times,
  parallel = 3
)
y = Sys.time()
y - x
x = Sys.time()
shadow2 = shadowHeight(
  location = r,
  obstacles = build,
  obstacles_height_field = "BLDG_HT",
  solar_pos = solarpos2(r, times),
  parallel = 3
)
y = Sys.time()
y - x
shadow = shadow1
opar = par(mfrow = c(1, 3))
for(i in 1:raster::nlayers(shadow)) {
  plot(shadow[[i]], col = grey(seq(0.9, 0.2, -0.01)), main = raster::getZ(shadow)[i])
  raster::contour(shadow[[i]], add = TRUE)
  plot(build, border = "red", add = TRUE)
}
par(opar)

# }