AeroEvapR_Tutorial

## Step 1: Load libraries and data

# Load in necessary libraries
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(ggplot2)
library(lubridate)
#> 
#> Attaching package: 'lubridate'
#> The following objects are masked from 'package:base':
#> 
#>     date, intersect, setdiff, union
library(data.table)
#> 
#> Attaching package: 'data.table'
#> The following objects are masked from 'package:lubridate':
#> 
#>     hour, isoweek, mday, minute, month, quarter, second, wday, week,
#>     yday, year
#> The following objects are masked from 'package:dplyr':
#> 
#>     between, first, last
library(AeroEvapR)

# Read in data
file_path <- system.file("extdata", "PadreBay_QAQC_2019_Sep_Oct.csv", package = "AeroEvapR")
df <- read.csv(file_path)

# View data and column names
str(df)
#> 'data.frame':    2928 obs. of  6 variables:
#>  $ date             : chr  "9/1/2019 0:00" "9/1/2019 0:30" "9/1/2019 1:00" "9/1/2019 1:30" ...
#>  $ AirT_C_Avg       : num  28.9 29.2 28.6 28.4 27.8 ...
#>  $ RH_Avg           : num  39.5 32.3 38.1 42.1 40.9 ...
#>  $ Wind_Speed_ms_Avg: num  2.367 2.732 0.891 1.017 1.932 ...
#>  $ SkinT_C_Corr_Avg : num  26.4 26.2 25.9 26.7 27 ...
#>  $ BP_kPa_Avg       : num  88.1 88.1 88.1 88.1 88.1 88.1 88.1 88.1 88.2 88.2 ...

## Step 2: Clean up and format the data for the function

# Select the variables needed for the calculation
old_var_names = c('SkinT_C_Corr_Avg', 'Wind_Speed_ms_Avg', 
                            'BP_kPa_Avg', 'AirT_C_Avg',
                            'RH_Avg', 'date')
# List the new variables that correspond to the function requirements
new_var_names <- c('T_skin', 'WS', 'P', 'T_air', 'RH', 'date')

# Replace names column names with names required for package
setnames(df, old = old_var_names, new = new_var_names)
  
# Convert date to datetime objects
df$date <- as.POSIXct(df$date, format = "%m/%d/%Y %H:%M")

str(df)
#> 'data.frame':    2928 obs. of  6 variables:
#>  $ date  : POSIXct, format: "2019-09-01 00:00:00" "2019-09-01 00:30:00" ...
#>  $ T_air : num  28.9 29.2 28.6 28.4 27.8 ...
#>  $ RH    : num  39.5 32.3 38.1 42.1 40.9 ...
#>  $ WS    : num  2.367 2.732 0.891 1.017 1.932 ...
#>  $ T_skin: num  26.4 26.2 25.9 26.7 27 ...
#>  $ P     : num  88.1 88.1 88.1 88.1 88.1 88.1 88.1 88.1 88.2 88.2 ...

# This data is 30-minute data. Therefore timestep = 1800 (1800 seconds in 30 minutes)
# Run aero_calc
evap = suppressWarnings(aero_calc(df=df, sensor_height=2, timestep=1800, out_file_format='csv',
          out_file_name = 'testing', 
          verbose = FALSE))

# View dataframe with new values 
str(evap)
#> 'data.frame':    2928 obs. of  12 variables:
#>  $ date     : POSIXct, format: "2019-09-01 00:00:00" "2019-09-01 00:30:00" ...
#>  $ T_air    : num  28.9 29.2 28.6 28.4 27.8 ...
#>  $ RH       : num  39.5 32.3 38.1 42.1 40.9 ...
#>  $ WS       : num  2.367 2.732 0.891 1.017 1.932 ...
#>  $ T_skin   : num  26.4 26.2 25.9 26.7 27 ...
#>  $ P        : num  88.1 88.1 88.1 88.1 88.1 88.1 88.1 88.1 88.2 88.2 ...
#>  $ SH       : num  2 2 2 2 2 2 2 2 2 2 ...
#>  $ dt       : num  1800 1800 1800 1800 1800 1800 1800 1800 1800 1800 ...
#>  $ E        : num  0.1452 0.1825 0.0659 0.0746 0.0979 ...
#>  $ Ce       : num  0.00218 0.00212 0.00266 0.00259 0.00162 ...
#>  $ VPD      : num  1.86 2.09 1.85 1.87 2.04 ...
#>  $ stability: num  0 0 0 0 0.334 ...

# Plot evaporation rate calculated from aero_calc function
plot(evap$date,evap$E,type='l',col='cadetblue3',
     ylab='Evaporation (mm/30-minute)',
     xlab = 'Date',axes=T)

# Hourly
evap_hourly <- evap %>%
  mutate(hour = floor_date(date, unit = "hour")) %>%
  group_by(hour) %>%
  summarise(across(where(is.numeric), \(x) sum(x, na.rm = TRUE)))

# Daily 
evap_daily <- evap %>%
  mutate(day = floor_date(date, unit = "day")) %>%
  group_by(day) %>%
  summarise(across(where(is.numeric), \(x) sum(x, na.rm = TRUE)))

old_par <- par(no.readonly = TRUE)
on.exit(par(old_par), add = TRUE)
par(mfrow = c(3, 1), mar = c(4, 4, 2, 1))
# Plot 30-minute
plot(evap$date,evap$E,type='l',col='cadetblue3',
     ylab='Evap (mm/30-min)',
     xlab = '',axes=T)
# Plot hourly
plot(evap_hourly$hour,evap_hourly$E,type='l',col='lightpink',
     ylab='Evap (mm/hr)',
     xlab = '',axes=T)
# Plot daily
plot(evap_daily$day,evap_daily$E,type='l',col='darkgreen',
     ylab='Evap (mm/day)',
     xlab = 'Date',axes=T)