derived.data <- function(lon, lat, cops.init, cops.raw) {
mymessage(paste(" ", "calculating derived data"))
if (str_detect(cops.raw$file, "_SB_")) {
print("Shadow band file format detected")
print("Change the init parameters")
number.of.fields.before.date.tmp <- number.of.fields.before.date
number.of.fields.before.date <- number.of.fields.before.date - 1
instruments.optics.tmp <- instruments.optics
instruments.optics <- "Ed0"
}
instruments <- c(instruments.optics)
d2r <- pi / 180
# tilt from Roll and Pitch
ret <- list()
for(instr in instruments) {
instr.anc <- paste(instr, "anc", sep = ".")
if("Roll" %in% names(cops.raw[[instr.anc]])) {
Roll <- cops.raw[[instr.anc]]$Roll
Pitch <- cops.raw[[instr.anc]]$Pitch
tilt <- atan(sqrt(tan(Roll*d2r)^2+tan(Pitch*d2r)^2))/d2r
} else {
tilt <- NULL
}
ret <- c(ret, list(tilt))
}
names(ret) <- paste(instruments, "tilt", sep = ".")
# dates
##### Added by Simon Belanger, 31 march 2016
if (str_detect(cops.raw$Others$DateTime[1], "PM") | str_detect(cops.raw$Others$DateTime[1], "AM")) {
cops.raw$Others$DateTime = convert.time.12h.to.24h.string(cops.raw$Others$DateTime, cops.init, 0)
}
####
dates <- as.POSIXct(strptime(cops.raw$Others$DateTime, format = format.date[1]))
dates.secs.from.beginning <- as.numeric(dates) - min(as.numeric(dates))
dates.good <- dates.secs.from.beginning >= time.window[1] & dates.secs.from.beginning <= time.window[2]
cops.duration.secs <- max(as.numeric(dates[dates.good])) - min(as.numeric(dates[dates.good]))
dte <- mean.POSIXct(dates[dates.good])
# dte <- dates[dates.good][1]
print(dates[1])
# longitude and latitude if BioGPS.Position column is present
if("BioGPS.anc" %in% names(cops.raw)) {
dummy <- cops.raw[["BioGPS.anc"]]$Position
position.signs <- sign(dummy)
dummy <- abs(dummy)
pos.deg <- floor(dummy / 100)
pos.min <- dummy - pos.deg * 100
position <- (pos.deg + pos.min / 60) * position.signs
i.lon <- abs(position - lon) < 0.5
i.lat <- abs(position - lat) < 0.5
if(length(which(i.lon)) < 0.5 | length(which(i.lat)) < 0.5) {
cat("if BioGPS.Position column is present in data file, you have to indicate\n")
cat(" approximate (+- 0.5 degree) longitude and latitude in file info.cops.dat\n")
cat("mismatch between longitude and latitude in file info.cops.dat\n")
cat(" ", lon, lat, "\n")
cat("AND longitude and latitude in data file (BioGPS.Position column)\n")
hh <- hist(position, plot = FALSE)
cat(" here is the histogram of the values\n")
cat(" breaks :", hh$breaks, "\n")
cat(" counts :", hh$counts, "\n")
stop()
}
longitude <- median(position[i.lon])
latitude <- median(position[i.lat])
change.position <- TRUE
} else {
# if not
if(is.na(lon) | is.na(lat)) {
# look for a gps file
gps.files <- paste(dirdat, cops.raw$potential.gps.file, sep = "/")
if(! file.exists(gps.files)) {
cat("longitude and(or) latitude not available\n")
cat(cops.raw$potential.gps.files, "\n")
cat("3 possibilities :\n")
cat("- PUT longitude and latitude in file info.cops.dat : 2nd and 3rd column\n")
cat("- USE BioGPS.Position column if present in data file; in this case,\n")
cat(" you have to indicate in file info.cops.dat (2nd and 3rd column)\n")
cat(" approximate (+- 0.5 degree) longitude and latitude\n")
cat("- PUT NA in file info.cops.dat (2nd and 3rd column)iin this case,\n")
cat(" a valid GPS file must be present in directory", dirdat,"\n")
cat(" such a file has a name with the same first fields as the name of the present data file\n")
cat(" until the time field, then add \"_gps\" and any extension\n")
stop()
} else {
# Added by Simon Belanger
if (length(gps.files) > 1) {
cat("Multiple gps files available\n")
xx = unlist(strsplit(cops.raw$file, "_data_"))
yy = unlist(strsplit(gps.files, "_gps_"))
ix.file = which(yy == xx[2]) / 2
gps.file=gps.files[ix.file]
} else gps.file=gps.files
#Added by Servet Cizmeli and Simon Belanger
print("Reading GPS file: ")
print(gps.file)
# Modified by Simon Belanger on Aug 2016 to process
# GPS data obtained with uprofile 1.9.10 and after
if (str_detect(gps.file, "GPS_")) {
ext = unlist(strsplit(gps.file, "[.]"))[2]
if (ext == "tsv" || ext =="txt") {
gps.data <- read.table(gps.file,sep="\t",
colClasses = c("character","character","numeric","character","numeric","numeric","numeric"),header=TRUE)
} else {
gps.data <- read.table(gps.file,sep=",",
colClasses = c("character","character","numeric","character","numeric","numeric","numeric"),header=TRUE)
}
names(gps.data)[4] <- "GpsTime"
names(gps.data)[1] <- "ComputerTime"
} else {
ext = unlist(strsplit(gps.file, "[.]"))[2]
if (ext == "tsv" || ext =="txt") {
gps.data <- read.table(gps.file,sep="\t",
colClasses = c("character","numeric","character","numeric","numeric","numeric"),header=TRUE)
} else {
gps.data <- read.table(gps.file,sep=",",
colClasses = c("character","numeric","character","numeric","numeric","numeric"),header=TRUE)
}
names(gps.data)[3] <- "GpsTime"
names(gps.data)[1] <- "ComputerTime"
}
### END
if (str_detect(gps.data$ComputerTime[1], "PM") | str_detect(gps.data$ComputerTime[1], "AM")) {
gps.data$ComputerTime = convert.time.12h.to.24h.string(gps.data$ComputerTime, cops.init, 0)
}
if (str_detect(gps.data$GpsTime[1], "PM") | str_detect(gps.data$GpsTime[1], "AM")) {
gps.data$GpsTime = convert.time.12h.to.24h.string(gps.data$GpsTime, cops.init, 0)
}
### Add by S Belanger in 2018 to deal with GPS having different
### DateTime format than the COPS profiles (in was the case in 2015).
### The user need to provide two format.date in the init.cops.dat file
### separated by coma.
if (length(cops.init$format.date)==2) {
gps.data$GpsTime <- as.POSIXct(strptime(gps.data$GpsTime, format = cops.init$format.date[2]))
gps.data$ComputerTime <- as.POSIXct(strptime(gps.data$ComputerTime, format = cops.init$format.date[2]))
} else {
gps.data$GpsTime <- as.POSIXct(strptime(gps.data$GpsTime, format = cops.init$format.date[1]))
gps.data$ComputerTime <- as.POSIXct(strptime(gps.data$ComputerTime, format = cops.init$format.date[1]))
}
if (is.na(gps.data$GpsTime[1])) {
print("WARNING: It seems that the GPS time format differs from the CAST files")
print("1. Check the time format inside the GPS file")
print("2. If it differs from the format of the DateTime column of the CAST files,")
print(" than add a second time format in the init.cops.dat at the ")
print(" format.date field (separated by coma).")
print(" ")
print(" Teminate processing of that profile.")
return(0)
}
valid_gps <- gps.data$GpsTime > min(dates) & gps.data$GpsTime < max(dates)
if (any(valid_gps)) {
longitude <- median(gps.data$Longitude[valid_gps], na.rm = TRUE)
latitude <- median(gps.data$Latitude[valid_gps], na.rm = TRUE)
change.position <- TRUE
} else {
valid_gps <- gps.data$ComputerTime > min(dates) & gps.data$ComputerTime < max(dates)
if (any(valid_gps)) {
longitude <- median(gps.data$Longitude[valid_gps], na.rm = TRUE)
latitude <- median(gps.data$Latitude[valid_gps], na.rm = TRUE)
change.position <- TRUE
} else {
print("Time found in gps file do not match time in")
print(cops.raw$file)
}
}
}
#END - Added by Servet Cizmeli
} else {
# keep longitude and latitude of file info.cops.dat
longitude <- lon
latitude <- lat
change.position <- FALSE
}
}
cat("longitude :", longitude, "\n")
cat("latitude :", latitude, "\n")
print(paste("date: ", dte))
# sun-zenith angle
day <- as.numeric(format(dte, format = "%d"))
month <- as.numeric(format(dte, format = "%m"))
year <- as.numeric(format(dte, format = "%Y"))
hour <- as.numeric(format(dte, format = "%H"))
minute <- as.numeric(format(dte, format = "%M"))
second <- as.numeric(format(dte, format = "%S"))
ah <- hour + minute / 60 + second / 3600
print(month)
sunzen <- possol(month, day, ah, longitude, latitude)[1]
if (str_detect(cops.raw$file, "_SB_")) {
Depth<-NULL
Depth.good<-NULL
depth.fitted <- NULL
ret <- c(ret, list(change.position = change.position,
longitude = longitude,
latitude = latitude,
dates = dates,
date.mean = dte,
cops.duration.secs = cops.duration.secs,
day = day,
month = month,
year = year,
sunzen = sunzen,
Depth = Depth,
Depth.good = Depth.good,
depth.fitted = depth.fitted))
} else {
# ELIMINATE A FEW INVALID DEPTHS
Depth <- cops.raw[[paste(depth.is.on, "anc", sep = ".")]]$Depth #+ delta.capteur.optics["LuZ"]
# median filter on Depth with big delta to remove a few inappropriate depths
delta <- max(Depth) / length(Depth) * 50
Depth.filtered <- filtre.mediane(3, Depth, delta = delta, fill = TRUE, replace = FALSE)
Depth.good <- !is.na(Depth.filtered)
# agglomerate dates.good with Depth.good
Depth.good <- Depth.good & dates.good
# depth where data will be fitted
depth.fitted <- NULL
for(i in seq(1, length(depth.discretization) - 2, 2)) {
depth.fitted <- append(depth.fitted, seq(depth.discretization[i], depth.discretization[i + 2] - depth.discretization[i + 1] / 2, depth.discretization[i + 1]))
}
maxdepth <- max(Depth[Depth.good]+ delta.capteur.optics["LuZ"])
depth.fitted <- depth.fitted[depth.fitted <= maxdepth]
change.interval <- exists("time.interval.for.smoothing.optics")
if (change.interval) {
depth.interval.for.smoothing.optics <- time.interval.for.smoothing.optics/cops.duration.secs* maxdepth
}
ret <- c(ret, list(change.position = change.position,
longitude = longitude,
latitude = latitude,
dates = dates,
date.mean = dte,
cops.duration.secs = cops.duration.secs,
day = day,
month = month,
year = year,
sunzen = sunzen,
Depth = Depth,
Depth.good = Depth.good,
depth.fitted = depth.fitted))
# PLOT
if(INTERACTIVE) x11(width = win.width, height = win.height)
plot(dates, Depth, ylim = rev(range(Depth)), xlab = "UTC Time (HH:MM:SS)", ylab = "Depth (m)", main = paste(dte, " duration =", cops.duration.secs, "s", " lon =", round(longitude, digits = 3), " lat =", round(latitude, digits = 3), " sun-zenith angle =", round(sunzen, digits = 2), "deg."), cex.main = 1, axes = FALSE, frame.plot = TRUE)
grid(col = 1)
axis.POSIXct(1, dates, format = "%H:%M:%S")
axis(2)
points(dates[!Depth.good], Depth[!Depth.good], pch = 20, col = 2, cex = 2)
legend("bottomleft", legend = paste(length(which(!Depth.good)), "points removed"), text.col = 2, cex = 2)
if(change.position) {
if (change.interval) {
legend("topright", legend = c("Longitude and Latitute found",
" (BioGPS column present) and",
" time.interval.for.smoothing",
"VALUES MODIFIED for depth.interval in meters",
signif(depth.interval.for.smoothing.optics,3)),
text.col = 3, cex = 1)
} else legend("topright", legend = c("Longitude and Latitute found",
" (BioGPS column present)",
"VALUES MODIFIED"),
text.col = 3, cex = 1)
} else {
if (change.interval) legend("topright", legend = c("time.interval.for.smoothing",
"VALUES MODIFIED for depth.interval in meters",
signif(depth.interval.for.smoothing.optics,3)),
text.col = 3, cex = 1)
}
if(INTERACTIVE) x11(width = win.width, height = win.height)
par(mfrow = c(1, length(instruments)))
for(instr in instruments) {
tilt <- ret[[paste(instr, "tilt", sep = ".")]]
tiltmax <- tiltmax.optics[instr]
if(!is.null(tilt)) {
plot(tilt[Depth.good], Depth[Depth.good],
ylim = rev(range(Depth)),
xlab = "Tilt (degree)", ylab = "Depth (m)", main = instr)
abline(v=tiltmax.optics[instr], col=2, lwd=2)
grid(col = 1)
points(tilt[!Depth.good], Depth[!Depth.good], col = 8) # grey
points(tilt[tilt > tiltmax], Depth[tilt > tiltmax], col = 2) # red
}
}
par(mfrow = c(1, 1))
}
ret
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.