#' Calculates operative body temperature
#' @param biophys.inputs Object created from running \code{\link[biophys]{biophys.prep}}
#'
#' @usage op.body.temp(biophys.inputs)
#'
#' @examples
#' # Create example data
#' r.max <- raster(ncol=10, nrow=10)
#'
#' # assign values to cells
#' r.max[] <- runif(ncell(r.max),10,20)
#' r.min <- r.max - 5
#' r.elev <- r.max * 10
#'
#' # Make vector of hours
#' night.hours <- c(19:24,1:7)
#'
#' biophys.input <- biophys.prep(r.max, r.min, r.elev, hour = night.hours, month = c("April"), Julian = c(106))
#' temp.data <- op.body.temp(biophys.input)
#'
#' @export
#' @author Bill Peterman <Bill.Peterman@@gmail.com>
#' @return By default, this function will return a list of hourly operative body temperatures for each raster pixel. Additionally, calculated parameters needed for subsequent functions are also saved. The list object generated by this function is formatted for input into \code{\link[biophys]{active.evap}}
op.body.temp <- function(biophys.inputs) {
# Unpack inputs
Tmax.dat = biophys.inputs$Tmax
type <- biophys.inputs$type
Tmin.dat = biophys.inputs$Tmin
elev.dat = biophys.inputs$elev
NoData.value = biophys.inputs$NoData.value
hours = biophys.inputs$hours
an.height = biophys.inputs$an.height
wind.height = biophys.inputs$wind.height
u = biophys.inputs$u
pCp = biophys.inputs$pCp
Julian = biophys.inputs$Julian
month = biophys.inputs$month
RH = biophys.inputs$RH
length.m = biophys.inputs$length.m
phi = biophys.inputs$phi
SM = biophys.inputs$SM
sigma = biophys.inputs$sigma
e = biophys.inputs$e
es = biophys.inputs$es
cp = biophys.inputs$cp
alphaS = biophys.inputs$alphaS
alphaL = biophys.inputs$alphaL
Fp = biophys.inputs$Fp
Fd.c = biophys.inputs$Fd.c
Fa.c = biophys.inputs$Fa.c
Fr.c = biophys.inputs$Fr.c
Fg.c = biophys.inputs$Fg.c
# Bioenergetic model for P. jordani (Energy Budget)
#=========================================================================
layers <- length(month)
all.results <- vector(mode = "list",length = layers)
for(i in 1:layers){
if(layers>1 & type=="list"){
Tx <- raster(Tmax.dat[i])
Tm <- raster(Tmin.dat[i])
raster.dat <- stack(Tx, Tm, elev.dat)
names(raster.dat) <- c("Tx","Tmin","elev")
} else if(layers>1 & type=="stack") {
Tx <- Tmax.dat[[i]]
Tm <- Tmin.dat[[i]]
raster.dat <- stack(Tx, Tm, elev.dat)
names(raster.dat) <- c("Tx","Tmin","elev")
} else {
raster.dat <- stack(Tmax.dat, Tmin.dat, elev.dat)
names(raster.dat) <- c("Tx","Tmin","elev")
}
NAvalue(raster.dat) <- NoData.value # Specify the NoData value
# plot(raster.dat) # Make sure things look right
# Convert rasters to vectors
Tx <- as.vector(raster.dat[[1]])
Tmin <- as.vector(raster.dat[[2]])
elev <- as.vector(raster.dat[[3]])
# Remove NA values
Tx <- Tx[!is.na(Tx)]
Tmin <- Tmin[!is.na(Tmin)]
elev <- elev[!is.na(elev)]
# Combine into single data frame, filter out NoData values
# vector.dat <- data.frame(Tx,Tmin=Tmin$Tmin,elev=elev$elev) %>% filter(!is.na(Tx),!is.na(Tmin),!is.na(elev))
Tn=Tx-Tmin #Tn=daily temperature range (Celcius)
# Tavg=(Tx+Tmin)/2
Tsa=(((1.39734+0.88841*Tx)+273.15)+((1.39734+0.88841*Tmin)+273.15))/2 #mean soil temperature from low and high points on transect (with adiabatic cooling rate of 0.575 degrees per 100m)
Ad=Tn/2 #variable used in one of the time functions below
remove(Tmin) #removes Tmin matrix from memory
#######################################################
#######################################################
# Start for-loop to calculate values for each hour
Te.hr <- vector(mode = "list",length = length(hours)) # Empty list to store hourly body temp results in
Ta.hr <- vector(mode = "list",length = length(hours)) # Empty list to store hourly air temp results in
for(j in seq_along(hours)){
hr <- hours[j]
# RH=0.95 #relative humidity (assumed to be constant); not a great assumption, but perhaps conservative
# d1=0.055 #characteristic dimension (i.e., SVL of an adult jordani, in meters)
# mass1=4599.1*d1^2.5297 #equation to predict the mass of the salamander from the SVL
#Day Length and Sun Angles
#=====================================
delta=asin(0.39795*cos(0.21631+2*atan(0.967*tan(0.0086*(-186+Julian))))) #delta=solar declination (in radians); Julian=calendar day (Jan 1 = 1)
#Day Length
#=====================================
#hd=day length
#phi=latitude (in radians)
phi=0.623
hd=24-(24/pi)*acos((sin(6*pi/180)+sin(phi)*sin(delta))/cos(phi)*cos(delta))
#time of solar noon
#=====================================
#t0=solar noon
f=(pi/180)*(279.575+0.9856*Julian)
#ET=Equation of Time
ET=((-104.7*sin(f))+(596.2*sin(2*f))+(4.3*sin(3*f))-(12.7*sin(4*f))-(429.3*cos(f)-(2.0*cos(2*f)+(19.3*cos(3*f)))))/3600
#LC=correction for each degree that a location is east of the standard meridian
#SM=standard meridian (0, 15, ...., 345 degrees) in radians??
# SM=1.5708
LC=0.1315*(1/15)+SM
t0=12-LC-ET
#Zenith angle, psi (radians), is the sun angle measured from vertical
#psi=Zenith Angle==========psi is actually cos(psi)
#hr=hour of day from above
psi=acos(sin(delta)*sin(phi)+cos(delta)*cos(phi)*(cos(pi/(12*(hr-t0)))))
#Air and Soil Temperature
#==========================================================================
#Measurement of air temperature at any particular hour
#w=pi/12
#Gamma=
w=pi/12
Gam=0.44-0.46*sin(0.9+(w*hr))+0.11*sin(0.9+(2*w*hr))
#Hourly air temperature
#Ta=air temperature
Ta=Tx-Tn*(1-Gam)
# remove(Tn) #remove temp range matrix from memory
# remove(Tx)
# remove(elev)
#Hourly soil surface temperature
#Ts=soil surface temperature (Kelvin)
Ts=Tsa+Ad*sin(w*(hr-8)) #The time variable in the sine function is phase adjusted by 8 hours
#Radiation and Environmental Temperature
# sigma=5.67e-8 #(Stefan-Boltzmann constant)
B=sigma*(Ta+273.15)^4 #emitted flux density (W m^2)
# Assume gray body emissivity (no wavelength dependency)
# e=0.93
# es=0.96
eac=9.2e-6*(Ta+273.15)^2 #clear sky emissivity (approximation from Swinbank 1963)
#Modeling wind speed at animal level (0.5 cm)
aheight=rep(an.height, length(Ta)) #height of animal in meters
z=rep(wind.height, length(Ta)) #height of wind speed measurement in meters
ds=0.65*aheight
zm=0.1*aheight
# u=1.0 #measured wind speed at measurement height (from instrument)
ustar=0.4*(u/(log((z-ds)/zm))) #friction velocity
u=(ustar/0.4)*log(aheight*(1-0.65)/zm) #wind speed at height of interest (m s^-1)
remove(z)
remove(aheight)
remove(ds)
remove(zm)
# cp=29.3 #specific heat of air
gHa=1.4*0.135*sqrt(u/length.m) #boundary conductance of air (mol/square meter/second)-species1
gr=(4*e*sigma*(Ta+273.15)^3)/cp #radiative conductance (mol/square meter/second)
#Longwave Component
##############################################
La=eac*sigma*(Ta+273.15)^4 #longwave flux density from atmosphere (W/square meter)
Lg=es*sigma*Ts^4 #longwave flux density from the ground (W/square meter)
# alphaS=0.93 #absorptivity in solar waveband
# alphaL=0.96 #absorptivity in thermal waveband
# Fp=0 # ctor (ratio of projected area perpendicular to solar beam) assumed zero because of nocturnal activity
Fd=rep(Fd.c, length(Ta)) #estimated for a standing lizard (Bartlett and Gates 1967)
Fr=rep(Fr.c, length(Ta)) # I made these matrices for the calculations (would need to change these to constants)
Fa=rep(Fa.c, length(Ta)) # these are just estimated proportions of the animal exposed to different types of temperature exchange
Fg=rep(Fg.c, length(Ta))
#Absorbed Energy
#################################################
Rabs=alphaL*(Fa*La+Fg*Lg) #Only conduction!!!!!!!
#Radiation emitted
Remit=(es*sigma*(Ta-273.15)^4)
#Operative body temperature
#############################################
# Te=Ta+((Rabs-Remit1)/cp*(gr+gHa))
# Store hourly results
Te.hr[[j]] <- Ta+((Rabs-Remit)/cp*(gr+gHa))
Ta.hr[[j]] <- Ta
names(Te.hr[[j]]) <- paste0(month[i],"_hr",hours[j])
names(Ta.hr[[j]]) <- paste0(month[i],"_hr",hours[j])
} # Close hours loop
#### Code below this line
# Get average of hourly values
# op.temp <- rowMeans(simplify2array(bt.hr))
# Put results back into raster
# op.temp_rast <- raster.dat[[1]]
# tmp <- as.vector(raster.dat$Tx)
# tmp <- replace(tmp,which(tmp!=NoData.value),op.temp)
# op.temp_rast <- setValues(op.temp_rast,tmp)
# Run function to package data
hr.dat <-package.dat(Te = Te.hr,
Ta = Ta.hr,
hours = hours,
month = month[i]
)
# Combine all data objects for export
# dat.out <- list(hr.dat = hr.dat,
# raster.dat = list(
# rast.extent = extent(raster.dat[[1]]),
# rast.res = res(raster.dat[[1]]),
# full.vector = as.vector(raster.dat[[1]]) # This object will be used to recompile a raster
# ))
# op.dat <- list(dat = dat.out,
# parmlist = parmlist)
all.results[[i]] <- hr.dat # Store results of monthly iteration
} # Close multiple layer for loop
names(all.results) <- month
return(all.results)
} # End function
#
#
#
#
#
#
#
#
package.dat <- function(Te, Ta, hours, month){
Te.list <- list()
Ta.list <- list()
for(i in seq_along(hours)){
hr <- hours[i]
Te.hr <- Te[[i]]
Ta.hr <- Ta[[i]]
# tmp.Te <- replace(vector.full,which(vector.full!=NoData.value),Te.hr)
# tmp.Ta <- replace(vector.full,which(vector.full!=NoData.value),Ta.hr)
Te[[i]] <- as.vector(Te.hr)
Ta[[i]] <- as.vector(Ta.hr)
# names(Te.list[[i]]) <- paste0(month,"_hr",hours[i])
# names(Ta.list[[i]]) <- paste0(month,"_hr",hours[i])
}
names(Te) <- paste0(month,"_hr",hours)
names(Ta) <- paste0(month,"_hr",hours)
Temp.dat <- list(Te=Te,
Ta=Ta)
return(Temp.dat)
}
#
#
#
#
#
# No longer using complex function below
# Function to export hourly temperature data:
# make.rast <- function(Te, Ta, hours, rast, NoData.value, vector.full, out, grid.type, month){
# if(out=='raster'){
# raster.stack <- stack()
# for(i in seq_along(hours)){
# hr <- hours[i]
# op.hr <- Te[[i]]
# tmp <- replace(vector.full,which(vector.full!=NoData.value),op.hr)
# op.temp_rast <- setValues(rast,tmp)
# raster.stack <- stack(raster.stack,op.temp_rast)
# names(raster.stack) <- paste0(month,"_hr",hours)
# }
# return(raster.stack)
# } else if(out=='list'){
# Te.list <- list()
# Ta.list <- list()
# for(i in seq_along(hours)){
# hr <- hours[i]
# Te.hr <- Te[[i]]
# Ta.hr <- Ta[[i]]
# tmp.Te <- replace(vector.full,which(vector.full!=NoData.value),Te.hr)
# tmp.Ta <- replace(vector.full,which(vector.full!=NoData.value),Ta.hr)
# Te.list[[i]] <- tmp.Te
# Ta.list[[i]] <- tmp.Ta
# # names(Te.list[[i]]) <- paste0(month,"_hr",hours[i])
# # names(Ta.list[[i]]) <- paste0(month,"_hr",hours[i])
# }
# }
# names(Te.list) <- paste0(month,"_hr",hours)
# names(Ta.list) <- paste0(month,"_hr",hours)
# return(hr.list)
# } else if(file.info(export)$isdir){
# for(i in seq_along(hours)){
# hr <- hours[i]
# op.hr <- results[[i]]
# tmp <- replace(vector.full,which(vector.full!=NoData.value),op.hr)
# op.temp_rast <- setValues(rast,tmp)
# writeRaster(op.temp_rast,paste0(out,month,"_hr",hr,grid.type),overwrite=TRUE)
# }
# } else {
# warning(paste0(out," is not recognized!"))
# }
# }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.