biophys: Spatially-explicit biophysical models

#' 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!"))
#   }
# }