runauto2: Function for automatically generating microclimate surfaces...

In ilyamaclean/microclima: microclimate modelling with R

Function for automatically generating microclimate surfaces for areas with variable habitat

Description

This function generating microclimate temperature surfaces from first principles. It can be used in place of runauto() when the region for which microclimate surfaces are wanted contains variable habitat. No mesoclimate adjustments are performed using this function. In contrast to runauto() bespoke parameterisations are performed for each pixel though the methods for computing soil and latent heat fluxes are more simplistic. If hourly weather data are not provided, it downloads coarse-resolution climate and radiation data from the NCEP-NCAR or NCEP–DOE Atmospheric Model Intercomparison Project (Kanamitso et al 2002) and interpolates these data to hourly. It runs the microclimate model in hourly time increments to generate an array of temperatures, relative humidities and wind speeds. Incoming solar and downward longwave radiation are also returned.

Usage

runauto2(
  dem,
  hourlydata = NA,
  dstart,
  dfinish,
  lat = NA,
  long = NA,
  hgt = 0.05,
  l,
  x,
  veghgt,
  alb,
  wind.agg = 10,
  merid = 0,
  dst = 0,
  plot.progress = TRUE
)

Arguments

dem	a raster object of elevations. Also defines extent and resolution of returned microclimate data. Supplied raster must have a projection such that the units of x, y and z are identical, and grid cells are square. NA values are assumed to be sea, with elevation zero
hourlydata	an optional dataframe of hourly climate forcing variables. If not supplied downloaded using hourlyNCEP() obs_time POSIXlt object of times in UTC temperature Temperatures at 2 m (ºC) humidity Specific humidity at 2m (Kg / Kg) pressure Surface pressure (Pa) windspeed Wind speed at 2 m (metres per second) winddir Wind direction (degrees from N) emissivity Emissivity of the atmosphere (0 - 1, downlong / uplong) netlong Net longwave radiation (MJ m-2 hr-1) uplong Upward longwave radiation (MJ m-2 hr-1) downlong Downward longwave radiation (MJ m-2 hr-1) rad_dni Direct radiation normal to the solar beam (MJ m-2 hr-1) rad_dif Diffuse radiation (MJ m-2 hr-1) szenith The zenith angle (degrees) cloudcover Cloud cover (Percentage)
dstart	start date as character of time period required in format DD/MM/YYYY. Ignored if hourly weather data supplied
dfinish	end date as character of time period required in format DD/MM/YYYY. Ignored if hourly weather data supplied
hgt	the height (in m) above or below ground for which temperature estimates are required. Below-ground not supported. Must be less than two metres.
l	a raster, matrix or array of leaf area index values, as for example returned by lai() (see details)
x	a raster, matrix or array of representing the ratio of vertical to horizontal projections of foliage as, for example, returned by leaf_geometry(). (see details)
veghgt	a raster, matrix or array of vegetation heights in metres (see details)
alb	a raster, matrix or array of surface albedo in range 0 to 1 (see details)
wind.agg	number of cells over which to average wind in calculation of vertical wind profiles must be lower than dimensions of dem
merid	an optional numeric value representing the longitude (decimal degrees) of the local time zone meridian (0 for GMT).
dst	an optional numeric value representing the time difference from the timezone meridian (hours, e.g. +1 for BST if merid = 0).
plot.progress	optional logial indicating whether to plot temperatures every 100 hours.

Details

This function uses K-theory to determine wind and temperature profiles above ground as a function of sensible heat flux, a roughness length and zero-plane displacement height. Sensible heat is determined from net radiation and simple approximations of latent heat and ground heat fluxes. Roughness lengths and zero-plane displacement heights are determined from vegetation height and leaf area. If hgt is below veghgt, the above vegetation temperature profile is extrapolated to a point 0.8 x the height of the vegetation as here extrapolated wind speeds, which ultimately govern the temperature profile, are approximately equivelent to the average wind speed within the canopy. Incoming radiation and outgoing radiation are adjusted by the proportion of l above hgt assuming uniform vertical density in vegtation. The returned wind profile values for below canopy are calculated using the method describedin Campbell & Norman (2012) Environmental Biophysics, Springer. Since wind profiles are determined by fetch as well as surface roughness, the vertical wind profiles are spatially averaged by an amount specified by wind.agg. if l, x, veghgt or alb are arrays and the third dimension of the array does not correspond to the number of hours in the time sequence, hourly values are derived by interpolation.

