inShadow: Logical shadow calculation (is given point shaded?) for 3D...

In michaeldorman/shadow: Geometric Shadow Calculations

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 signature 'SpatialPoints,Raster,missing,missing'
inShadow(
  location,
  shadowHeightRaster,
  obstacles,
  obstacles_height_field,
  solar_pos
)

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

## S4 method for signature '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

Note

For a correct geometric calculation, make sure that:

• The layers location and obstacles are projected and in same CRS

• The values in obstacles_height_field of obstacles are given in the same distance units as the CRS (e.g. meters when using UTM)

Examples

# Method for 3D points - Manually defined

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

# Ground level
location = sp::spsample(
  # rgeos::gBuffer(rgeos::gEnvelope(build), width = 20),
  as(sf::st_buffer(sf::st_as_sfc(sf::st_bbox(sf::st_as_sf(build))), dist = 20), "Spatial"),
  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),
  as(sf::st_buffer(sf::st_as_sfc(sf::st_bbox(sf::st_as_sf(build))), dist = 20), "Spatial"),
  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)


## End(Not run)

michaeldorman/shadow documentation built on Sept. 10, 2023, 4:17 a.m.