#' Prepares weather station data
#' @param WSdata csv file with weather station data or data.frame
#' @param ... additional parameter to pass to read.csv()
#' @param height weather station sensors height in meters
#' @param lat latitude of weather station in decimal degrees.
#' Negative values for south latitude
#' @param long longitude of weather station in decimal degrees.
#' Negative values for west longitude
#' @param elev elevation of weather station in meters
#' @param columns vector with the column numbers in WSdata for the
#' date, time, radiation, wind, RH, temperature and
#' rain. If date and time are in the same column, the
#' column number has to be the same. Names in this
#' vector are ignored and are presented on Usage and
#' examples only as a reference.
#' @param date.format date format. See strptime format argument.
#' @param time.format time format. See strptime format argument.
#' @param datetime.format datetime format. See strptime format argument.
#' @param tz timezone of the weather station dates. If not
#' present assumes the same timezone as the computer running the code. See
#' strptime for details.
#' @param cf conversion factor to convert radiation, wind, and
#' temperature to W/m2; m/s and Celsius. See Details.
#' @param MTL Metadata file. If not provided will look for one on
#' working directory. If provided or present will calculate weather conditions
#' on satellite overpass.
#' @details
#' For cf, if your data is in W/m2, km/h and Celsius (radiation, wind,
#' temperature), cf should be: cf = c(1,0.2777778,1)
#' @examples
#' csvfile <- system.file("extdata", "apples.csv", package="water")
#' MTLfile <- system.file("extdata", "L7.MTL.txt", package="water")
#' WS <- read.WSdata(WSdata = csvfile, date.format = "%d/%m/%Y",
#' lat=-35.42222, long= -71.38639, elev=201, height= 2.2,
#' columns=c("date" = 1, "time" = 2, "radiation" = 3,
#' "wind" = 4, "RH" = 6, "temp" = 7, "rain" = 8),
#' MTL = MTLfile, tz = "America/Santiago")
#' print(WS)
#' plot(WS, alldata=FALSE)
#' plot(WS, alldata=TRUE)
#' @author Guillermo Federico Olmedo
#' @return waterWeatherStation object, with data.frames with all data, hourly
#' data and conditions at satellite flyby.
#' @family Weather station related functions
#' @seealso \code{\link{read.WSdata2}} for the equivalent using read.csv2()
#' @references
#' Landsat 7 Metadata example file available from the U.S. Geological Survey.
#' Weather Station example file courtesy of CITRA, Universidad de Talca, Chile
#' @export
read.WSdata <- function(WSdata, ..., height = 2.2, lat, long, elev,
columns=c("date" = 1, "time" = 2, "radiation" = 3,
"wind" = 4, "RH" = 6, "temp" = 7, "rain" = 8),
date.format = "%Y-%m-%d", time.format = "%H:%M:%S",
datetime.format = "%Y-%m-%d %H:%M:%S", tz = "",
cf = c(1, 1, 1), MTL){
if(tz==""){warning(paste('As tz = "", assuming the weather station time zone is', Sys.timezone()))}
if("temp" %in% columns){warning("The parameter columns has changed and the old
sintaxis is now deprecated. In future versions of water package ONLY the
new (numeric and simple!) sintaxis will work.")}
if("pp" %in% columns){columns[columns == "pp"] <- "rain"} ## TODO: deprecated "pp"
if(class(WSdata) == "character"){WSdata <- utils::read.csv(WSdata, ...)}
else if(class(WSdata) == "data.frame"){WSdata <- WSdata}
else stop("WSdata should be a character string with the name of the csv file
or a data.frame with the WS data")
result <- list()
result$location <- data.frame(lat=lat, long=long, elev=elev, height=height)
### Old method for columns : the character method!
if(class(columns) == "character"){
if("date" %in% columns & "time" %in% columns){
datetime <- paste(WSdata[, which(columns == "date")],
WSdata[, which(columns == "time")])
datetime <- strptime(datetime, format = paste(date.format, time.format),
tz = tz)
} else {
if("datetime" %in% columns){
datetime <- WSdata[, which(columns == "datetime")]
datetime <- strptime(datetime, format = datetime.format, tz = tz)
} else {message("ERROR: date and time or datetime are needed columns")}
}
if("rain" %in% columns){rain = WSdata[, which(columns == "rain")]}
else rain = NA
radiation = WSdata[, which(columns == "radiation")] * cf[1]
wind = WSdata[, which(columns == "wind")] * cf[2]
RH = WSdata[, which(columns == "RH")]
temp = WSdata[, which(columns == "temp")] * cf[3]
}
### New method for columns : the numeric method!
# vector with the column numbers in WSdata for the
# date, time, radiation, wind, RH, temp and rain. If
# date and time are in the same column, the column
# number has to be the same.
if(class(columns)== "numeric"){
if(columns[1] != columns[2]){
datetime <- paste(WSdata[, columns[1]],
WSdata[, columns[2]])
datetime <- strptime(datetime, format = paste(date.format, time.format),
tz = tz)
} else {
datetime <- WSdata[, columns[1]]
datetime <- strptime(datetime, format = datetime.format, tz = tz)
}
if(!is.na(columns[7])){rain = WSdata[, columns[7]]}
else rain = NA
radiation = WSdata[, columns[3]] * cf[1]
wind = WSdata[, columns[4]] * cf[2]
RH = WSdata[, columns[5]]
temp = WSdata[, columns[6]] * cf[3]
}
ea = (RH/100)*0.6108*exp((17.27*temp)/(temp+237.3))
WSdata <- data.frame(datetime=datetime, radiation=radiation, wind=wind,
RH=RH, ea=ea, temp=temp, rain=rain)
result$alldata <- WSdata
## Daily
WSdata$date <- as.Date(WSdata$datetime, tz = tz)
result$daily <- data.frame(date=unique(WSdata$date),
radiation_sum=tapply(WSdata$radiation, WSdata$date, sum),
wind_mean=tapply(WSdata$wind, WSdata$date, mean),
RH_mean=tapply(WSdata$RH, WSdata$date, mean),
temp_mean=tapply(WSdata$temp, WSdata$date, mean),
temp_max=tapply(WSdata$temp, WSdata$date, max),
temp_min=tapply(WSdata$temp, WSdata$date, min),
ea_mean=tapply(WSdata$ea, WSdata$date, mean),
rain_sum=tapply(WSdata$rain, WSdata$date, sum))
## Hourly
result$hourly <- list()
datetime <- as.POSIXlt(WSdata$datetime)
if(datetime[1]$min==0 & datetime[1]$sec==0){
first = datetime[1]
} else {
first <- datetime[1] + 3600 - datetime$min[1] * 60 - datetime$sec[1]
}
sequence <- seq.POSIXt(from=first, to = tail(datetime,1),by= "1 hour")
# Time interpolation
radiationFun <- stats::approxfun(result$alldata$datetime, result$alldata$radiation)
windFun <- stats::approxfun(result$alldata$datetime, result$alldata$wind)
RHFun <- stats::approxfun(result$alldata$datetime, result$alldata$RH)
eaFun <- stats::approxfun(result$alldata$datetime, result$alldata$ea)
tempFun <- stats::approxfun(result$alldata$datetime, result$alldata$temp)
if(!all(is.na(rain))){
rainFun <- stats::approxfun(result$alldata$datetime, result$alldata$rain)
} else { rainFun <- function(seq){seq <- NA}}
result$hourly <- data.frame(datetime = sequence,
radiation = radiationFun(sequence),
wind = windFun(sequence),
RH = RHFun(sequence),
ea = eaFun(sequence),
temp = tempFun(sequence),
rain = rainFun(sequence))
## Join with satellite data
if(missing(MTL)){MTL <- list.files(pattern = "MTL.txt", full.names = T)}
if(length(MTL)!=0){
MTL <- readLines(MTL, warn=FALSE)
time.line <- grep("SCENE_CENTER_TIME",MTL,value=TRUE)
date.line <- grep("DATE_ACQUIRED",MTL,value=TRUE)
sat.time <-regmatches(time.line,regexec(text=time.line,
pattern="([0-9]{2})(:)([0-9]{2})(:)([0-9]{2})(.)([0-9]{2})"))[[1]][1]
sat.date <-regmatches(date.line,regexec(text=date.line,
pattern="([0-9]{4})(-)([0-9]{2})(-)([0-9]{2})"))[[1]][1]
sat.datetime <- strptime(paste(sat.date, sat.time),
format = "%Y-%m-%d %H:%M:%S", tz="Europe/London")
WS.prev<-suppressWarnings(WSdata[WSdata$datetime == tail(datetime[datetime <
sat.datetime],1),])
WS.after <- suppressWarnings(WSdata[WSdata$datetime == datetime[datetime >
sat.datetime][1],])
delta1 <- as.numeric(difftime(WS.after$datetime,
WS.prev$datetime, units="secs"))
delta2 <- as.numeric(difftime(sat.datetime, WS.prev$datetime, units="secs"))
sat.data <- WS.prev + (WS.after - WS.prev)/delta1 * delta2
sat.data[,2:6] <- round(sat.data[,2:6],2)
result$at.sat <- sat.data
}
class(result) <- "waterWeatherStation"
return(result)
}
#' Prepares weather station data 2
#' @param WSdata csv file with weather station data
#' @param ... additional parameter to pass to read.csv()
#' @param height weather station sensors height in meters
#' @param lat latitude of weather station in decimal degrees.
#' Negative values for south latitude
#' @param long longitude of weather station in decimal degrees.
#' Negative values for west longitude
#' @param elev elevation of weather station in meters
#' @param columns columns order of needed data. Vector containing
#' "date", "time", "radiation", "wind", "RH", "temp" and "rain". Other values are
#' ignored. If you have a column with date and time in the same column, you can
#' include "datetime" and "date" and "time" are no longer needed.
#' @param date.format date format. See strptime format argument.
#' @param time.format time format. See strptime format argument.
#' @param datetime.format datetime format. See strptime format argument.
#' @param tz timezone of the weather station dates. If not
#' present assumes the same timezone as the computer running the code. See
#' strptime for details.
#' @param cf conversion factor to convert radiation, wind, and
#' temperature to W/m2; m/s and Celsius. See Details.
#' @param MTL Metadata file. If not provided will look for one on
#' working directory. If provided or present will calculate weather conditions
#' on satellite flyby.
#' @details
#' For cf, if your data is in W/m2, km/h and Celsius (radiation, wind,
#' temperature), cf should be: cf = c(1,0.2777778,1)
#' @family Weather station related functions
#' @author Guillermo Federico Olmedo
#' @export
read.WSdata2 <- function(WSdata, ..., height = 2.2, lat, long, elev,
columns=c("date" = 1, "time" = 2, "radiation" = 3,
"wind" = 4, "RH" = 6, "temp" = 7, "rain" = 8),
date.format = "%d/%m/%Y", time.format = "%H:%M:%S",
datetime.format = "%Y-%m-%d %H:%M:%S", tz = "",
cf = c(1, 3.6, 1, 1), MTL){
read.WSdata(WSdata, ..., height = 2.2, lat = lat, long = long, elev = elev,
columns = columns, date.format = date.format,
time.format = time.format, datetime.format = datetime.format,
tz = tz, cf = cf, MTL)
}
#' Plot method for waterWeatherStation S3 class
#' @param x waterWeatherStation object. See read.WSdata()
#' @param hourly If TRUE will plot only hourly data, instead of all data
#' from the waterWeatherStation object
#' @param sat If TRUE, and if the waterWeatherStation object was created
#' using a Landsat Metadata File, will plot only data from day of the satellite
#' overpass
#' @param date When sat = FALSE, the date to plot in a multi-day object
#' @param ... additional parameters to pass to plot()
#' @author Guillermo Federico Olmedo
#' @importFrom utils read.csv
#' @importFrom graphics abline axis axis.POSIXct mtext par points
#' @export
#' @family Weather station related functions
#' @method plot waterWeatherStation
plot.waterWeatherStation <- function(x, hourly=FALSE, sat=TRUE, date, ...){
# Based on http://evolvingspaces.blogspot.cl/2011/05/multiple-y-axis-in-r-plot.html
WSp <- x$hourly
if(hourly == FALSE) {WSp <- x$alldata}
long=FALSE
daily_mean = ""
if(sat == TRUE & "datetime" %in% names(x$at.sat)){
atsat <- as.POSIXlt(x$at.sat$datetime)
WSp <- WSp[as.POSIXlt(WSp$datetime)$yday == atsat$yday
& as.POSIXlt(WSp$datetime)$year == atsat$year,]
} else {
if(nrow(x$daily)>5){WSp <- x$daily
names(WSp) <- c("datetime", "radiation", "wind", "RH", "temp",
"temp_max", "temp_min", "ea", "rain")
long=TRUE
daily_mean = "(mean daily value)"}
}
if(!missing(date)){
date <- as.POSIXlt(date)
WSp <- WSp[as.POSIXlt(WSp$datetime)$yday == date$yday
& as.POSIXlt(WSp$datetime)$year == date$year,]
}
WSp <- WSp[order(WSp$datetime),]
time <- WSp$datetime
graphics::par(mar=c(5, 7, 4, 9.5) + 0.1)
plot(time, WSp$radiation, axes=F, ylim=c(0,max(WSp$radiation)), xlab="",
ylab="",type="l",col="red", xlim=range(time), ...)
if(sat == TRUE & "datetime" %in% names(x$at.sat)){
graphics::abline(v=as.POSIXct(atsat), lwd=5, col="gray")
graphics::text(atsat, max(WSp$radiation)*0.85, "satellite overpass", cex=0.7,
adj=c(NA, -0.5), srt=90, col="gray")
}
graphics::points(time,WSp$radiation,pch=20,col="red")
graphics::axis(2, ylim=c(0,max(WSp$radiation)),col="red",lwd=1, cex.axis=0.5, tcl=-0.25)
graphics::mtext(2,text=paste("Solar radiation (W.m-2)", daily_mean),line=1.7, cex=0.7)
# Rain
#WSp[WSp$rain == 0] <- NA
if(!all(is.na(WSp$rain))){
graphics::par(new=T)
plot(time, WSp$rain, axes=F, ylim=c(0,max(WSp$rain)*1.7), xlab="", ylab="",
type="h",col="light blue",lty=1, main="",xlim=range(time),lwd=3)
graphics::axis(4, ylim=c(0,max(WSp$rain)*1.7),col="light blue", line = 6.7, lwd=1, cex.axis=0.5, tcl=-0.25)
graphics::mtext(4,text=paste("Rain (mm)", daily_mean), cex=0.7, line = 8.4)
}
#
graphics::par(new=T)
plot(time, WSp$wind, axes=F, ylim=c(0,max(WSp$wind)), xlab="", ylab="",
type="l",col="green",lty=2, main="",xlim=range(time),lwd=2)
graphics::axis(2, ylim=c(0,max(WSp$wind)),col="green",lwd=1,line=3.5, cex.axis=0.5, tcl=-0.25)
graphics::points(time, WSp$wind,pch=20,col="green")
graphics::mtext(2,text=paste("Wind speed (m.s-1)", daily_mean),line=5.5, cex=0.7)
graphics::par(new=T)
plot(time, WSp$temp, axes=F, ylim=c(0,max(WSp$temp)), xlab="", ylab="",
type="l",col="black",lty=2, main="",xlim=range(time),lwd=2)
graphics::axis(4, ylim=c(0,max(WSp$temp)),col="black",lwd=1, cex.axis=0.5, tcl=-0.25)
graphics::points(time, WSp$temp,pch=20,col="black")
graphics::mtext(4,text=paste("Temperature (C)", daily_mean),line=1.7, cex=0.7)
graphics::par(new=T)
plot(time, WSp$ea, axes=F, ylim=c(0,max(WSp$ea)), xlab="", ylab="",
type="l",col="blue",lty=2, main="",xlim=range(time),lwd=2)
graphics::axis(4, ylim=c(0,max(WSp$ea)),col="blue",lwd=1,line=3.5, cex.axis=0.5, tcl=-0.25)
graphics::points(time, WSp$ea,pch=20,col="blue")
graphics::mtext(4,text=paste("vapor pressure (kPa)", daily_mean),line=5.5, cex=0.7)
if(long==FALSE){r <- as.POSIXct(round(range(time), "hours"))
graphics::axis.POSIXct(1, at = seq(r[1], r[2], by = "hour"), format = "%H:%M")}
if(long==TRUE){r <- range(time)
graphics::axis.Date(1, at = seq(r[1], r[2], by = "day"))}
graphics::mtext("Time",side=1,col="black",line=2)
}
#' Print method for waterWeatherStation S3 class
#' @param x waterWeatherStation object. See read.WSdata()
#' @param ... additional parameters to pass to print()
#' @author Guillermo Federico Olmedo
#' @export
#' @family Weather station related functions
#' @method print waterWeatherStation
print.waterWeatherStation <- function(x, ...){
cat("Weather Station @ lat:", round(x$location$lat, 2), "long:",
round(x$location$long, 2), "elev:", round(x$location$elev, 2), "\n")
if(!is.null(x$at.sat)){
cat("Summary:\n")
print(summary(x$alldata[format(x$alldata$datetime, format = "%d-%m-%Y") ==
format(x$at.sat$datetime, format = "%d-%m-%Y"), 2:7]))
cat("\n Conditions at satellite overpass:\n")
print(x$at.sat)} else {
cat("Summary:\n")
print(summary(x$alldata[,1:7]))
}
}
#' Export data.frame from waterWeatherStation Object
#' @description
#' Export weather conditions at satellite flyby from waterWeatherStation object
#' @param WeatherStation waterWeatherStation object. See read.WSdata()
#' @author Guillermo Federico Olmedo
#' @export
getDataWS <- function(WeatherStation){
date <- as.POSIXlt(WeatherStation$at.sat$datetime, format="%Y-%m-%d %H:%M:%S")
WeatherStation2 <- data.frame(wind=WeatherStation$at.sat$wind,
RH=WeatherStation$at.sat$RH,
temp=WeatherStation$at.sat$temp,
radiation=WeatherStation$at.sat$radiation,
height=WeatherStation$location$height,
lat=WeatherStation$location$lat,
long=WeatherStation$location$long,
elev=WeatherStation$location$elev,
DOY=date$yday+1, hours=date$hour + date$min/60 +
date$sec/3600)
return(WeatherStation2)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.