Value

a list with the following objects: (1) temperature: an array of temperatures for each pixel of dem and hour of the time sequence (degrees C. (2) relhum: an array of temperatures for each pixel of dem and hour of the time sequence (Percentage). (3) windspeed an array of temperatures for each pixel of dem and hour of the time sequence (m / s). (4) swdown an array of incoming solar radiation for each pixel of dem and hour of the time sequence (MJ / m^2 / hr). (5) lwdown an array of downward longwave radiation for each pixel of dem and hour of the time sequence (MJ / m^2 / hr).

Examples

library(raster)
# =========================================================================== #
# ~~~~~~~~~~~~~~~~~ Create spatial datasets needed for input ~~~~~~~~~~~~~~~~ #
# =========================================================================== #
e<-extent(c(169450,169550,12450,12550))
dem <- crop(dtm1m,e) # elevation
alb <- albedo(aerial_image[,,1], aerial_image[,,2], aerial_image[,,3],
              aerial_image[,,4]) # albedo
plot(if_raster(alb, dtm1m), col = gray(0:255/255))
# Adjust albedo using modis input
img <- if_raster(alb, dtm1m)
mds <- raster(modis, xmn = 169000, xmx = 170000, ymn = 12000, ymx = 13000)
alb <- albedo_adjust(img, mds)
alb <- crop(alb, e)
plot(alb, col = gray(0:255/255))
# Leaf area index
leaf <- lai(aerial_image[,,3], aerial_image[,,4])
leaf<- raster(leaf,template=dtm1m)
l<-crop(leaf,dem)
# Leaf angle coefficient
x <- leaf_geometry(veg_hgt)
x <-crop(x,e)
veghgt <- crop(veg_hgt, e) # vegetation height
# =========================================================================== #
# ~~~~~~~~~~~ Run model in hourly time-steps ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
# ~~~~~~~~~~~ Weather data automatically downloaded as hourlydata = NA ~~~~~~ #
# =========================================================================== #
microout <- runauto2(dem, hourlydata = NA, dstart = "01/05/2019",
                     dfinish = "31/05/2019", hgt = 0.05, l = l, x = x,
                     veghgt = veghgt, alb = alb)
# =========================================================================== #
# ~~~~~~~~~~~ Extract temperature data and plot results ~~~~~~~~~~~~~~~~~~~~~ #
# =========================================================================== #
temps <- microout$temperature
tmin <- apply(temps,c(1,2),min) # spatial min
tmax <- apply(temps,c(1,2),max) # spatial max
tmin2 <- apply(temps,3,min, na.rm = T) # time series min
tmax2 <- apply(temps,3,max, na.rm = T) # time series max
# Create palette:
mypal <- colorRampPalette(c("darkblue", "blue", "green", "yellow", "orange", "red"))(255)
# Raster plots
par(mfrow=(c(1,2)))
plot(if_raster(tmin,dem), col = mypal)
plot(if_raster(tmax,dem), col = mypal)
# Time trend plots
ymn<-min(tmin2); ymx <- max(tmax2) # limits for plot
tme <- as.POSIXct(c(0:743)*3600,origin="2019-05-01 00:00", tz = "GMT")
par(mfrow = c(1,1))
plot(tmin2 ~ tme, type = "l", col = "blue", ylim = c(ymn, ymx),
     xlab = "", ylab = "")
par(new=T)
plot(tmax2 ~ tme, type = "l", col = "red", ylim = c(ymn, ymx),
     xlab = "Day", ylab = "Temperature")

ilyamaclean/microclima documentation built on Oct. 31, 2023, 11:41 p.m.