Nothing
R/data_preparation.R
In FLightR: Reconstruct Animal Paths from Solar Geolocation Loggers Data
Defines functions
geologger.sampler.create.arrays
create.proposal
make.processed.light.object
make.result.list
correct.hours
create.calibration
make_likelihood_correction_function
suggest.irrad.boundaries
get.calibration.parameters
plot_slopes
get.calib.param
logger.template.calibration
logger.template.calibrarion.internal
get.Irradiance
process.twilights
make.prerun.object
find.stationary.location
make.calibration
plot_slopes_by_location
Documented in
find.stationary.location
make.calibration
make.prerun.object
plot_slopes_by_location
# data_preparation functions
#' plots log of observed versus expected slope by time for a known location
#'
#' The function calculates and plots calibration slopes for sunsets and sunrises for every day of the tracking period, based on the assumption that the tag remained in the same (calibration) location all the time.
#'
#' @param Proc.data processed data object generated by \code{\link{get.tags.data}}
#' @param location vector with longitude and latitude of calibration location (degrees).
#' @param log.light.borders numeric vector with length of 2 for minimum and maximum log(light) levels to use. Default value 'auto', will take these values from the Proc.data object.
#' @param log.irrad.borders numeric vector with length of 2 for minimum and maximum log(irradiance) values to use. Default value 'auto', will take these values from the Proc.data object.
#' @param ylim the y limits of the plot. The default value, NULL, indicates that the range of the finite values to be plotted should be used.
#' @param xlim the x limits of the plot. The default value, NULL, otherwise can be POSIXct or character in a form readable by \code{\link{as.POSIXct}}.
#' @details The plot of calibration slopes is used for finding start and end dates of a calibration period (the time period, during which the tag remained in the calibration location with coordinates (x,y)). During the calibration period, the calibration slopes vary little both, between the twilight events (sunrises and sunsets) and in time. When the tag changes location, the slopes for sunrises and sunsets start to deviate. There may potentially be several calibration periods for the same location (if the bird returned to the same location several times). The boundaries (start and end dates) of each of these periods are captured visually. If there were more than one calibration location, the procedure is repeated, once for each location.
#' All the obtained calibration periods can be entered in a data frame 'Calibration.periods', for further analysis. Each line of the data frame contains start and end dates (if applicable) of the calibration period and geographic coordinates of the location.
#' @return 'NULL'
#' @examples
#' File<-system.file("extdata", "Godwit_TAGS_format.csv", package = "FLightR")
#' Proc.data<-get.tags.data(File)
#' plot_slopes_by_location(Proc.data=Proc.data, location=c(5.43, 52.93))
#' abline(v=as.POSIXct("2013-08-20", tz='GMT')) # end of first calibration period
#' abline(v=as.POSIXct("2014-05-05", tz='GMT')) # start of the second calibration period
#'
#' @author Eldar Rakhimberdiev
#' @export plot_slopes_by_location
plot_slopes_by_location<-function(Proc.data, location, log.light.borders='auto', log.irrad.borders='auto', ylim=NULL, xlim=NULL) {
old.par <- graphics::par(no.readonly = TRUE)
Calibration.period<-data.frame(
calibration.start=as.POSIXct("1900-01-01", tz='GMT'),
calibration.stop=as.POSIXct("2050-01-01", tz='GMT'),
lon=location[1], lat=location[2])
if (log.light.borders[1]=='auto') log.light.borders<-Proc.data$log.light.borders
if (log.irrad.borders[1]=='auto') log.irrad.borders<-Proc.data$log.irrad.borders
calibration.parameters<-suppressWarnings(get.calibration.parameters(Calibration.period,
Proc.data, model.ageing=FALSE,
log.light.borders=log.light.borders,
log.irrad.borders=log.irrad.borders,
plot.each = FALSE, plot.final = FALSE, suggest.irrad.borders=FALSE))
tmp<-suppressWarnings(plot_slopes(calibration.parameters$All.slopes, ylim=ylim, xlim=xlim))
#graphics::par(old.par)
return(tmp)
}
#' Creates a calibration object, further used for calculation of coordinates in the \code{\link{run.particle.filter}}.
#'
#' Function estimates all necessary parameters from the calibration data logged in a known location or locations.
#'
#' @param Proc.data processed data object generated by \code{\link{get.tags.data}}
#' @param Calibration.periods a data frame containing start and end dates of all the calibration periods (POSIXct) and geographic coordinates of the corresponding calibration locations.
#' @param model.ageing if set to TRUE, accounts for the tag ageing (with opacification of its transparent shell of a light sensor), resulting into decreasing sensitivity of the device. This option is useful only if there were several calibration periods or if calibration period was very long (~ longer than a month).
#' @param plot.each Do you want every twilight to be plotted while processing
#' @param plot.final Do you want final calibration graph to be plotted. On the graph you can see all the observed versus expected light levels. All slopes should be similar.
#' @param likelihood.correction will estimate correction of likelihood for the current calibration parameters. Highly recommended not to be change from 'auto'. In this case FLightR will switch it to FALSE in case tag saved data on 10 minutes or longer period.
#' @param fixed.logSlope these are mean (1) and SD (2) for distribution of slopes. Should normally be estimated from the data (and thus default is c(NA, NA)). Change any of these two finite values if you want them to be predetermined and not estimated from the calibration data.
#' @param suggest.irrad.borders experimental parameter! If set to TRUE function will try to find the best values for the log.irrad.borders
#' @param return.slopes if true function will return estimated individual twilight slopes.
#' @return calibration object to be uses in the \code{\link{make.prerun.object}}
#' @examples
#' File<-system.file("extdata", "Godwit_TAGS_format.csv", package = "FLightR")
#' Proc.data<-get.tags.data(File, end.date=as.POSIXct('2013-08-20', tz='GMT'))
#' Calibration.periods<-data.frame(
#' calibration.start=NA,
#' calibration.stop=as.POSIXct("2013-08-20", tz='GMT'),
#' lon=5.43, lat=52.93)
#' #use c() also for the geographic coordinates, if you have more than one calibration location
#' # (e. g., lon=c(5.43, 6.00), lat=c(52.93,52.94))
#' print(Calibration.periods)
#'
#' # NB Below likelihood.correction is set to FALSE for fast run!
#' # Leave it as default 'auto' for real examples
#' Calibration<-make.calibration(Proc.data, Calibration.periods, likelihood.correction=FALSE)
#'
#' @author Eldar Rakhimberdiev
#' @export
make.calibration<-function(Proc.data, Calibration.periods, model.ageing=FALSE, plot.each=FALSE, plot.final=FALSE, likelihood.correction='auto', fixed.logSlope=c(NA,NA), suggest.irrad.borders=FALSE, return.slopes=FALSE) {
Calibration.periods[,1]<-as.POSIXct(Calibration.periods[,1], tz='GMT')
Calibration.periods[,2]<-as.POSIXct(Calibration.periods[,2], tz='GMT')
Calibration.periods$calibration.start[is.na(Calibration.periods$calibration.start)]<-as.POSIXct("1900-01-01", tz='GMT')
Calibration.periods$calibration.stop[is.na(Calibration.periods$calibration.stop)]<-as.POSIXct("2100-01-01", tz='GMT')
for (i in 1:nrow(Calibration.periods)) {
if (Calibration.periods[i,1]>=Calibration.periods[i,2]) stop('Calibration start in some period is later than calibration end')
if (i>1) if ( Calibration.periods[i,1]<=Calibration.periods[(i-1),2]) stop('Calibration periods overlap')
}
calibration.parameters<-suppressWarnings(get.calibration.parameters(
Calibration.periods, Proc.data,
model.ageing=model.ageing,
log.light.borders=Proc.data$log.light.borders,
log.irrad.borders=Proc.data$log.irrad.borders,
plot.each= plot.each, plot.final= plot.final,
suggest.irrad.borders=suggest.irrad.borders))
# here I check whether tag measures with 10 minutes or more inteval and if so change likelihood correction to FALSE
if (likelihood.correction=='auto') {
if (Proc.data$saving.period>=550) {
likelihood.correction<-FALSE
message('likelihood correction switched to FALSE\n')
} else {
likelihood.correction<-TRUE
message('likelihood correction switched to TRUE\n')
}
}
Proc.data$log.irrad.borders=calibration.parameters$log.irrad.borders
if (length(calibration.parameters$calib_outliers)>0) {
Proc.data$FLightR.data$twilights$excluded[which(sapply(Proc.data$FLightR.data$twilights$datetime,
FUN=function(x) min(abs(calibration.parameters$calib_outliers-as.numeric(x))))<Proc.data$saving.period*24)]<-1
Proc.data_tmp<-process.twilights(Proc.data$FLightR.data$Data,
Proc.data$FLightR.data$twilights[Proc.data$FLightR.data$twilights$excluded==0,],
measurement.period=Proc.data$measurement.period, saving.period=Proc.data$saving.period,
impute.on.boundaries=Proc.data$impute.on.boundaries)
calibration.parameters<-suppressWarnings(get.calibration.parameters(Calibration.periods, Proc.data_tmp,
model.ageing=model.ageing, log.light.borders=Proc.data$log.light.borders,
log.irrad.borders=Proc.data$log.irrad.borders,
suggest.irrad.borders=suggest.irrad.borders))
plot_slopes(calibration.parameters$All.slopes)
Proc.data$Twilight.time.mat.dusk<-Proc.data_tmp$Twilight.time.mat.dusk
Proc.data$Twilight.time.mat.dawn<-Proc.data_tmp$Twilight.time.mat.dawn
Proc.data$Twilight.log.light.mat.dusk<-Proc.data_tmp$Twilight.log.light.mat.dusk
Proc.data$Twilight.log.light.mat.dawn<-Proc.data_tmp$Twilight.log.light.mat.dawn
}
# second time search for outliers..
if (length(calibration.parameters$calib_outliers)>0) {
Proc.data$FLightR.data$twilights$excluded[which(sapply(Proc.data$FLightR.data$twilights$datetime,
FUN=function(x) min(abs(calibration.parameters$calib_outliers-as.numeric(x))))<Proc.data$saving.period*24)]<-1
Proc.data_tmp<-process.twilights(Proc.data$FLightR.data$Data,
Proc.data$FLightR.data$twilights[Proc.data$FLightR.data$twilights$excluded==0,],
measurement.period=Proc.data$measurement.period, saving.period=Proc.data$saving.period,
impute.on.boundaries=Proc.data$impute.on.boundaries)
calibration.parameters<-suppressWarnings(get.calibration.parameters(
Calibration.periods, Proc.data_tmp,
model.ageing=model.ageing, log.light.borders=Proc.data$log.light.borders,
log.irrad.borders=Proc.data$log.irrad.borders,
suggest.irrad.borders=suggest.irrad.borders))
plot_slopes(calibration.parameters$All.slopes)
Proc.data$Twilight.time.mat.dusk<-Proc.data_tmp$Twilight.time.mat.dusk
Proc.data$Twilight.time.mat.dawn<-Proc.data_tmp$Twilight.time.mat.dawn
Proc.data$Twilight.log.light.mat.dusk<-Proc.data_tmp$Twilight.log.light.mat.dusk
Proc.data$Twilight.log.light.mat.dawn<-Proc.data_tmp$Twilight.log.light.mat.dawn
}
if (!is.na(fixed.logSlope[1])) calibration.parameters$All.slopes$Parameters$LogSlope[1]<-fixed.logSlope[1]
if (!is.na(fixed.logSlope[2])) calibration.parameters$All.slopes$Parameters$LogSlope[2]<-fixed.logSlope[2]
# clean the LogSigma values not to use it under the fixed slope scenario
if (!is.na(fixed.logSlope[1])) calibration.parameters$All.slopes$Parameters$LogSigma<-c(NA,NA)
Calibration<-create.calibration(calibration.parameters$All.slopes,
Proc.data,
Proc.data$FLightR.data,
log.light.borders=Proc.data$log.light.borders, log.irrad.borders=Proc.data$log.irrad.borders,
ageing.model=calibration.parameters$ageing.model,
location=NA, likelihood.correction=likelihood.correction)
Calibration$Calibration.periods<-Calibration.periods
if (return.slopes) Calibration$Calibration.Slopes<-calibration.parameters$All.slopes$Slopes
return(Calibration)
}
#' find unknown calibration location
#'
#' Functions attempts to find a location where
#' The function attempts to find a location for a time period assuming animal was not moving. Does not work well will shaded data!
#'
#' @param Proc.data processed data object generated by \code{\link{get.tags.data}}
#' @param calibration.start POSIXct time when stationary period started
#' @param calibration.stop POSIXct time when stationary period ended
#' @param plot plots every iteration
#' @param initial.coords location vector with initial values for location (longitude and latitude). Should be close (+-2000 km from the real location)
#' @param print.optimization do you want every optimization iteration to be printed? If TRUE - Lon, Lat, calibration mean and calibration sd are being printed. Optimization tries to minimize the latter.
#' @param reltol tolerance for optimization, see \code{\link{optim}} for more details
#' @details The idea behind the function is that it tries to minimize variance between slopes for the whole period by optimizing location. It can be seen as an extension of Hill-Ekstrom calibration idea.
#' @return vector with coordinates - longitude and latitude.
#' @examples
#' #this example takes about 15 minutes to run
#' \dontrun{
#' File<-system.file("extdata", "Godwit_TAGS_format.csv", package = "FLightR")
#' Proc.data<-get.tags.data(File)
#' plot_slopes_by_location(Proc.data=Proc.data, location=c(5.43, 52.93))
#' abline(v=as.POSIXct("2013-08-20", tz='GMT')) # end of first calibration period
#' abline(v=as.POSIXct("2014-05-05", tz='GMT')) # start of the second calibration period
#' Location<-find.stationary.location(Proc.data, '2013-07-20', '2013-08-20', initial.coords=c(10, 50))
#' }
#' @author Eldar Rakhimberdiev
#' @export
find.stationary.location<-function(Proc.data, calibration.start, calibration.stop, plot=TRUE, initial.coords=NULL, print.optimization=TRUE, reltol=1e-4) {
if (is.null(initial.coords)) stop('current function vesrion requires some inital coordinates to start search, they should not be very close but within few thousand km!')
ll_function<-function(initial.coords, Proc.data, calibration.start, calibration.stop, plot=TRUE, stage=1) {
sink("flightr_tmp")
Calibration.period<-data.frame(
calibration.start=as.POSIXct(calibration.start, tz='GMT'),
calibration.stop=as.POSIXct(calibration.stop, tz='GMT'),
lon=initial.coords[1], lat=initial.coords[2])
#use c() also for the geographic coordinates, if you have more than one calibration location
# (e. g., lon=c(5.43, 6.00), lat=c(52.93,52.94))
log.light.borders<-Proc.data$log.light.borders
log.irrad.borders<-Proc.data$log.irrad.borders
calibration.parameters<-invisible(suppressWarnings(get.calibration.parameters(Calibration.period,
Proc.data, model.ageing=FALSE,
log.light.borders=log.light.borders,
log.irrad.borders=log.irrad.borders,
plot.each = FALSE, plot.final = FALSE,
suggest.irrad.borders=FALSE)))
Twilights_total<-calibration.parameters$Twilights_total
if (length(calibration.parameters$calib_outliers)>0) {
Proc.data$FLightR.data$twilights$excluded[which(sapply(Proc.data$FLightR.data$twilights$datetime,
FUN=function(x) min(abs(calibration.parameters$calib_outliers-as.numeric(x))))<Proc.data$saving.period*24)]<-1
Proc.data_tmp<-process.twilights(Proc.data$FLightR.data$Data,
Proc.data$FLightR.data$twilights[Proc.data$FLightR.data$twilights$excluded==0,],
measurement.period=Proc.data$measurement.period, saving.period=Proc.data$saving.period,
impute.on.boundaries=Proc.data$impute.on.boundaries)
calibration.parameters<-suppressWarnings(get.calibration.parameters(Calibration.period,
Proc.data_tmp, model.ageing=FALSE,
log.light.borders=log.light.borders,
log.irrad.borders=log.irrad.borders,
plot.each = FALSE, plot.final = FALSE,
suggest.irrad.borders=FALSE))
}
if (plot) plot_slopes(calibration.parameters$All.slopes)
suppressWarnings(sink())
percent_excluded<-1-(sum(is.finite(calibration.parameters$All.slopes$Slopes$logSlope))/Twilights_total)
if (stage==1) {
if (print.optimization) message(paste(initial.coords[1], initial.coords[2], calibration.parameters$All.slopes$Parameters$LogSlope[1], calibration.parameters$All.slopes$Parameters$LogSlope[2], percent_excluded), '\n')
#print(table(calibration.parameters$All.slopes$Slopes$Type))
if (length(table(calibration.parameters$All.slopes$Slopes$Type))==1) {
warning('only_one_twilight_type_left!\n')
return(sum(is.finite(calibration.parameters$All.slopes$Slopes$logSlope)))
} else {
return(calibration.parameters$All.slopes$Parameters$LogSlope[2]^2+percent_excluded^2)
}
} else {
Dat<-calibration.parameters$All.slopes$Slopes[is.finite(calibration.parameters$All.slopes$Slopes$logSlope),]
if (length(table(Dat$Type))==1 | min(table(Dat$Type))<=2) {
warning('only_one_twilight_type_left!\n')
Val<-nrow(Dat)
} else {
#Dat<-subset(calibration.parameters$All.slopes$Slopes, is.finite(logSlope), select=c(logSlope, Time, Type))
#print(Dat)
Val<-log(1/(stats::anova(stats::lm(logSlope~Time+I(Time^2)+Type, data=Dat),stats::lm(logSlope~Time, data=Dat))[2,6])) +percent_excluded^2 #-log(1/(percent_excluded+0.001)
}
if (print.optimization) message(paste(initial.coords[1], initial.coords[2], calibration.parameters$All.slopes$Parameters$LogSlope[1], calibration.parameters$All.slopes$Parameters$LogSlope[2]), Val, percent_excluded, '\n')
return(Val)
}
}
message('stage 1...\n')
tryCatch(Res<-stats::optim(initial.coords, fn=ll_function, Proc.data=Proc.data, calibration.start=calibration.start, calibration.stop=calibration.stop, plot=plot, control=list(reltol=1e-2)), finally=try(suppressWarnings(sink())))
message('stage 2...\n')
tryCatch(Res<-stats::optim(Res$par, fn=ll_function, Proc.data=Proc.data, calibration.start=calibration.start, calibration.stop=calibration.stop, plot=plot, stage=2,control=list(reltol=reltol)), finally=try(suppressWarnings(sink())))
return(Res$par)
}
#' combines data, calibration and sets up priors
#'
#' This function is one step before \code{\link{run.particle.filter}}. It combines data, calibration, spatial extent and movement priors and estimates spatial likelihoods that used later in the particle filter.
#' @param Proc.data Processed data object created by \code{\link{get.tags.data}}.
#' @param Grid Spatial grid created by \code{\link{make.grid}}.
#' @param start release location (lat, lon).
#' @param end end of the track location. Will use \code{start} by default. Use NA in case of unknown end point.
#' @param Calibration Calibration object created by \code{\link{make.calibration}}.
#' @param threads number of parallel threads to use. default is -1, which means FLightR will use all available threads except 1. Value 1 will force sequential evaluation
#' @param Decision prior for migration probability values from 0 to 1 are allowed
#' @param Direction Direction prior for direction of migration (in degrees) with 0 pointing to the North
#' @param Kappa concentration parameter for vonMises distribution, 0 means uniform or even distribution. Will set some prior for direction for all the track, so is not recommended to be changed
#' @param M.mean Prior for mean distance travelled between consecutive twilights, km
#' @param M.sd Prior for sd of distance travelled between consecutive twilights, the higher the value is the wider is the the distribution
#' @param likelihood.correction Should likelihood correction estimated during \code{\link{make.calibration}} run be used?
#' @return Object to be uses in the \code{\link{run.particle.filter}}
#' @examples
#' File<-system.file("extdata", "Godwit_TAGS_format.csv", package = "FLightR")
#' # to run example fast we will cut the real data file by 2013 Aug 20
#' Proc.data<-get.tags.data(File, end.date=as.POSIXct('2013-07-02', tz='GMT'))
#' Calibration.periods<-data.frame(
#' calibration.start=NA,
#' calibration.stop=as.POSIXct("2013-08-20", tz='GMT'),
#' lon=5.43, lat=52.93)
#' #use c() also for the geographic coordinates, if you have more than one calibration location
#' # (e. g., lon=c(5.43, 6.00), lat=c(52.93,52.94))
#' print(Calibration.periods)
#'
#' # NB Below likelihood.correction is set to FALSE for fast run!
#' # Leave it as default TRUE for real examples
#' Calibration<-make.calibration(Proc.data, Calibration.periods, likelihood.correction=FALSE)
#'
#' Grid<-make.grid(left=0, bottom=50, right=10, top=56,
#' distance.from.land.allowed.to.use=c(-Inf, Inf),
#' distance.from.land.allowed.to.stay=c(-Inf, Inf))
#'
#' all.in<-make.prerun.object(Proc.data, Grid, start=c(5.43, 52.93),
#' Calibration=Calibration, threads=2)
#'
#' @author Eldar Rakhimberdiev
#' @export
make.prerun.object<-function(Proc.data, Grid, start, end=start, Calibration, threads=-1, Decision=0.05, Direction=0, Kappa=0, M.mean=300, M.sd=500, likelihood.correction=TRUE) {
if (length(Decision)>1) stop("Decision has to have length of 1, to sepcify it per twilight, change it in the result object of this function")
if (length(Direction)>1) stop("Direction has to have length of 1, to sepcify it per twilight, change it in the result object of this function")
if (length(Kappa)>1) stop("Kappa has to have length of 1, to sepcify it per twilight, change it in the result object of this function")
if (length(M.mean)>1) stop("M.mean has to have length of 1, to sepcify it per twilight, change it in the result object of this function")
if (length(M.sd)>1) stop("M.sd has to have length of 1, to sepcify it per twilight, change it in the result object of this function")
if (any(is.na(Proc.data$FLightR.data$twilights$datetime))) {
warning('NAs found and deleted in Proc.data$FLightR.data$twilights$datetime\n')
Proc.data$FLightR.data$twilights<- stats::na.omit(Proc.data$FLightR.data$twilights)
}
Processed.light<-make.processed.light.object(Proc.data$FLightR.data)
Index.tab<-create.proposal(Processed.light, start=start, Grid=Grid)
Index.tab$Decision<-Decision # prob of migration
Index.tab$Direction<- Direction # direction 0 - North
Index.tab$Kappa<-Kappa # distr concentration 0 means even
Index.tab$M.mean<- M.mean # distance mu
Index.tab$M.sd<- M.sd # distance sd
all.in<-geologger.sampler.create.arrays(Index.tab, Grid, start=start, stop=end)
all.in$Calibration<-Calibration
all.in$Data<-Proc.data$FLightR.data
if (threads!=1 ) {
Possible.threads<-parallel::detectCores()
if (threads<=0) Threads<-max(Possible.threads+threads, 1)
if (threads>0) Threads<-min(Possible.threads,threads)
} else {
Threads<-1
}
Phys.Mat<-get.Phys.Mat.parallel(all.in, Proc.data$Twilight.time.mat.dusk,
Proc.data$Twilight.log.light.mat.dusk,
Proc.data$Twilight.time.mat.dawn,
Proc.data$Twilight.log.light.mat.dawn,
threads=Threads, calibration=all.in$Calibration, likelihood.correction=likelihood.correction)
all.in$Spatial$Phys.Mat<-Phys.Mat
return(all.in)
}
##############
# functions below were developed before FLightR 0.3.8
# process.data
process.twilights<-function(All.p, Filtered_tw, measurement.period=60, saving.period=NULL, impute.on.boundaries=FALSE) {
# this function just prepares data for the next steps..
##########
## Dusk
# processing Dusk
Dusk.all<-Filtered_tw$datetime[Filtered_tw$type==2]
Twilight.index.mat.dusk<-sapply(which(All.p$gmt %in% Dusk.all & All.p$type==2), FUN=function(x) (x-24):(x+24))
Twilight.index.mat.dusk<-apply(Twilight.index.mat.dusk, c(1,2), FUN=function(x) ifelse (x>0, x, NA))
Max.Index<-nrow(All.p)
Twilight.index.mat.dusk<-apply(Twilight.index.mat.dusk, c(1,2), FUN=function(x) ifelse (x>Max.Index, NA, x))
Twilight.time.mat.dusk<-apply(Twilight.index.mat.dusk, c(1,2), FUN=function(x) as.numeric(All.p$gmt[x]))
Twilight.time.mat.dusk<-apply(Twilight.time.mat.dusk, c(1,2), FUN=function(x) ifelse(is.finite(x), x, 0))
Twilight.log.light.mat.dusk<-apply(Twilight.index.mat.dusk, c(1,2), FUN=function(x) log(All.p$light[x]))
#Twilight.log.light.mat.dusk<-apply(Twilight.index.mat.dusk, c(1,2), FUN=function(x) All.p$light[x])
Twilight.log.light.mat.dusk<-apply(Twilight.log.light.mat.dusk, c(1,2), FUN=function(x) ifelse(is.finite(x), x, -1))
Twilight.time.mat.dusk<-Twilight.time.mat.dusk-(saving.period-measurement.period)
# processing Dawn
Dawn.all<-Filtered_tw$datetime[Filtered_tw$type==1]
Twilight.index.mat.dawn<-sapply(which(All.p$gmt %in% Dawn.all & All.p$type==1), FUN=function(x) (x-24):(x+24))
Twilight.index.mat.dawn<-apply(Twilight.index.mat.dawn, c(1,2), FUN=function(x) ifelse (x>0, x, NA))
Max.Index<-nrow(All.p)
Twilight.index.mat.dawn<-apply(Twilight.index.mat.dawn, c(1,2), FUN=function(x) ifelse (x>Max.Index, NA, x))
Twilight.time.mat.dawn<-apply(Twilight.index.mat.dawn, c(1,2), FUN=function(x) as.numeric(All.p$gmt[x]))
Twilight.time.mat.dawn<-apply(Twilight.time.mat.dawn, c(1,2), FUN=function(x) ifelse(is.finite(x), x, 0))
Twilight.log.light.mat.dawn<-apply(Twilight.index.mat.dawn, c(1,2), FUN=function(x) log(All.p$light[x]))
#Twilight.log.light.mat.dawn<-apply(Twilight.index.mat.dawn, c(1,2), FUN=function(x) All.p$light[x])
Twilight.log.light.mat.dawn<-apply(Twilight.log.light.mat.dawn, c(1,2), FUN=function(x) ifelse(is.finite(x), x, -1))
Res<-list(Twilight.time.mat.dusk=Twilight.time.mat.dusk, Twilight.log.light.mat.dusk=Twilight.log.light.mat.dusk, Twilight.time.mat.dawn=Twilight.time.mat.dawn, Twilight.log.light.mat.dawn=Twilight.log.light.mat.dawn, measurement.period=measurement.period, saving.period=saving.period, impute.on.boundaries=impute.on.boundaries)
return(Res)
}
get.Irradiance<-function(alpha, r=6378, s=6.9, intigeo.template.correction=FALSE) {
# function from Ekstrom 2007
erf <- function(x) 2 * stats::pnorm(x * sqrt(2)) - 1
## (see Abramowitz and Stegun 29.2.29)
## and the so-called 'complementary error function'
erfc <- function(x) 2 * stats::pnorm(x * sqrt(2), lower = FALSE)
## and the inverses
#erfinv <- function (x) qnorm((1 + x)/2)/sqrt(2)
#erfcinv <- function (x) qnorm(x/2, lower = FALSE)/sqrt(2)
Res<-alpha
u<-sqrt(r/(2*s))*sin(alpha)
Res[(u<=0)]<-(exp(-u^2)/(1+erf(-u)))[(u<=0)]
Res[(u>0)]<-(exp(-u^2)/(erfc(u)))[(u>0)]
if (intigeo.template.correction) Res[Res>0.001]=exp(log(Res[Res>0.001])*1.279+log(Res[Res>0.001])^2*0.091)
return(Res)
}
logger.template.calibrarion.internal<-function( Twilight.time.mat.Calib.dawn, Twilight.log.light.mat.Calib.dawn, Twilight.time.mat.Calib.dusk, Twilight.log.light.mat.Calib.dusk, positions=NA, plot.each=TRUE, plot.final=TRUE,log.light.borders=NA, log.irrad.borders=c(-8, 1.3), impute.on.boundaries=FALSE) {
oldpar <- graphics::par(no.readonly = TRUE)
on.exit(graphics::par(oldpar))
# =================
# in this function I'll add a new lnorm calibration...
#
# now I'll try to do them both..
if (!is.list(positions)) {
if (length (positions) !=2) stop("the positions should be a list with $dawn and $dusk dataframes with pairs of coordiantes OR just one pair\n")
dawn=matrix(ncol=2, nrow=dim(Twilight.time.mat.Calib.dawn)[2])
dawn[,1]<-positions[1]
dawn[,2]<-positions[2]
#
dusk=matrix(ncol=2, nrow=dim(Twilight.time.mat.Calib.dusk)[2])
dusk[,1]<-positions[1]
dusk[,2]<-positions[2]
#
positions=list(dawn=dawn, dusk=dusk)
}
Twilight.time.mat.Calib.dawn<-Twilight.time.mat.Calib.dawn[-25,]
Twilight.log.light.mat.Calib.dawn<-Twilight.log.light.mat.Calib.dawn[-25,]
Twilight.time.mat.Calib.dusk<-Twilight.time.mat.Calib.dusk[-25,]
Twilight.log.light.mat.Calib.dusk<-Twilight.log.light.mat.Calib.dusk[-25,]
# let's try to create Calib.data.all first!!
Calib.data.dawn<-data.frame()
if (plot.each) graphics::par(ask=FALSE)
for (dawn in 1:dim(Twilight.time.mat.Calib.dawn)[2]) {
message("\r checking dawn", dawn)
#Twilight.solar.vector<-solar(as.POSIXct(Twilight.time.mat.Calib.dawn[, dawn], tz="gmt", origin="1970-01-01"))
Data<-check.boundaries(positions$dawn[dawn,], Twilight.solar.vector=NULL, Twilight.log.light.vector = Twilight.log.light.mat.Calib.dawn[,dawn], plot=plot.each, verbose=FALSE, log.light.borders=log.light.borders, log.irrad.borders=log.irrad.borders, dusk=FALSE, Twilight.time.vector=Twilight.time.mat.Calib.dawn[, dawn], impute.on.boundaries=impute.on.boundaries)
if (length(Data)==0) {
message("\r dawn", dawn, "was excluded from the calibration")
} else {
Calib.data.dawn<-rbind(Calib.data.dawn, cbind(LogLight=Data[,1], LogIrrad=Data[,2], Day=dawn, Time=Data[,3], Elevs=Data[,4]))
}
}
if (nrow(Calib.data.dawn)>0) {Calib.data.dawn$type<-"Dawn"}
Calib.data.dusk<-data.frame()
for (dusk in 1:dim(Twilight.time.mat.Calib.dusk)[2]) {
message("\r checking dusk", dusk )
#Twilight.solar.vector<-solar(as.POSIXct(Twilight.time.mat.Calib.dusk[, dusk], tz="gmt", origin="1970-01-01"))
Data<-check.boundaries(positions$dusk[dusk,], Twilight.solar.vector=NULL, Twilight.log.light.vector=Twilight.log.light.mat.Calib.dusk[,dusk], plot=plot.each, verbose=FALSE, log.light.borders=log.light.borders, log.irrad.borders=log.irrad.borders, dusk=TRUE, Twilight.time.vector=Twilight.time.mat.Calib.dusk[, dusk], impute.on.boundaries=impute.on.boundaries)
if (length(Data)==0) {
message("\r dusk", dusk, "was excluded from the calibration")
} else {
Calib.data.dusk<-rbind(Calib.data.dusk, cbind(LogLight=Data[,1], LogIrrad=Data[,2], Day=dim(Twilight.time.mat.Calib.dawn)[2]+dusk, Time=Data[,3], Elevs=Data[,4]))
}
}
if (nrow(Calib.data.dusk)>0) {Calib.data.dusk$type<-"Dusk"}
message('\r')
if (nrow(Calib.data.dusk)==0 & nrow(Calib.data.dawn)==0) stop('Error: None of twilights could happen in the location... Something went wrong at the previous stages')
Calib.data.all<-rbind(Calib.data.dawn, Calib.data.dusk)
if (plot.final) {
graphics::plot(LogLight~LogIrrad,data=Calib.data.all, type="n")
for (Day in unique(Calib.data.all$Day)) {
graphics::lines(LogLight~LogIrrad, data=Calib.data.all[Calib.data.all$Day==Day,], type="b", pch="+", col=grDevices::rainbow(length(unique(Calib.data.all$Day)))[Day])
}
}
Calib.data.all$fDay<-as.factor(Calib.data.all$Day)
message('\r\n')
message('\r\n')
return(Calib.data.all)
}
# ok now we need template.calibration.function..
logger.template.calibration<-function(Twilight.time.mat.Calib.dawn, Twilight.log.light.mat.Calib.dawn, Twilight.time.mat.Calib.dusk, Twilight.log.light.mat.Calib.dusk, time.shift=0, positions, log.light.borders=NA, log.irrad.borders=c(-15, 50), plot.each=TRUE, plot.final=TRUE, impute.on.boundaries=FALSE) {
Calibration.original<-logger.template.calibrarion.internal(Twilight.time.mat.Calib.dawn+time.shift, Twilight.log.light.mat.Calib.dawn, Twilight.time.mat.Calib.dusk+time.shift, Twilight.log.light.mat.Calib.dusk, positions=positions, log.light.borders=log.light.borders, log.irrad.borders= log.irrad.borders, plot.each=plot.each, plot.final=plot.final,impute.on.boundaries=impute.on.boundaries)
return(Calibration.original)
}
get.calib.param<-function(Calib.data.all, plot=FALSE, calibration.type=NULL) {
if (is.null(calibration.type)) calibration.type="parametric.slope"
message("calibration method used:", calibration.type, "\n")
Calib.data.all$fTwilight<-Calib.data.all$fDay
#============================
# we should do some outlier tests...
cur.data<-Calib.data.all
p.lm1<-stats::lm(LogLight~fTwilight/LogIrrad,data=cur.data)
cur.slope<-matrix(stats::coef(p.lm1),ncol=2) # foo[,2] - slope
cur.slope<-as.data.frame(cur.slope)
Diag<-diag(stats::vcov(p.lm1))
Diag.slopes<-Diag[which(sapply(strsplit(names(Diag), ":"), length)==2)]
cur.slope$sd[!is.na(cur.slope[,2])]<-sqrt(Diag.slopes)
cur.slope$type<-stats::aggregate(cur.data[,"type"],by=list(Day=cur.data$fTwilight),FUN=function(x) unique(x))[,2]
cur.slope$time<-stats::aggregate(cur.data[,"Time"],by=list(Day=cur.data$fTwilight),FUN=function(x) x[1])[,2]
cur.slope$Duration<-stats::aggregate(cur.data[,"Time"],by=list(Day=cur.data$fTwilight),FUN=function(x) max(x)-min(x))[,2]
names(cur.slope)[2]<-"slope"
#===========================
# ok I've decided to make separate lm
Slopes<-c()
Slopes.sd<-c()
Intercept=c()
Sigma=c()
Type=c()
Time=c()
Elevs<-c()
Day<-c()
for (i in (unique(cur.data$fTwilight))) {
# lm
#plot(LogLight~LogIrrad, data=cur.data[cur.data$fTwilight==i,])
Data<-cur.data[cur.data$fTwilight==i,]
Data<-Data[Data$LogLight>0,]
Lm<-stats::lm(LogLight~LogIrrad,data=Data)
Cur_slope_mean<-stats::coef(Lm)[2]
Cur_slope_sd<-sqrt(stats::vcov(Lm)[4])
if (calibration.type=="nonparametric.slope") {
# Here I want to physically estimate slopes
All_slopes_cur<-diff(Data$LogLight)/diff(Data$LogIrrad)
All_slopes_mean<-mean(All_slopes_cur)
All_slopes_sd<-stats::sd(All_slopes_cur)
# exclude outliers
All_slopes_diff<-abs(All_slopes_cur-All_slopes_mean)
Slopes_outliers<-which(All_slopes_diff>3*All_slopes_sd)
if (length(Slopes_outliers)>0) {
All_slopes_cur<-All_slopes_cur[-Slopes_outliers]
All_slopes_mean<-mean(All_slopes_cur)
All_slopes_sd<-stats::sd(All_slopes_cur)
}
Cur_slope_mean<-All_slopes_mean
Cur_slope_sd<-All_slopes_sd
}
Slopes<-c(Slopes, Cur_slope_mean)
Slopes.sd<-c(Slopes.sd, Cur_slope_sd)
Intercept=c(Intercept, stats::coef(Lm)[1])
Sigma<-c(Sigma, summary(Lm)$sigma)
Type=c(Type, cur.data$type[cur.data$fTwilight==i][1])
Time<-c(Time,ifelse(cur.data$type[cur.data$fTwilight==i][1]=="Dusk", max(cur.data$Time[cur.data$fTwilight==i]), min(cur.data$Time[cur.data$fTwilight==i])))
Elevs=c(Elevs, mean(cur.data$Elevs[cur.data$fTwilight==i], na.rm=TRUE))
Day=c(Day, cur.data$Day[cur.data$fTwilight==i][1])
# nonparametric slope estimation
}
#plot(cur.slope$slope~Slopes)
#plot(cur.slope$sd~Slopes.sd)
# ok, we now want to replace old slopes with new ones.
# but first we want to add 0.5 instead of NA - at least this is how it works now inside the functions and how it was done in a simulations
cur.slope$slope<-Slopes
cur.slope$Intercept<-Intercept
cur.slope$Sigma<-Sigma
#Slopes.sd[is.na(Slopes.sd)]<-sd.fast
cur.slope$sd<-Slopes.sd
names(cur.slope)[2]<-"slope"
#=====================
if (plot) graphics::hist(log(cur.slope$slope))
cur.slope.int<-cur.slope[is.finite(log(cur.slope$slope)),]
Parameters<-list(Intercept=c(mean(cur.slope.int$Intercept, na.rm=TRUE), stats::sd(cur.slope.int$Intercept, na.rm=TRUE)), LogSlope=c(mean(log(cur.slope.int$slope), na.rm=TRUE), stats::sd(log(cur.slope.int$slope), na.rm=TRUE)), LogSigma=c(mean(log(cur.slope.int$Sigma[!is.na(cur.slope.int$sd)])), stats::sd(log(cur.slope.int$Sigma[!is.na(cur.slope.int$sd)]))), mean.of.individual.slope.sigma=mean(cur.slope.int$sd, na.rm=TRUE), calibration.type=calibration.type)
#cur.slope$time<-aggregate(cur.data[,"Time"],by=list(Day=cur.data$fTwilight),FUN=function(x) x[1])[,2]
#cur.slope$time<-aggregate(cur.data[,"Time"],by=list(Day=cur.data$fTwilight),FUN=mean)[,2]
Res<-list(Parameters=Parameters, Slopes=data.frame(Slope=cur.slope$slope, Time=Time, Intercept=cur.slope$Intercept, Sigma=cur.slope$Sigma, Slopes.sd=Slopes.sd, Type=Type, Elevs=Elevs, Day=Day))
return(Res)
}
plot_slopes<-function(all.slopes, ylim=NULL, xlim=NULL) {
oldpar <- graphics::par(no.readonly = TRUE)
all.slopes$Slopes<-all.slopes$Slopes[all.slopes$Slopes$Slope>0,]
if (is.null(xlim)) {
graphics::plot(log(all.slopes$Slopes$Slope)~as.POSIXct(all.slopes$Slopes$Time, tz="GMT", origin="1970-01-01"), type="n", main="red - dawn, black - dusk", xlab="time", ylab="log(Slope)", ylim=ylim, las=1)
} else {
graphics::plot(log(all.slopes$Slopes$Slope)~ as.POSIXct(all.slopes$Slopes$Time, tz="GMT", origin="1970-01-01"), type="n", main="red - dawn, black - dusk", xlab="time", ylab="log(Slope)", ylim=ylim, las=1, xlim= as.POSIXct(xlim, tz='GMT'))
}
graphics::lines(log(Slope)~ as.POSIXct(Time, tz="GMT", origin="1970-01-01"), data=all.slopes$Slopes[all.slopes$Slopes$Type=="Dusk",])
graphics::points(log(Slope)~ as.POSIXct(Time, tz="GMT", origin="1970-01-01"), data=all.slopes$Slopes[all.slopes$Slopes$Type=="Dusk",], pch="+")
graphics::points(log(Slope)~ as.POSIXct(Time, tz="GMT", origin="1970-01-01"), data=all.slopes$Slopes[all.slopes$Slopes$Type=="Dawn",], pch="+", col="red")
graphics::lines(log(Slope)~ as.POSIXct(Time, tz="GMT", origin="1970-01-01"), data=all.slopes$Slopes[all.slopes$Slopes$Type=="Dawn",], col="red")
oldpar$usr<-graphics::par()$usr
on.exit(graphics::par(oldpar))
#invisible()
return(NULL)
}
get.calibration.parameters<-function(Calibration.periods, Proc.data, model.ageing=TRUE, log.light.borders=log(c(2, 63)), log.irrad.borders=c(-7,1.5), plot.each=FALSE, plot.final=TRUE, calibration.type=NULL, suggest.irrad.borders=FALSE) {
# ageing.model this option should be used only in case there are two periods of calibration or there is one but very long.
if (nrow(Calibration.periods)==1 & model.ageing) warning("you have only one calibration period ageing estimation is unreliable and should be turned to F!!!")
Calibration.periods.int<-Calibration.periods
Calibration.periods.int$calibration.start<-as.numeric(Calibration.periods.int$calibration.start)
Calibration.periods.int$calibration.stop<-as.numeric(Calibration.periods.int$calibration.stop)
Dusk.calib.days<-c()
Dawn.calib.days<-c()
Positions=list(dawn=c(), dusk=c())
#sum(apply(Calibration.periods.int, 1, FUN=function(x) x[2]-x[1]))/3600/24
for (i in 1:nrow(Calibration.periods)) {
Dusk.calib.days.cur<-NULL
Dusk.calib.days.cur<-which(Proc.data$Twilight.time.mat.dusk[25,]>Calibration.periods.int$calibration.start[i] & Proc.data$Twilight.time.mat.dusk[25,] < Calibration.periods.int$calibration.stop[i])
if (length(Dusk.calib.days.cur) > 0) {
Dusk.calib.days<-unique(c(Dusk.calib.days, Dusk.calib.days.cur))
dusk.cur=matrix(ncol=2, nrow=length(Dusk.calib.days.cur))
dusk.cur[,1]<-Calibration.periods.int$lon[i]
dusk.cur[,2]<-Calibration.periods.int$lat[i]
Positions$dusk<-rbind(Positions$dusk, dusk.cur)
}
Dawn.calib.days.cur<-NULL
Dawn.calib.days.cur<-which(Proc.data$Twilight.time.mat.dawn[25,]>Calibration.periods.int$calibration.start[i] & Proc.data$Twilight.time.mat.dawn[25,] < Calibration.periods.int$calibration.stop[i] )
if (length(Dawn.calib.days.cur) > 0) {
Dawn.calib.days<-unique(c(Dawn.calib.days,Dawn.calib.days.cur ))
dawn.cur=matrix(ncol=2, nrow=length(Dawn.calib.days.cur))
dawn.cur[,1]<-Calibration.periods.int$lon[i]
dawn.cur[,2]<-Calibration.periods.int$lat[i]
Positions$dawn<-rbind(Positions$dawn, dawn.cur)
}
}
Twilight.time.mat.Calib.dusk<-Proc.data$Twilight.time.mat.dusk[,Dusk.calib.days]
Twilight.log.light.mat.Calib.dusk<-Proc.data$Twilight.log.light.mat.dusk[,Dusk.calib.days]
## Dawn
Twilight.time.mat.Calib.dawn<-Proc.data$Twilight.time.mat.dawn[,Dawn.calib.days]
Twilight.log.light.mat.Calib.dawn<-Proc.data$Twilight.log.light.mat.dawn[,Dawn.calib.days ]
Twilights_total<-length(c(Dusk.calib.days,Dawn.calib.days))
#-------------------#
if (suggest.irrad.borders) {
Calib.data.all<-logger.template.calibration(Twilight.time.mat.Calib.dawn, Twilight.log.light.mat.Calib.dawn, Twilight.time.mat.Calib.dusk, Twilight.log.light.mat.Calib.dusk, positions=Positions, log.light.borders=log.light.borders, log.irrad.borders=c(-10,10), plot.each=plot.each, plot.final=plot.final, impute.on.boundaries=Proc.data$impute.on.boundaries)
log.irrad.borders<-suggest.irrad.boundaries(Calib.data.all)
warning('!!! log.irrad.borders were updated to new values:', log.irrad.borders[1], log.irrad.borders[2], '\n')
}
Calib.data.all<-logger.template.calibration(Twilight.time.mat.Calib.dawn, Twilight.log.light.mat.Calib.dawn, Twilight.time.mat.Calib.dusk, Twilight.log.light.mat.Calib.dusk, positions=Positions, log.light.borders=log.light.borders, log.irrad.borders=log.irrad.borders, plot.each=plot.each, plot.final=plot.final, impute.on.boundaries=Proc.data$impute.on.boundaries)
All.slopes<-suppressWarnings(get.calib.param(Calib.data.all, plot=FALSE, calibration.type=calibration.type))
All.slopes$Slopes$logSlope<-suppressWarnings(log(All.slopes$Slopes$Slope))
if (model.ageing) {
All.slopes.int<-All.slopes
All.slopes.int$Slopes<-All.slopes.int$Slopes[is.finite(All.slopes.int$Slopes$logSlope),]
Time.start=min(c(Twilight.time.mat.Calib.dusk, Twilight.time.mat.Calib.dawn))
Model=stats::lm(logSlope~I(Time-Time.start), data=All.slopes.int$Slopes)
Model$Time.start<-Time.start
calib_outliers<-All.slopes.int$Slopes$Time[which(abs(stats::residuals(Model))>3*stats::sd(stats::residuals(Model))*sqrt(length(stats::residuals(Model))))]
Res<-list(calib_outliers=calib_outliers, ageing.model=Model, All.slopes=All.slopes,log.irrad.borders=log.irrad.borders)
} else {
All.slopes.int<-All.slopes
All.slopes.int$Slopes<-All.slopes.int$Slopes[is.finite(All.slopes.int$Slopes$logSlope),]
calib_outliers<-All.slopes.int$Slopes$Time[which(abs(All.slopes.int$Slopes$logSlope-mean(All.slopes.int$Slopes$logSlope, na.rm=TRUE))>3*stats::sd(All.slopes.int$Slopes$logSlope, na.rm=TRUE))]
Res<-list(calib_outliers=calib_outliers, All.slopes=All.slopes, log.irrad.borders=log.irrad.borders, Twilights_total=Twilights_total)
}
return(Res)
}
suggest.irrad.boundaries<-function(Calib.data.all, leave.out=0.01) {
Model<-stats::lm(LogLight~LogIrrad + fDay, data=Calib.data.all)
MinLight<-round(min(Calib.data.all$LogLight),1)
Left.border<- -(stats::coef(Model)[1] + stats::quantile(stats::coef(Model)[3:length(stats::coef(Model))], leave.out) - MinLight)*stats::coef(Model)[2]
MaxLight<-round(max(Calib.data.all$LogLight),1)
Right.border<-(MaxLight- (stats::coef(Model)[1] + stats::quantile(stats::coef(Model)[3:length(stats::coef(Model))], 1-leave.out))) *stats::coef(Model)[2]
return(c(Left.border, Right.border))
}
make_likelihood_correction_function<-function(calib_log_mean, calib_log_sd, cur_sd_range, npoints=300, plot=FALSE) {
cur_mean_range<-c(exp(calib_log_mean)-5*exp(calib_log_sd), exp(calib_log_mean)+5*exp(calib_log_sd))
Res<-c()
for (i in 1:npoints) {
message('\r simulation: ', round(i/npoints*100), '%', sep='')
message('\r')
cur_mean<-stats::runif(300, cur_mean_range[1], cur_mean_range[2])
cur_sd<-stats::runif(1, cur_sd_range[1], cur_sd_range[2])
Cur_max_real<- stats::optimize(f=function(x)
mean(stats::dlnorm(stats::rnorm(200000,x , cur_sd),calib_log_mean, calib_log_sd)),
interval=cur_mean_range, maximum = TRUE)$maximum
Res<-rbind(Res, cbind(calib_log_mean, calib_log_sd, cur_sd, cur_mean_max=Cur_max_real))
}
Res<-as.data.frame(Res)
Res$Corr <- exp(Res$calib_log_mean)-Res$cur_mean_max
message('\r estimating correction function...')
message('\r')
MvExp<-mgcv::gamm(Corr~s(cur_sd), data=Res, weights=nlme::varExp(form =~ cur_sd))
# now we check for the outliers..
Resid<-stats::resid(MvExp$lme, type='normalized')
Outliers<-NULL
Outliers<-which(abs(Resid)>3*stats::sd(Resid))
if (length(Outliers)>0) {
Res<-Res[-Outliers,]
MvExp<-mgcv::gamm(Corr~s(cur_sd), data=Res, weights=nlme::varExp(form =~ cur_sd))
}
XX<-seq(cur_sd_range[1], cur_sd_range[2], by=0.01)
c_fun<-stats::approxfun(mgcv::predict.gam(MvExp$gam, newdata=data.frame(cur_sd=XX)) ~XX)
if (plot) {
graphics::plot(Res[,4]~Res[,3], pch='+', ylab='cur_mean', xlab='cur_sd')
graphics::lines(exp(Res$calib_log_mean[1])-mgcv::predict.gam(MvExp$gam, newdata=data.frame(cur_sd=XX))~XX, col='red')
}
Out<-list(c_fun=c_fun, Res=Res)
message('\r')
message('\n')
return(Out)
}
create.calibration<-function( All.slopes, Proc.data, FLightR.data, location, log.light.borders, log.irrad.borders, ageing.model=NULL, likelihood.correction=TRUE) {
# Now we create 'parameters' object that will have all the details about the calibration
Parameters<-All.slopes$Parameters # LogSlope #
Parameters$measurement.period<-Proc.data$measurement.period
Parameters$saving.period<-Proc.data$saving.period
Parameters$log.light.borders<-log.light.borders # these are the boundaries in which one should use the BAS tag.. for the other types of tags they will be different.
Parameters$min.max.values<-c(min(FLightR.data$Data$light), max(FLightR.data$Data$light))
Parameters$log.irrad.borders=log.irrad.borders
Parameters$impute.on.boundaries=Proc.data$impute.on.boundaries
Parameters$location=location
Parameters$LogSlope_1_minute<-Parameters$LogSlope
if (is.null(ageing.model)) {
time_correction_fun= eval(parse(text=paste("function (x,y) return(", Parameters$LogSlope[1], ")")))
lat_correction_fun<-function(x, y, z) return(0)
} else {
time_correction_fun= function(x, y) return(0)
lat_correction_fun<-eval(parse(text=paste("function (x,y) return(", stats::coef(ageing.model)[1], "+", stats::coef(ageing.model)[2], "* (y-",ageing.model$Time.start," ))")))
}
c_fun=NULL
if (likelihood.correction) {
if (is.na(All.slopes$Parameters$LogSigma[1])) {
cur_sd_range=c(0, 1)
} else {
cur_sd_range=c(0, max(exp(All.slopes$Parameters$LogSigma[1]+5*All.slopes$Parameters$LogSigma[2]),1))
}
c_fun=make_likelihood_correction_function(Parameters$LogSlope[1], Parameters$LogSlope[2], cur_sd_range=cur_sd_range)$c_fun
}
Calibration<-list(Parameters=Parameters, time_correction_fun=time_correction_fun, lat_correction_fun=lat_correction_fun, c_fun=c_fun)
return(Calibration)
}
correct.hours<-function(datetime) {
# this function is supposed to correct hours by adding some amount of them.
hours <- as.numeric(format(datetime,"%H"))+as.numeric(format(datetime,"%M"))/60
hours[as.POSIXlt(datetime, tz='GMT')$isdst==1] <-hours[as.POSIXlt(datetime, tz='GMT')$isdst==1]-1
cor <- rep(NA, 24)
for(i in 0:23){
cor[i+1] <- max(abs((c(hours[1],hours)+i)%%24 -
(c(hours,hours[length(hours)])+i)%%24),na.rm=TRUE)
}
hours <- (hours + (which.min(round(cor,2)))-1)%%24
return(hours)
}
make.result.list<-function(Data, raw.X, raw.Y) {
Int<-c(min(raw.Y), max(raw.Y))
Index<-sapply(raw.X, FUN=function(x) which.min(abs(as.numeric(Data$d$gmt-x))))
Res<-Data$d[Index,]
# now I want to adjust time and light
Res$light<-stats::approx(x=Data$d$gmt, y=Data$d$light, xout=raw.X)$y
Res$gmt<-raw.X
Res$Hour<-raw.Y
Result<-list(Data=Res)
Result$Data$gmt.adj<-Result$Data$gmt
Result$Data$gmt<- as.POSIXct(Result$Data$gmt, tz="GMT", origin="1970-01-01")
Result$Data$gmt.adj<- as.POSIXct(Result$Data$gmt.adj, tz="GMT", origin="1970-01-01")
return(Result)
}
make.processed.light.object<-function(FLightR.data) {
# dusk
raw.Y.dusk<-correct.hours(FLightR.data$twilights$datetime[FLightR.data$twilights$type==2 & FLightR.data$twilights$excluded==0])
raw.X.dusk<-as.numeric(FLightR.data$twilights$datetime[FLightR.data$twilights$type==2 & FLightR.data$twilights$excluded==0])
Data_tmp<-list(d=FLightR.data$Data)
Result.Dusk<-make.result.list(Data_tmp, raw.X.dusk, raw.Y.dusk)
# dusk
raw.Y.dawn<-correct.hours(FLightR.data$twilights$datetime[FLightR.data$twilights$type==1 & FLightR.data$twilights$excluded==0])
raw.X.dawn<-as.numeric(FLightR.data$twilights$datetime[FLightR.data$twilights$type==1 & FLightR.data$twilights$excluded==0])
Result.Dawn<-make.result.list(Data_tmp, raw.X.dawn, raw.Y.dawn)
# combine
processed.light<-list(Final.dusk=Result.Dusk, Final.dawn=Result.Dawn)
return(processed.light)
}
# this create proposal is corrected for the missing loess
create.proposal<-function(processed.light, start=c(-98.7, 34.7), end=NA, Grid) {
if (is.na(end)) end=start
All.Days<-seq(min(min(processed.light$Final.dusk$Data$gmt), min(processed.light$Final.dawn$Data$gmt)),max(max(processed.light$Final.dusk$Data$gmt), max(processed.light$Final.dawn$Data$gmt)), by="days")
## Probabilty.of.migration
All.Days.extended<-seq(min(All.Days)-86400, max(All.Days)+86400, by="days")
# ver 0.4.2
# correct twilights for polar periods..
start_no_polar<-start
start_no_polar[2]<-min(abs(start_no_polar[2]), 65.73) * sign(start_no_polar[2])
# these are all potential twilights that we would have for the start point...
Potential.twilights<-sort(c(suntools::sunriset(matrix(start_no_polar, nrow=1), All.Days.extended, direction="sunrise", POSIXct.out=TRUE)[,2], suntools::sunriset(matrix(start_no_polar, nrow=1), All.Days.extended, direction="sunset", POSIXct.out=TRUE)[,2]))
# now we need to add dask or dawn to it..
Sunrise<-suntools::sunriset(matrix(start_no_polar, nrow=1), Potential.twilights[1], direction="sunrise", POSIXct.out=TRUE)[,2]
Sunrises<-c(Sunrise-(3600*24), Sunrise, Sunrise+(3600*24))
Sunset<-suntools::sunriset(matrix(start_no_polar, nrow=1), Potential.twilights[1], direction="sunset", POSIXct.out=TRUE)[,2]
Sunsets<-c(Sunset-(3600*24), Sunset, Sunset+(3600*24))
First.twilight<-ifelse(which.min(c(min(abs(difftime(Sunrises,Potential.twilights[1], units="mins"))), min(abs(difftime(Sunsets,Potential.twilights[1], units="mins")))))==1, "dawn", "dusk")
Index.tab<-data.frame(Date=Potential.twilights)
if (First.twilight=="dusk") {
Index.tab$Dusk<-rep(c(TRUE,FALSE), times=length(All.Days.extended))
} else {
Index.tab$Dusk<-rep(c(FALSE,TRUE), times=length(All.Days.extended))
}
Index.tab$Curr.mat<-NA # this will be NA if no data and row number if there are..
Index.tab$Real.time<- as.POSIXct(NA, tz="GMT")
Index.tab$time<- as.POSIXct(NA, tz="GMT")
######################################
# !!!! this will not work if bird will move for over 12 time zones!!!
# Dusk
for (i in 1:length(processed.light$Final.dusk$Data$gmt)) {
Row2write<-which.min(abs(difftime(processed.light$Final.dusk$Data$gmt[i], Index.tab$Date[Index.tab$Dusk==TRUE], units="mins")))
Index.tab$time[Index.tab$Dusk==TRUE][Row2write]<-processed.light$Final.dusk$Data$gmt[i]
Index.tab$Real.time[Index.tab$Dusk==TRUE][Row2write]<-processed.light$Final.dusk$Data$gmt.adj[i]
Index.tab$Curr.mat[Index.tab$Dusk==TRUE][Row2write]<-i
}
# Dawn
for (i in 1:length(processed.light$Final.dawn$Data$gmt)) {
Row2write<-which.min(abs(difftime(processed.light$Final.dawn$Data$gmt[i], Index.tab$Date[Index.tab$Dusk==FALSE], units="mins")))
Index.tab$time[Index.tab$Dusk==FALSE][Row2write]<-processed.light$Final.dawn$Data$gmt[i]
Index.tab$Real.time[Index.tab$Dusk==FALSE][Row2write]<-processed.light$Final.dawn$Data$gmt.adj[i]
Index.tab$Curr.mat[Index.tab$Dusk==FALSE][Row2write]<-i
}
# cutting empty ends..
while (is.na(Index.tab$Curr.mat[1])) Index.tab<-Index.tab[-1,]
while (is.na(Index.tab$Curr.mat[nrow(Index.tab)])) Index.tab<-Index.tab[-nrow(Index.tab),]
Index.tab$Point<-NA
#First.Point<-which.min(sp::spDistsN1(Grid[,1:2], start, longlat=TRUE))
Grid_sf<-Grid[,1:2] |> sf::st_multipoint() |> sf::st_sfc() |> sf::st_cast('POINT')
First.Point<-sf::st_nearest_feature(sf::st_point(start),Grid_sf, longlat = TRUE)
Index.tab$Point[1]<-First.Point
#
# I decided that Curr.mat is not needed anymore
#Index.tab$Main.Index[,-which(names(Index.tab$Main.Index)=="Curr.mat")]
# this need to double checked
Index.tab$yday<-as.POSIXlt(Index.tab$Date, tz="GMT")$yday
return(Index.tab)
}
geologger.sampler.create.arrays<-function(Index.tab, Grid, start, stop=start) {
# this function wil take in Index.table with all the proposals and will create an array for case of missed data..
# the main feature is that it can account for missing data..
Index.tab.old<-Index.tab
Index.tab$proposal.index<-NA
Index.tab$proposal.index[Index.tab$Dusk==TRUE]<-"Dusk"
Index.tab$proposal.index[Index.tab$Dusk==FALSE]<-"Dawn"
Missed.twilights<-which(is.na(Index.tab$Curr.mat))
if (length(Missed.twilights)>0) {
Deleted.Rows.Count<-0
for (i in Missed.twilights) {
i=i-Deleted.Rows.Count
Index.tab$proposal.index[i-1]<-"Comb"
Index.tab$Decision[i-1]<-(Index.tab$Decision[i-1]+Index.tab$Decision[i])/2 #
if (Index.tab$Direction[i-1] != Index.tab$Direction[i]) {
Index.tab$Direction[i-1]<-0
Index.tab$Kappa[i-1]<-0
}
Index.tab$M.mean[i-1]<-Index.tab$M.mean[i-1]+Index.tab$M.mean[i]
Index.tab$M.sd[i-1]<-sqrt((Index.tab$M.sd[i-1])^2+(Index.tab$M.sd[i])^2)
Index.tab<-Index.tab[-i,]
Deleted.Rows.Count=Deleted.Rows.Count+1
}
}
output<-list()
Matrix.Index.Table<- Index.tab[,which(names(Index.tab) %in% c("Decision", "Direction", "Kappa", "M.mean", "M.sd", "yday", "Real.time", "time", "Dusk", "Loess.se.fit", "Loess.n")) ]
# I also remove last line as bird was not flying after last twilight
Matrix.Index.Table<-Matrix.Index.Table[-nrow(Matrix.Index.Table),]
# main index will have the same amount of rows we have in twilight matrices without first..
Main.Index<-c()
Main.Index$Biol.Prev<-1:nrow(Matrix.Index.Table)
Main.Index$proposal.index<-Index.tab$proposal.index[-nrow(Index.tab)]
Main.Index<-as.data.frame(Main.Index)
Indices<-list(Matrix.Index.Table=Matrix.Index.Table, Main.Index=Main.Index)
output$Indices=Indices
# ok now we need to add to output all other parts..
output$Spatial<-list(Grid=Grid) #[,1:2]
output$Spatial$Behav.mask<-as.integer(Grid[,3])
#output$Spatial$start.point<-which.min(sp::spDistsN1(Grid[,1:2], start, longlat=TRUE))
Grid_sf<-Grid[,1:2] |> sf::st_multipoint() |> sf::st_sfc() |> sf::st_cast('POINT')
output$Spatial$start.point<-sf::st_nearest_feature(sf::st_point(start),Grid_sf, longlat = TRUE)
output$Spatial$start.location<-start
if (!is.na(stop[1]) ) {
#output$Spatial$stop.point<-which.min(sp::spDistsN1(Grid[,1:2], stop, longlat=TRUE))
output$Spatial$stop.point<- output$Spatial$start.point<-sf::st_nearest_feature(sf::st_point(start),Grid_sf, longlat = TRUE)
output$Spatial$stop.location<-stop
} else {
output$Spatial$stop.point<-NA
output$Spatial$stop.location<-NA
}
return(output)
}
Try the FLightR package in your browser
Any scripts or data that you put into this service are public.
FLightR documentation built on Oct. 3, 2023, 1:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions geologger.sampler.create.arrays create.proposal make.processed.light.object make.result.list correct.hours create.calibration make_likelihood_correction_function suggest.irrad.boundaries get.calibration.parameters plot_slopes get.calib.param logger.template.calibration logger.template.calibrarion.internal get.Irradiance process.twilights make.prerun.object find.stationary.location make.calibration plot_slopes_by_location
Documented in find.stationary.location make.calibration make.prerun.object plot_slopes_by_location
# data_preparation functions
#' plots log of observed versus expected slope by time for a known location
#'
#' The function calculates and plots calibration slopes for sunsets and sunrises for every day of the tracking period, based on the assumption that the tag remained in the same (calibration) location all the time.
#'
#' @param Proc.data processed data object generated by \code{\link{get.tags.data}}
#' @param location vector with longitude and latitude of calibration location (degrees).
#' @param log.light.borders numeric vector with length of 2 for minimum and maximum log(light) levels to use. Default value 'auto', will take these values from the Proc.data object.
#' @param log.irrad.borders numeric vector with length of 2 for minimum and maximum log(irradiance) values to use. Default value 'auto', will take these values from the Proc.data object.
#' @param ylim the y limits of the plot. The default value, NULL, indicates that the range of the finite values to be plotted should be used.
#' @param xlim the x limits of the plot. The default value, NULL, otherwise can be POSIXct or character in a form readable by \code{\link{as.POSIXct}}.
#' @details The plot of calibration slopes is used for finding start and end dates of a calibration period (the time period, during which the tag remained in the calibration location with coordinates (x,y)). During the calibration period, the calibration slopes vary little both, between the twilight events (sunrises and sunsets) and in time. When the tag changes location, the slopes for sunrises and sunsets start to deviate. There may potentially be several calibration periods for the same location (if the bird returned to the same location several times). The boundaries (start and end dates) of each of these periods are captured visually. If there were more than one calibration location, the procedure is repeated, once for each location.
#' All the obtained calibration periods can be entered in a data frame 'Calibration.periods', for further analysis. Each line of the data frame contains start and end dates (if applicable) of the calibration period and geographic coordinates of the location.
#' @return 'NULL'
#' @examples
#' File<-system.file("extdata", "Godwit_TAGS_format.csv", package = "FLightR")
#' Proc.data<-get.tags.data(File)
#' plot_slopes_by_location(Proc.data=Proc.data, location=c(5.43, 52.93))
#' abline(v=as.POSIXct("2013-08-20", tz='GMT')) # end of first calibration period
#' abline(v=as.POSIXct("2014-05-05", tz='GMT')) # start of the second calibration period
#'
#' @author Eldar Rakhimberdiev
#' @export plot_slopes_by_location
plot_slopes_by_location<-function(Proc.data, location, log.light.borders='auto', log.irrad.borders='auto', ylim=NULL, xlim=NULL) {
old.par <- graphics::par(no.readonly = TRUE)
Calibration.period<-data.frame(
calibration.start=as.POSIXct("1900-01-01", tz='GMT'),
calibration.stop=as.POSIXct("2050-01-01", tz='GMT'),
lon=location[1], lat=location[2])
if (log.light.borders[1]=='auto') log.light.borders<-Proc.data$log.light.borders
if (log.irrad.borders[1]=='auto') log.irrad.borders<-Proc.data$log.irrad.borders
calibration.parameters<-suppressWarnings(get.calibration.parameters(Calibration.period,
Proc.data, model.ageing=FALSE,
log.light.borders=log.light.borders,
log.irrad.borders=log.irrad.borders,
plot.each = FALSE, plot.final = FALSE, suggest.irrad.borders=FALSE))
tmp<-suppressWarnings(plot_slopes(calibration.parameters$All.slopes, ylim=ylim, xlim=xlim))
#graphics::par(old.par)
return(tmp)
}
#' Creates a calibration object, further used for calculation of coordinates in the \code{\link{run.particle.filter}}.
#'
#' Function estimates all necessary parameters from the calibration data logged in a known location or locations.
#'
#' @param Proc.data processed data object generated by \code{\link{get.tags.data}}
#' @param Calibration.periods a data frame containing start and end dates of all the calibration periods (POSIXct) and geographic coordinates of the corresponding calibration locations.
#' @param model.ageing if set to TRUE, accounts for the tag ageing (with opacification of its transparent shell of a light sensor), resulting into decreasing sensitivity of the device. This option is useful only if there were several calibration periods or if calibration period was very long (~ longer than a month).
#' @param plot.each Do you want every twilight to be plotted while processing
#' @param plot.final Do you want final calibration graph to be plotted. On the graph you can see all the observed versus expected light levels. All slopes should be similar.
#' @param likelihood.correction will estimate correction of likelihood for the current calibration parameters. Highly recommended not to be change from 'auto'. In this case FLightR will switch it to FALSE in case tag saved data on 10 minutes or longer period.
#' @param fixed.logSlope these are mean (1) and SD (2) for distribution of slopes. Should normally be estimated from the data (and thus default is c(NA, NA)). Change any of these two finite values if you want them to be predetermined and not estimated from the calibration data.
#' @param suggest.irrad.borders experimental parameter! If set to TRUE function will try to find the best values for the log.irrad.borders
#' @param return.slopes if true function will return estimated individual twilight slopes.
#' @return calibration object to be uses in the \code{\link{make.prerun.object}}
#' @examples
#' File<-system.file("extdata", "Godwit_TAGS_format.csv", package = "FLightR")
#' Proc.data<-get.tags.data(File, end.date=as.POSIXct('2013-08-20', tz='GMT'))
#' Calibration.periods<-data.frame(
#' calibration.start=NA,
#' calibration.stop=as.POSIXct("2013-08-20", tz='GMT'),
#' lon=5.43, lat=52.93)
#' #use c() also for the geographic coordinates, if you have more than one calibration location
#' # (e. g., lon=c(5.43, 6.00), lat=c(52.93,52.94))
#' print(Calibration.periods)
#'
#' # NB Below likelihood.correction is set to FALSE for fast run!
#' # Leave it as default 'auto' for real examples
#' Calibration<-make.calibration(Proc.data, Calibration.periods, likelihood.correction=FALSE)
#'
#' @author Eldar Rakhimberdiev
#' @export
make.calibration<-function(Proc.data, Calibration.periods, model.ageing=FALSE, plot.each=FALSE, plot.final=FALSE, likelihood.correction='auto', fixed.logSlope=c(NA,NA), suggest.irrad.borders=FALSE, return.slopes=FALSE) {
Calibration.periods[,1]<-as.POSIXct(Calibration.periods[,1], tz='GMT')
Calibration.periods[,2]<-as.POSIXct(Calibration.periods[,2], tz='GMT')
Calibration.periods$calibration.start[is.na(Calibration.periods$calibration.start)]<-as.POSIXct("1900-01-01", tz='GMT')
Calibration.periods$calibration.stop[is.na(Calibration.periods$calibration.stop)]<-as.POSIXct("2100-01-01", tz='GMT')
for (i in 1:nrow(Calibration.periods)) {
if (Calibration.periods[i,1]>=Calibration.periods[i,2]) stop('Calibration start in some period is later than calibration end')
if (i>1) if ( Calibration.periods[i,1]<=Calibration.periods[(i-1),2]) stop('Calibration periods overlap')
}
calibration.parameters<-suppressWarnings(get.calibration.parameters(
Calibration.periods, Proc.data,
model.ageing=model.ageing,
log.light.borders=Proc.data$log.light.borders,
log.irrad.borders=Proc.data$log.irrad.borders,
plot.each= plot.each, plot.final= plot.final,
suggest.irrad.borders=suggest.irrad.borders))
# here I check whether tag measures with 10 minutes or more inteval and if so change likelihood correction to FALSE
if (likelihood.correction=='auto') {
if (Proc.data$saving.period>=550) {
likelihood.correction<-FALSE
message('likelihood correction switched to FALSE\n')
} else {
likelihood.correction<-TRUE
message('likelihood correction switched to TRUE\n')
}
}
Proc.data$log.irrad.borders=calibration.parameters$log.irrad.borders
if (length(calibration.parameters$calib_outliers)>0) {
Proc.data$FLightR.data$twilights$excluded[which(sapply(Proc.data$FLightR.data$twilights$datetime,
FUN=function(x) min(abs(calibration.parameters$calib_outliers-as.numeric(x))))<Proc.data$saving.period*24)]<-1
Proc.data_tmp<-process.twilights(Proc.data$FLightR.data$Data,
Proc.data$FLightR.data$twilights[Proc.data$FLightR.data$twilights$excluded==0,],
measurement.period=Proc.data$measurement.period, saving.period=Proc.data$saving.period,
impute.on.boundaries=Proc.data$impute.on.boundaries)
calibration.parameters<-suppressWarnings(get.calibration.parameters(Calibration.periods, Proc.data_tmp,
model.ageing=model.ageing, log.light.borders=Proc.data$log.light.borders,
log.irrad.borders=Proc.data$log.irrad.borders,
suggest.irrad.borders=suggest.irrad.borders))
plot_slopes(calibration.parameters$All.slopes)
Proc.data$Twilight.time.mat.dusk<-Proc.data_tmp$Twilight.time.mat.dusk
Proc.data$Twilight.time.mat.dawn<-Proc.data_tmp$Twilight.time.mat.dawn
Proc.data$Twilight.log.light.mat.dusk<-Proc.data_tmp$Twilight.log.light.mat.dusk
Proc.data$Twilight.log.light.mat.dawn<-Proc.data_tmp$Twilight.log.light.mat.dawn
}
# second time search for outliers..
if (length(calibration.parameters$calib_outliers)>0) {
Proc.data$FLightR.data$twilights$excluded[which(sapply(Proc.data$FLightR.data$twilights$datetime,
FUN=function(x) min(abs(calibration.parameters$calib_outliers-as.numeric(x))))<Proc.data$saving.period*24)]<-1
Proc.data_tmp<-process.twilights(Proc.data$FLightR.data$Data,
Proc.data$FLightR.data$twilights[Proc.data$FLightR.data$twilights$excluded==0,],
measurement.period=Proc.data$measurement.period, saving.period=Proc.data$saving.period,
impute.on.boundaries=Proc.data$impute.on.boundaries)
calibration.parameters<-suppressWarnings(get.calibration.parameters(
Calibration.periods, Proc.data_tmp,
model.ageing=model.ageing, log.light.borders=Proc.data$log.light.borders,
log.irrad.borders=Proc.data$log.irrad.borders,
suggest.irrad.borders=suggest.irrad.borders))
plot_slopes(calibration.parameters$All.slopes)
Proc.data$Twilight.time.mat.dusk<-Proc.data_tmp$Twilight.time.mat.dusk
Proc.data$Twilight.time.mat.dawn<-Proc.data_tmp$Twilight.time.mat.dawn
Proc.data$Twilight.log.light.mat.dusk<-Proc.data_tmp$Twilight.log.light.mat.dusk
Proc.data$Twilight.log.light.mat.dawn<-Proc.data_tmp$Twilight.log.light.mat.dawn
}
if (!is.na(fixed.logSlope[1])) calibration.parameters$All.slopes$Parameters$LogSlope[1]<-fixed.logSlope[1]
if (!is.na(fixed.logSlope[2])) calibration.parameters$All.slopes$Parameters$LogSlope[2]<-fixed.logSlope[2]
# clean the LogSigma values not to use it under the fixed slope scenario
if (!is.na(fixed.logSlope[1])) calibration.parameters$All.slopes$Parameters$LogSigma<-c(NA,NA)
Calibration<-create.calibration(calibration.parameters$All.slopes,
Proc.data,
Proc.data$FLightR.data,
log.light.borders=Proc.data$log.light.borders, log.irrad.borders=Proc.data$log.irrad.borders,
ageing.model=calibration.parameters$ageing.model,
location=NA, likelihood.correction=likelihood.correction)
Calibration$Calibration.periods<-Calibration.periods
if (return.slopes) Calibration$Calibration.Slopes<-calibration.parameters$All.slopes$Slopes
return(Calibration)
}
#' find unknown calibration location
#'
#' Functions attempts to find a location where
#' The function attempts to find a location for a time period assuming animal was not moving. Does not work well will shaded data!
#'
#' @param Proc.data processed data object generated by \code{\link{get.tags.data}}
#' @param calibration.start POSIXct time when stationary period started
#' @param calibration.stop POSIXct time when stationary period ended
#' @param plot plots every iteration
#' @param initial.coords location vector with initial values for location (longitude and latitude). Should be close (+-2000 km from the real location)
#' @param print.optimization do you want every optimization iteration to be printed? If TRUE - Lon, Lat, calibration mean and calibration sd are being printed. Optimization tries to minimize the latter.
#' @param reltol tolerance for optimization, see \code{\link{optim}} for more details
#' @details The idea behind the function is that it tries to minimize variance between slopes for the whole period by optimizing location. It can be seen as an extension of Hill-Ekstrom calibration idea.
#' @return vector with coordinates - longitude and latitude.
#' @examples
#' #this example takes about 15 minutes to run
#' \dontrun{
#' File<-system.file("extdata", "Godwit_TAGS_format.csv", package = "FLightR")
#' Proc.data<-get.tags.data(File)
#' plot_slopes_by_location(Proc.data=Proc.data, location=c(5.43, 52.93))
#' abline(v=as.POSIXct("2013-08-20", tz='GMT')) # end of first calibration period
#' abline(v=as.POSIXct("2014-05-05", tz='GMT')) # start of the second calibration period
#' Location<-find.stationary.location(Proc.data, '2013-07-20', '2013-08-20', initial.coords=c(10, 50))
#' }
#' @author Eldar Rakhimberdiev
#' @export
find.stationary.location<-function(Proc.data, calibration.start, calibration.stop, plot=TRUE, initial.coords=NULL, print.optimization=TRUE, reltol=1e-4) {
if (is.null(initial.coords)) stop('current function vesrion requires some inital coordinates to start search, they should not be very close but within few thousand km!')
ll_function<-function(initial.coords, Proc.data, calibration.start, calibration.stop, plot=TRUE, stage=1) {
sink("flightr_tmp")
Calibration.period<-data.frame(
calibration.start=as.POSIXct(calibration.start, tz='GMT'),
calibration.stop=as.POSIXct(calibration.stop, tz='GMT'),
lon=initial.coords[1], lat=initial.coords[2])
#use c() also for the geographic coordinates, if you have more than one calibration location
# (e. g., lon=c(5.43, 6.00), lat=c(52.93,52.94))
log.light.borders<-Proc.data$log.light.borders
log.irrad.borders<-Proc.data$log.irrad.borders
calibration.parameters<-invisible(suppressWarnings(get.calibration.parameters(Calibration.period,
Proc.data, model.ageing=FALSE,
log.light.borders=log.light.borders,
log.irrad.borders=log.irrad.borders,
plot.each = FALSE, plot.final = FALSE,
suggest.irrad.borders=FALSE)))
Twilights_total<-calibration.parameters$Twilights_total
if (length(calibration.parameters$calib_outliers)>0) {
Proc.data$FLightR.data$twilights$excluded[which(sapply(Proc.data$FLightR.data$twilights$datetime,
FUN=function(x) min(abs(calibration.parameters$calib_outliers-as.numeric(x))))<Proc.data$saving.period*24)]<-1
Proc.data_tmp<-process.twilights(Proc.data$FLightR.data$Data,
Proc.data$FLightR.data$twilights[Proc.data$FLightR.data$twilights$excluded==0,],
measurement.period=Proc.data$measurement.period, saving.period=Proc.data$saving.period,
impute.on.boundaries=Proc.data$impute.on.boundaries)
calibration.parameters<-suppressWarnings(get.calibration.parameters(Calibration.period,
Proc.data_tmp, model.ageing=FALSE,
log.light.borders=log.light.borders,
log.irrad.borders=log.irrad.borders,
plot.each = FALSE, plot.final = FALSE,
suggest.irrad.borders=FALSE))
}
if (plot) plot_slopes(calibration.parameters$All.slopes)
suppressWarnings(sink())
percent_excluded<-1-(sum(is.finite(calibration.parameters$All.slopes$Slopes$logSlope))/Twilights_total)
if (stage==1) {
if (print.optimization) message(paste(initial.coords[1], initial.coords[2], calibration.parameters$All.slopes$Parameters$LogSlope[1], calibration.parameters$All.slopes$Parameters$LogSlope[2], percent_excluded), '\n')
#print(table(calibration.parameters$All.slopes$Slopes$Type))
if (length(table(calibration.parameters$All.slopes$Slopes$Type))==1) {
warning('only_one_twilight_type_left!\n')
return(sum(is.finite(calibration.parameters$All.slopes$Slopes$logSlope)))
} else {
return(calibration.parameters$All.slopes$Parameters$LogSlope[2]^2+percent_excluded^2)
}
} else {
Dat<-calibration.parameters$All.slopes$Slopes[is.finite(calibration.parameters$All.slopes$Slopes$logSlope),]
if (length(table(Dat$Type))==1 | min(table(Dat$Type))<=2) {
warning('only_one_twilight_type_left!\n')
Val<-nrow(Dat)
} else {
#Dat<-subset(calibration.parameters$All.slopes$Slopes, is.finite(logSlope), select=c(logSlope, Time, Type))
#print(Dat)
Val<-log(1/(stats::anova(stats::lm(logSlope~Time+I(Time^2)+Type, data=Dat),stats::lm(logSlope~Time, data=Dat))[2,6])) +percent_excluded^2 #-log(1/(percent_excluded+0.001)
}
if (print.optimization) message(paste(initial.coords[1], initial.coords[2], calibration.parameters$All.slopes$Parameters$LogSlope[1], calibration.parameters$All.slopes$Parameters$LogSlope[2]), Val, percent_excluded, '\n')
return(Val)
}
}
message('stage 1...\n')
tryCatch(Res<-stats::optim(initial.coords, fn=ll_function, Proc.data=Proc.data, calibration.start=calibration.start, calibration.stop=calibration.stop, plot=plot, control=list(reltol=1e-2)), finally=try(suppressWarnings(sink())))
message('stage 2...\n')
tryCatch(Res<-stats::optim(Res$par, fn=ll_function, Proc.data=Proc.data, calibration.start=calibration.start, calibration.stop=calibration.stop, plot=plot, stage=2,control=list(reltol=reltol)), finally=try(suppressWarnings(sink())))
return(Res$par)
}
#' combines data, calibration and sets up priors
#'
#' This function is one step before \code{\link{run.particle.filter}}. It combines data, calibration, spatial extent and movement priors and estimates spatial likelihoods that used later in the particle filter.
#' @param Proc.data Processed data object created by \code{\link{get.tags.data}}.
#' @param Grid Spatial grid created by \code{\link{make.grid}}.
#' @param start release location (lat, lon).
#' @param end end of the track location. Will use \code{start} by default. Use NA in case of unknown end point.
#' @param Calibration Calibration object created by \code{\link{make.calibration}}.
#' @param threads number of parallel threads to use. default is -1, which means FLightR will use all available threads except 1. Value 1 will force sequential evaluation
#' @param Decision prior for migration probability values from 0 to 1 are allowed
#' @param Direction Direction prior for direction of migration (in degrees) with 0 pointing to the North
#' @param Kappa concentration parameter for vonMises distribution, 0 means uniform or even distribution. Will set some prior for direction for all the track, so is not recommended to be changed
#' @param M.mean Prior for mean distance travelled between consecutive twilights, km
#' @param M.sd Prior for sd of distance travelled between consecutive twilights, the higher the value is the wider is the the distribution
#' @param likelihood.correction Should likelihood correction estimated during \code{\link{make.calibration}} run be used?
#' @return Object to be uses in the \code{\link{run.particle.filter}}
#' @examples
#' File<-system.file("extdata", "Godwit_TAGS_format.csv", package = "FLightR")
#' # to run example fast we will cut the real data file by 2013 Aug 20
#' Proc.data<-get.tags.data(File, end.date=as.POSIXct('2013-07-02', tz='GMT'))
#' Calibration.periods<-data.frame(
#' calibration.start=NA,
#' calibration.stop=as.POSIXct("2013-08-20", tz='GMT'),
#' lon=5.43, lat=52.93)
#' #use c() also for the geographic coordinates, if you have more than one calibration location
#' # (e. g., lon=c(5.43, 6.00), lat=c(52.93,52.94))
#' print(Calibration.periods)
#'
#' # NB Below likelihood.correction is set to FALSE for fast run!
#' # Leave it as default TRUE for real examples
#' Calibration<-make.calibration(Proc.data, Calibration.periods, likelihood.correction=FALSE)
#'
#' Grid<-make.grid(left=0, bottom=50, right=10, top=56,
#' distance.from.land.allowed.to.use=c(-Inf, Inf),
#' distance.from.land.allowed.to.stay=c(-Inf, Inf))
#'
#' all.in<-make.prerun.object(Proc.data, Grid, start=c(5.43, 52.93),
#' Calibration=Calibration, threads=2)
#'
#' @author Eldar Rakhimberdiev
#' @export
make.prerun.object<-function(Proc.data, Grid, start, end=start, Calibration, threads=-1, Decision=0.05, Direction=0, Kappa=0, M.mean=300, M.sd=500, likelihood.correction=TRUE) {
if (length(Decision)>1) stop("Decision has to have length of 1, to sepcify it per twilight, change it in the result object of this function")
if (length(Direction)>1) stop("Direction has to have length of 1, to sepcify it per twilight, change it in the result object of this function")
if (length(Kappa)>1) stop("Kappa has to have length of 1, to sepcify it per twilight, change it in the result object of this function")
if (length(M.mean)>1) stop("M.mean has to have length of 1, to sepcify it per twilight, change it in the result object of this function")
if (length(M.sd)>1) stop("M.sd has to have length of 1, to sepcify it per twilight, change it in the result object of this function")
if (any(is.na(Proc.data$FLightR.data$twilights$datetime))) {
warning('NAs found and deleted in Proc.data$FLightR.data$twilights$datetime\n')
Proc.data$FLightR.data$twilights<- stats::na.omit(Proc.data$FLightR.data$twilights)
}
Processed.light<-make.processed.light.object(Proc.data$FLightR.data)
Index.tab<-create.proposal(Processed.light, start=start, Grid=Grid)
Index.tab$Decision<-Decision # prob of migration
Index.tab$Direction<- Direction # direction 0 - North
Index.tab$Kappa<-Kappa # distr concentration 0 means even
Index.tab$M.mean<- M.mean # distance mu
Index.tab$M.sd<- M.sd # distance sd
all.in<-geologger.sampler.create.arrays(Index.tab, Grid, start=start, stop=end)
all.in$Calibration<-Calibration
all.in$Data<-Proc.data$FLightR.data
if (threads!=1 ) {
Possible.threads<-parallel::detectCores()
if (threads<=0) Threads<-max(Possible.threads+threads, 1)
if (threads>0) Threads<-min(Possible.threads,threads)
} else {
Threads<-1
}
Phys.Mat<-get.Phys.Mat.parallel(all.in, Proc.data$Twilight.time.mat.dusk,
Proc.data$Twilight.log.light.mat.dusk,
Proc.data$Twilight.time.mat.dawn,
Proc.data$Twilight.log.light.mat.dawn,
threads=Threads, calibration=all.in$Calibration, likelihood.correction=likelihood.correction)
all.in$Spatial$Phys.Mat<-Phys.Mat
return(all.in)
}
##############
# functions below were developed before FLightR 0.3.8
# process.data
process.twilights<-function(All.p, Filtered_tw, measurement.period=60, saving.period=NULL, impute.on.boundaries=FALSE) {
# this function just prepares data for the next steps..
##########
## Dusk
# processing Dusk
Dusk.all<-Filtered_tw$datetime[Filtered_tw$type==2]
Twilight.index.mat.dusk<-sapply(which(All.p$gmt %in% Dusk.all & All.p$type==2), FUN=function(x) (x-24):(x+24))
Twilight.index.mat.dusk<-apply(Twilight.index.mat.dusk, c(1,2), FUN=function(x) ifelse (x>0, x, NA))
Max.Index<-nrow(All.p)
Twilight.index.mat.dusk<-apply(Twilight.index.mat.dusk, c(1,2), FUN=function(x) ifelse (x>Max.Index, NA, x))
Twilight.time.mat.dusk<-apply(Twilight.index.mat.dusk, c(1,2), FUN=function(x) as.numeric(All.p$gmt[x]))
Twilight.time.mat.dusk<-apply(Twilight.time.mat.dusk, c(1,2), FUN=function(x) ifelse(is.finite(x), x, 0))
Twilight.log.light.mat.dusk<-apply(Twilight.index.mat.dusk, c(1,2), FUN=function(x) log(All.p$light[x]))
#Twilight.log.light.mat.dusk<-apply(Twilight.index.mat.dusk, c(1,2), FUN=function(x) All.p$light[x])
Twilight.log.light.mat.dusk<-apply(Twilight.log.light.mat.dusk, c(1,2), FUN=function(x) ifelse(is.finite(x), x, -1))
Twilight.time.mat.dusk<-Twilight.time.mat.dusk-(saving.period-measurement.period)
# processing Dawn
Dawn.all<-Filtered_tw$datetime[Filtered_tw$type==1]
Twilight.index.mat.dawn<-sapply(which(All.p$gmt %in% Dawn.all & All.p$type==1), FUN=function(x) (x-24):(x+24))
Twilight.index.mat.dawn<-apply(Twilight.index.mat.dawn, c(1,2), FUN=function(x) ifelse (x>0, x, NA))
Max.Index<-nrow(All.p)
Twilight.index.mat.dawn<-apply(Twilight.index.mat.dawn, c(1,2), FUN=function(x) ifelse (x>Max.Index, NA, x))
Twilight.time.mat.dawn<-apply(Twilight.index.mat.dawn, c(1,2), FUN=function(x) as.numeric(All.p$gmt[x]))
Twilight.time.mat.dawn<-apply(Twilight.time.mat.dawn, c(1,2), FUN=function(x) ifelse(is.finite(x), x, 0))
Twilight.log.light.mat.dawn<-apply(Twilight.index.mat.dawn, c(1,2), FUN=function(x) log(All.p$light[x]))
#Twilight.log.light.mat.dawn<-apply(Twilight.index.mat.dawn, c(1,2), FUN=function(x) All.p$light[x])
Twilight.log.light.mat.dawn<-apply(Twilight.log.light.mat.dawn, c(1,2), FUN=function(x) ifelse(is.finite(x), x, -1))
Res<-list(Twilight.time.mat.dusk=Twilight.time.mat.dusk, Twilight.log.light.mat.dusk=Twilight.log.light.mat.dusk, Twilight.time.mat.dawn=Twilight.time.mat.dawn, Twilight.log.light.mat.dawn=Twilight.log.light.mat.dawn, measurement.period=measurement.period, saving.period=saving.period, impute.on.boundaries=impute.on.boundaries)
return(Res)
}
get.Irradiance<-function(alpha, r=6378, s=6.9, intigeo.template.correction=FALSE) {
# function from Ekstrom 2007
erf <- function(x) 2 * stats::pnorm(x * sqrt(2)) - 1
## (see Abramowitz and Stegun 29.2.29)
## and the so-called 'complementary error function'
erfc <- function(x) 2 * stats::pnorm(x * sqrt(2), lower = FALSE)
## and the inverses
#erfinv <- function (x) qnorm((1 + x)/2)/sqrt(2)
#erfcinv <- function (x) qnorm(x/2, lower = FALSE)/sqrt(2)
Res<-alpha
u<-sqrt(r/(2*s))*sin(alpha)
Res[(u<=0)]<-(exp(-u^2)/(1+erf(-u)))[(u<=0)]
Res[(u>0)]<-(exp(-u^2)/(erfc(u)))[(u>0)]
if (intigeo.template.correction) Res[Res>0.001]=exp(log(Res[Res>0.001])*1.279+log(Res[Res>0.001])^2*0.091)
return(Res)
}
logger.template.calibrarion.internal<-function( Twilight.time.mat.Calib.dawn, Twilight.log.light.mat.Calib.dawn, Twilight.time.mat.Calib.dusk, Twilight.log.light.mat.Calib.dusk, positions=NA, plot.each=TRUE, plot.final=TRUE,log.light.borders=NA, log.irrad.borders=c(-8, 1.3), impute.on.boundaries=FALSE) {
oldpar <- graphics::par(no.readonly = TRUE)
on.exit(graphics::par(oldpar))
# =================
# in this function I'll add a new lnorm calibration...
#
# now I'll try to do them both..
if (!is.list(positions)) {
if (length (positions) !=2) stop("the positions should be a list with $dawn and $dusk dataframes with pairs of coordiantes OR just one pair\n")
dawn=matrix(ncol=2, nrow=dim(Twilight.time.mat.Calib.dawn)[2])
dawn[,1]<-positions[1]
dawn[,2]<-positions[2]
#
dusk=matrix(ncol=2, nrow=dim(Twilight.time.mat.Calib.dusk)[2])
dusk[,1]<-positions[1]
dusk[,2]<-positions[2]
#
positions=list(dawn=dawn, dusk=dusk)
}
Twilight.time.mat.Calib.dawn<-Twilight.time.mat.Calib.dawn[-25,]
Twilight.log.light.mat.Calib.dawn<-Twilight.log.light.mat.Calib.dawn[-25,]
Twilight.time.mat.Calib.dusk<-Twilight.time.mat.Calib.dusk[-25,]
Twilight.log.light.mat.Calib.dusk<-Twilight.log.light.mat.Calib.dusk[-25,]
# let's try to create Calib.data.all first!!
Calib.data.dawn<-data.frame()
if (plot.each) graphics::par(ask=FALSE)
for (dawn in 1:dim(Twilight.time.mat.Calib.dawn)[2]) {
message("\r checking dawn", dawn)
#Twilight.solar.vector<-solar(as.POSIXct(Twilight.time.mat.Calib.dawn[, dawn], tz="gmt", origin="1970-01-01"))
Data<-check.boundaries(positions$dawn[dawn,], Twilight.solar.vector=NULL, Twilight.log.light.vector = Twilight.log.light.mat.Calib.dawn[,dawn], plot=plot.each, verbose=FALSE, log.light.borders=log.light.borders, log.irrad.borders=log.irrad.borders, dusk=FALSE, Twilight.time.vector=Twilight.time.mat.Calib.dawn[, dawn], impute.on.boundaries=impute.on.boundaries)
if (length(Data)==0) {
message("\r dawn", dawn, "was excluded from the calibration")
} else {
Calib.data.dawn<-rbind(Calib.data.dawn, cbind(LogLight=Data[,1], LogIrrad=Data[,2], Day=dawn, Time=Data[,3], Elevs=Data[,4]))
}
}
if (nrow(Calib.data.dawn)>0) {Calib.data.dawn$type<-"Dawn"}
Calib.data.dusk<-data.frame()
for (dusk in 1:dim(Twilight.time.mat.Calib.dusk)[2]) {
message("\r checking dusk", dusk )
#Twilight.solar.vector<-solar(as.POSIXct(Twilight.time.mat.Calib.dusk[, dusk], tz="gmt", origin="1970-01-01"))
Data<-check.boundaries(positions$dusk[dusk,], Twilight.solar.vector=NULL, Twilight.log.light.vector=Twilight.log.light.mat.Calib.dusk[,dusk], plot=plot.each, verbose=FALSE, log.light.borders=log.light.borders, log.irrad.borders=log.irrad.borders, dusk=TRUE, Twilight.time.vector=Twilight.time.mat.Calib.dusk[, dusk], impute.on.boundaries=impute.on.boundaries)
if (length(Data)==0) {
message("\r dusk", dusk, "was excluded from the calibration")
} else {
Calib.data.dusk<-rbind(Calib.data.dusk, cbind(LogLight=Data[,1], LogIrrad=Data[,2], Day=dim(Twilight.time.mat.Calib.dawn)[2]+dusk, Time=Data[,3], Elevs=Data[,4]))
}
}
if (nrow(Calib.data.dusk)>0) {Calib.data.dusk$type<-"Dusk"}
message('\r')
if (nrow(Calib.data.dusk)==0 & nrow(Calib.data.dawn)==0) stop('Error: None of twilights could happen in the location... Something went wrong at the previous stages')
Calib.data.all<-rbind(Calib.data.dawn, Calib.data.dusk)
if (plot.final) {
graphics::plot(LogLight~LogIrrad,data=Calib.data.all, type="n")
for (Day in unique(Calib.data.all$Day)) {
graphics::lines(LogLight~LogIrrad, data=Calib.data.all[Calib.data.all$Day==Day,], type="b", pch="+", col=grDevices::rainbow(length(unique(Calib.data.all$Day)))[Day])
}
}
Calib.data.all$fDay<-as.factor(Calib.data.all$Day)
message('\r\n')
message('\r\n')
return(Calib.data.all)
}
# ok now we need template.calibration.function..
logger.template.calibration<-function(Twilight.time.mat.Calib.dawn, Twilight.log.light.mat.Calib.dawn, Twilight.time.mat.Calib.dusk, Twilight.log.light.mat.Calib.dusk, time.shift=0, positions, log.light.borders=NA, log.irrad.borders=c(-15, 50), plot.each=TRUE, plot.final=TRUE, impute.on.boundaries=FALSE) {
Calibration.original<-logger.template.calibrarion.internal(Twilight.time.mat.Calib.dawn+time.shift, Twilight.log.light.mat.Calib.dawn, Twilight.time.mat.Calib.dusk+time.shift, Twilight.log.light.mat.Calib.dusk, positions=positions, log.light.borders=log.light.borders, log.irrad.borders= log.irrad.borders, plot.each=plot.each, plot.final=plot.final,impute.on.boundaries=impute.on.boundaries)
return(Calibration.original)
}
get.calib.param<-function(Calib.data.all, plot=FALSE, calibration.type=NULL) {
if (is.null(calibration.type)) calibration.type="parametric.slope"
message("calibration method used:", calibration.type, "\n")
Calib.data.all$fTwilight<-Calib.data.all$fDay
#============================
# we should do some outlier tests...
cur.data<-Calib.data.all
p.lm1<-stats::lm(LogLight~fTwilight/LogIrrad,data=cur.data)
cur.slope<-matrix(stats::coef(p.lm1),ncol=2) # foo[,2] - slope
cur.slope<-as.data.frame(cur.slope)
Diag<-diag(stats::vcov(p.lm1))
Diag.slopes<-Diag[which(sapply(strsplit(names(Diag), ":"), length)==2)]
cur.slope$sd[!is.na(cur.slope[,2])]<-sqrt(Diag.slopes)
cur.slope$type<-stats::aggregate(cur.data[,"type"],by=list(Day=cur.data$fTwilight),FUN=function(x) unique(x))[,2]
cur.slope$time<-stats::aggregate(cur.data[,"Time"],by=list(Day=cur.data$fTwilight),FUN=function(x) x[1])[,2]
cur.slope$Duration<-stats::aggregate(cur.data[,"Time"],by=list(Day=cur.data$fTwilight),FUN=function(x) max(x)-min(x))[,2]
names(cur.slope)[2]<-"slope"
#===========================
# ok I've decided to make separate lm
Slopes<-c()
Slopes.sd<-c()
Intercept=c()
Sigma=c()
Type=c()
Time=c()
Elevs<-c()
Day<-c()
for (i in (unique(cur.data$fTwilight))) {
# lm
#plot(LogLight~LogIrrad, data=cur.data[cur.data$fTwilight==i,])
Data<-cur.data[cur.data$fTwilight==i,]
Data<-Data[Data$LogLight>0,]
Lm<-stats::lm(LogLight~LogIrrad,data=Data)
Cur_slope_mean<-stats::coef(Lm)[2]
Cur_slope_sd<-sqrt(stats::vcov(Lm)[4])
if (calibration.type=="nonparametric.slope") {
# Here I want to physically estimate slopes
All_slopes_cur<-diff(Data$LogLight)/diff(Data$LogIrrad)
All_slopes_mean<-mean(All_slopes_cur)
All_slopes_sd<-stats::sd(All_slopes_cur)
# exclude outliers
All_slopes_diff<-abs(All_slopes_cur-All_slopes_mean)
Slopes_outliers<-which(All_slopes_diff>3*All_slopes_sd)
if (length(Slopes_outliers)>0) {
All_slopes_cur<-All_slopes_cur[-Slopes_outliers]
All_slopes_mean<-mean(All_slopes_cur)
All_slopes_sd<-stats::sd(All_slopes_cur)
}
Cur_slope_mean<-All_slopes_mean
Cur_slope_sd<-All_slopes_sd
}
Slopes<-c(Slopes, Cur_slope_mean)
Slopes.sd<-c(Slopes.sd, Cur_slope_sd)
Intercept=c(Intercept, stats::coef(Lm)[1])
Sigma<-c(Sigma, summary(Lm)$sigma)
Type=c(Type, cur.data$type[cur.data$fTwilight==i][1])
Time<-c(Time,ifelse(cur.data$type[cur.data$fTwilight==i][1]=="Dusk", max(cur.data$Time[cur.data$fTwilight==i]), min(cur.data$Time[cur.data$fTwilight==i])))
Elevs=c(Elevs, mean(cur.data$Elevs[cur.data$fTwilight==i], na.rm=TRUE))
Day=c(Day, cur.data$Day[cur.data$fTwilight==i][1])
# nonparametric slope estimation
}
#plot(cur.slope$slope~Slopes)
#plot(cur.slope$sd~Slopes.sd)
# ok, we now want to replace old slopes with new ones.
# but first we want to add 0.5 instead of NA - at least this is how it works now inside the functions and how it was done in a simulations
cur.slope$slope<-Slopes
cur.slope$Intercept<-Intercept
cur.slope$Sigma<-Sigma
#Slopes.sd[is.na(Slopes.sd)]<-sd.fast
cur.slope$sd<-Slopes.sd
names(cur.slope)[2]<-"slope"
#=====================
if (plot) graphics::hist(log(cur.slope$slope))
cur.slope.int<-cur.slope[is.finite(log(cur.slope$slope)),]
Parameters<-list(Intercept=c(mean(cur.slope.int$Intercept, na.rm=TRUE), stats::sd(cur.slope.int$Intercept, na.rm=TRUE)), LogSlope=c(mean(log(cur.slope.int$slope), na.rm=TRUE), stats::sd(log(cur.slope.int$slope), na.rm=TRUE)), LogSigma=c(mean(log(cur.slope.int$Sigma[!is.na(cur.slope.int$sd)])), stats::sd(log(cur.slope.int$Sigma[!is.na(cur.slope.int$sd)]))), mean.of.individual.slope.sigma=mean(cur.slope.int$sd, na.rm=TRUE), calibration.type=calibration.type)
#cur.slope$time<-aggregate(cur.data[,"Time"],by=list(Day=cur.data$fTwilight),FUN=function(x) x[1])[,2]
#cur.slope$time<-aggregate(cur.data[,"Time"],by=list(Day=cur.data$fTwilight),FUN=mean)[,2]
Res<-list(Parameters=Parameters, Slopes=data.frame(Slope=cur.slope$slope, Time=Time, Intercept=cur.slope$Intercept, Sigma=cur.slope$Sigma, Slopes.sd=Slopes.sd, Type=Type, Elevs=Elevs, Day=Day))
return(Res)
}
plot_slopes<-function(all.slopes, ylim=NULL, xlim=NULL) {
oldpar <- graphics::par(no.readonly = TRUE)
all.slopes$Slopes<-all.slopes$Slopes[all.slopes$Slopes$Slope>0,]
if (is.null(xlim)) {
graphics::plot(log(all.slopes$Slopes$Slope)~as.POSIXct(all.slopes$Slopes$Time, tz="GMT", origin="1970-01-01"), type="n", main="red - dawn, black - dusk", xlab="time", ylab="log(Slope)", ylim=ylim, las=1)
} else {
graphics::plot(log(all.slopes$Slopes$Slope)~ as.POSIXct(all.slopes$Slopes$Time, tz="GMT", origin="1970-01-01"), type="n", main="red - dawn, black - dusk", xlab="time", ylab="log(Slope)", ylim=ylim, las=1, xlim= as.POSIXct(xlim, tz='GMT'))
}
graphics::lines(log(Slope)~ as.POSIXct(Time, tz="GMT", origin="1970-01-01"), data=all.slopes$Slopes[all.slopes$Slopes$Type=="Dusk",])
graphics::points(log(Slope)~ as.POSIXct(Time, tz="GMT", origin="1970-01-01"), data=all.slopes$Slopes[all.slopes$Slopes$Type=="Dusk",], pch="+")
graphics::points(log(Slope)~ as.POSIXct(Time, tz="GMT", origin="1970-01-01"), data=all.slopes$Slopes[all.slopes$Slopes$Type=="Dawn",], pch="+", col="red")
graphics::lines(log(Slope)~ as.POSIXct(Time, tz="GMT", origin="1970-01-01"), data=all.slopes$Slopes[all.slopes$Slopes$Type=="Dawn",], col="red")
oldpar$usr<-graphics::par()$usr
on.exit(graphics::par(oldpar))
#invisible()
return(NULL)
}
get.calibration.parameters<-function(Calibration.periods, Proc.data, model.ageing=TRUE, log.light.borders=log(c(2, 63)), log.irrad.borders=c(-7,1.5), plot.each=FALSE, plot.final=TRUE, calibration.type=NULL, suggest.irrad.borders=FALSE) {
# ageing.model this option should be used only in case there are two periods of calibration or there is one but very long.
if (nrow(Calibration.periods)==1 & model.ageing) warning("you have only one calibration period ageing estimation is unreliable and should be turned to F!!!")
Calibration.periods.int<-Calibration.periods
Calibration.periods.int$calibration.start<-as.numeric(Calibration.periods.int$calibration.start)
Calibration.periods.int$calibration.stop<-as.numeric(Calibration.periods.int$calibration.stop)
Dusk.calib.days<-c()
Dawn.calib.days<-c()
Positions=list(dawn=c(), dusk=c())
#sum(apply(Calibration.periods.int, 1, FUN=function(x) x[2]-x[1]))/3600/24
for (i in 1:nrow(Calibration.periods)) {
Dusk.calib.days.cur<-NULL
Dusk.calib.days.cur<-which(Proc.data$Twilight.time.mat.dusk[25,]>Calibration.periods.int$calibration.start[i] & Proc.data$Twilight.time.mat.dusk[25,] < Calibration.periods.int$calibration.stop[i])
if (length(Dusk.calib.days.cur) > 0) {
Dusk.calib.days<-unique(c(Dusk.calib.days, Dusk.calib.days.cur))
dusk.cur=matrix(ncol=2, nrow=length(Dusk.calib.days.cur))
dusk.cur[,1]<-Calibration.periods.int$lon[i]
dusk.cur[,2]<-Calibration.periods.int$lat[i]
Positions$dusk<-rbind(Positions$dusk, dusk.cur)
}
Dawn.calib.days.cur<-NULL
Dawn.calib.days.cur<-which(Proc.data$Twilight.time.mat.dawn[25,]>Calibration.periods.int$calibration.start[i] & Proc.data$Twilight.time.mat.dawn[25,] < Calibration.periods.int$calibration.stop[i] )
if (length(Dawn.calib.days.cur) > 0) {
Dawn.calib.days<-unique(c(Dawn.calib.days,Dawn.calib.days.cur ))
dawn.cur=matrix(ncol=2, nrow=length(Dawn.calib.days.cur))
dawn.cur[,1]<-Calibration.periods.int$lon[i]
dawn.cur[,2]<-Calibration.periods.int$lat[i]
Positions$dawn<-rbind(Positions$dawn, dawn.cur)
}
}
Twilight.time.mat.Calib.dusk<-Proc.data$Twilight.time.mat.dusk[,Dusk.calib.days]
Twilight.log.light.mat.Calib.dusk<-Proc.data$Twilight.log.light.mat.dusk[,Dusk.calib.days]
## Dawn
Twilight.time.mat.Calib.dawn<-Proc.data$Twilight.time.mat.dawn[,Dawn.calib.days]
Twilight.log.light.mat.Calib.dawn<-Proc.data$Twilight.log.light.mat.dawn[,Dawn.calib.days ]
Twilights_total<-length(c(Dusk.calib.days,Dawn.calib.days))
#-------------------#
if (suggest.irrad.borders) {
Calib.data.all<-logger.template.calibration(Twilight.time.mat.Calib.dawn, Twilight.log.light.mat.Calib.dawn, Twilight.time.mat.Calib.dusk, Twilight.log.light.mat.Calib.dusk, positions=Positions, log.light.borders=log.light.borders, log.irrad.borders=c(-10,10), plot.each=plot.each, plot.final=plot.final, impute.on.boundaries=Proc.data$impute.on.boundaries)
log.irrad.borders<-suggest.irrad.boundaries(Calib.data.all)
warning('!!! log.irrad.borders were updated to new values:', log.irrad.borders[1], log.irrad.borders[2], '\n')
}
Calib.data.all<-logger.template.calibration(Twilight.time.mat.Calib.dawn, Twilight.log.light.mat.Calib.dawn, Twilight.time.mat.Calib.dusk, Twilight.log.light.mat.Calib.dusk, positions=Positions, log.light.borders=log.light.borders, log.irrad.borders=log.irrad.borders, plot.each=plot.each, plot.final=plot.final, impute.on.boundaries=Proc.data$impute.on.boundaries)
All.slopes<-suppressWarnings(get.calib.param(Calib.data.all, plot=FALSE, calibration.type=calibration.type))
All.slopes$Slopes$logSlope<-suppressWarnings(log(All.slopes$Slopes$Slope))
if (model.ageing) {
All.slopes.int<-All.slopes
All.slopes.int$Slopes<-All.slopes.int$Slopes[is.finite(All.slopes.int$Slopes$logSlope),]
Time.start=min(c(Twilight.time.mat.Calib.dusk, Twilight.time.mat.Calib.dawn))
Model=stats::lm(logSlope~I(Time-Time.start), data=All.slopes.int$Slopes)
Model$Time.start<-Time.start
calib_outliers<-All.slopes.int$Slopes$Time[which(abs(stats::residuals(Model))>3*stats::sd(stats::residuals(Model))*sqrt(length(stats::residuals(Model))))]
Res<-list(calib_outliers=calib_outliers, ageing.model=Model, All.slopes=All.slopes,log.irrad.borders=log.irrad.borders)
} else {
All.slopes.int<-All.slopes
All.slopes.int$Slopes<-All.slopes.int$Slopes[is.finite(All.slopes.int$Slopes$logSlope),]
calib_outliers<-All.slopes.int$Slopes$Time[which(abs(All.slopes.int$Slopes$logSlope-mean(All.slopes.int$Slopes$logSlope, na.rm=TRUE))>3*stats::sd(All.slopes.int$Slopes$logSlope, na.rm=TRUE))]
Res<-list(calib_outliers=calib_outliers, All.slopes=All.slopes, log.irrad.borders=log.irrad.borders, Twilights_total=Twilights_total)
}
return(Res)
}
suggest.irrad.boundaries<-function(Calib.data.all, leave.out=0.01) {
Model<-stats::lm(LogLight~LogIrrad + fDay, data=Calib.data.all)
MinLight<-round(min(Calib.data.all$LogLight),1)
Left.border<- -(stats::coef(Model)[1] + stats::quantile(stats::coef(Model)[3:length(stats::coef(Model))], leave.out) - MinLight)*stats::coef(Model)[2]
MaxLight<-round(max(Calib.data.all$LogLight),1)
Right.border<-(MaxLight- (stats::coef(Model)[1] + stats::quantile(stats::coef(Model)[3:length(stats::coef(Model))], 1-leave.out))) *stats::coef(Model)[2]
return(c(Left.border, Right.border))
}
make_likelihood_correction_function<-function(calib_log_mean, calib_log_sd, cur_sd_range, npoints=300, plot=FALSE) {
cur_mean_range<-c(exp(calib_log_mean)-5*exp(calib_log_sd), exp(calib_log_mean)+5*exp(calib_log_sd))
Res<-c()
for (i in 1:npoints) {
message('\r simulation: ', round(i/npoints*100), '%', sep='')
message('\r')
cur_mean<-stats::runif(300, cur_mean_range[1], cur_mean_range[2])
cur_sd<-stats::runif(1, cur_sd_range[1], cur_sd_range[2])
Cur_max_real<- stats::optimize(f=function(x)
mean(stats::dlnorm(stats::rnorm(200000,x , cur_sd),calib_log_mean, calib_log_sd)),
interval=cur_mean_range, maximum = TRUE)$maximum
Res<-rbind(Res, cbind(calib_log_mean, calib_log_sd, cur_sd, cur_mean_max=Cur_max_real))
}
Res<-as.data.frame(Res)
Res$Corr <- exp(Res$calib_log_mean)-Res$cur_mean_max
message('\r estimating correction function...')
message('\r')
MvExp<-mgcv::gamm(Corr~s(cur_sd), data=Res, weights=nlme::varExp(form =~ cur_sd))
# now we check for the outliers..
Resid<-stats::resid(MvExp$lme, type='normalized')
Outliers<-NULL
Outliers<-which(abs(Resid)>3*stats::sd(Resid))
if (length(Outliers)>0) {
Res<-Res[-Outliers,]
MvExp<-mgcv::gamm(Corr~s(cur_sd), data=Res, weights=nlme::varExp(form =~ cur_sd))
}
XX<-seq(cur_sd_range[1], cur_sd_range[2], by=0.01)
c_fun<-stats::approxfun(mgcv::predict.gam(MvExp$gam, newdata=data.frame(cur_sd=XX)) ~XX)
if (plot) {
graphics::plot(Res[,4]~Res[,3], pch='+', ylab='cur_mean', xlab='cur_sd')
graphics::lines(exp(Res$calib_log_mean[1])-mgcv::predict.gam(MvExp$gam, newdata=data.frame(cur_sd=XX))~XX, col='red')
}
Out<-list(c_fun=c_fun, Res=Res)
message('\r')
message('\n')
return(Out)
}
create.calibration<-function( All.slopes, Proc.data, FLightR.data, location, log.light.borders, log.irrad.borders, ageing.model=NULL, likelihood.correction=TRUE) {
# Now we create 'parameters' object that will have all the details about the calibration
Parameters<-All.slopes$Parameters # LogSlope #
Parameters$measurement.period<-Proc.data$measurement.period
Parameters$saving.period<-Proc.data$saving.period
Parameters$log.light.borders<-log.light.borders # these are the boundaries in which one should use the BAS tag.. for the other types of tags they will be different.
Parameters$min.max.values<-c(min(FLightR.data$Data$light), max(FLightR.data$Data$light))
Parameters$log.irrad.borders=log.irrad.borders
Parameters$impute.on.boundaries=Proc.data$impute.on.boundaries
Parameters$location=location
Parameters$LogSlope_1_minute<-Parameters$LogSlope
if (is.null(ageing.model)) {
time_correction_fun= eval(parse(text=paste("function (x,y) return(", Parameters$LogSlope[1], ")")))
lat_correction_fun<-function(x, y, z) return(0)
} else {
time_correction_fun= function(x, y) return(0)
lat_correction_fun<-eval(parse(text=paste("function (x,y) return(", stats::coef(ageing.model)[1], "+", stats::coef(ageing.model)[2], "* (y-",ageing.model$Time.start," ))")))
}
c_fun=NULL
if (likelihood.correction) {
if (is.na(All.slopes$Parameters$LogSigma[1])) {
cur_sd_range=c(0, 1)
} else {
cur_sd_range=c(0, max(exp(All.slopes$Parameters$LogSigma[1]+5*All.slopes$Parameters$LogSigma[2]),1))
}
c_fun=make_likelihood_correction_function(Parameters$LogSlope[1], Parameters$LogSlope[2], cur_sd_range=cur_sd_range)$c_fun
}
Calibration<-list(Parameters=Parameters, time_correction_fun=time_correction_fun, lat_correction_fun=lat_correction_fun, c_fun=c_fun)
return(Calibration)
}
correct.hours<-function(datetime) {
# this function is supposed to correct hours by adding some amount of them.
hours <- as.numeric(format(datetime,"%H"))+as.numeric(format(datetime,"%M"))/60
hours[as.POSIXlt(datetime, tz='GMT')$isdst==1] <-hours[as.POSIXlt(datetime, tz='GMT')$isdst==1]-1
cor <- rep(NA, 24)
for(i in 0:23){
cor[i+1] <- max(abs((c(hours[1],hours)+i)%%24 -
(c(hours,hours[length(hours)])+i)%%24),na.rm=TRUE)
}
hours <- (hours + (which.min(round(cor,2)))-1)%%24
return(hours)
}
make.result.list<-function(Data, raw.X, raw.Y) {
Int<-c(min(raw.Y), max(raw.Y))
Index<-sapply(raw.X, FUN=function(x) which.min(abs(as.numeric(Data$d$gmt-x))))
Res<-Data$d[Index,]
# now I want to adjust time and light
Res$light<-stats::approx(x=Data$d$gmt, y=Data$d$light, xout=raw.X)$y
Res$gmt<-raw.X
Res$Hour<-raw.Y
Result<-list(Data=Res)
Result$Data$gmt.adj<-Result$Data$gmt
Result$Data$gmt<- as.POSIXct(Result$Data$gmt, tz="GMT", origin="1970-01-01")
Result$Data$gmt.adj<- as.POSIXct(Result$Data$gmt.adj, tz="GMT", origin="1970-01-01")
return(Result)
}
make.processed.light.object<-function(FLightR.data) {
# dusk
raw.Y.dusk<-correct.hours(FLightR.data$twilights$datetime[FLightR.data$twilights$type==2 & FLightR.data$twilights$excluded==0])
raw.X.dusk<-as.numeric(FLightR.data$twilights$datetime[FLightR.data$twilights$type==2 & FLightR.data$twilights$excluded==0])
Data_tmp<-list(d=FLightR.data$Data)
Result.Dusk<-make.result.list(Data_tmp, raw.X.dusk, raw.Y.dusk)
# dusk
raw.Y.dawn<-correct.hours(FLightR.data$twilights$datetime[FLightR.data$twilights$type==1 & FLightR.data$twilights$excluded==0])
raw.X.dawn<-as.numeric(FLightR.data$twilights$datetime[FLightR.data$twilights$type==1 & FLightR.data$twilights$excluded==0])
Result.Dawn<-make.result.list(Data_tmp, raw.X.dawn, raw.Y.dawn)
# combine
processed.light<-list(Final.dusk=Result.Dusk, Final.dawn=Result.Dawn)
return(processed.light)
}
# this create proposal is corrected for the missing loess
create.proposal<-function(processed.light, start=c(-98.7, 34.7), end=NA, Grid) {
if (is.na(end)) end=start
All.Days<-seq(min(min(processed.light$Final.dusk$Data$gmt), min(processed.light$Final.dawn$Data$gmt)),max(max(processed.light$Final.dusk$Data$gmt), max(processed.light$Final.dawn$Data$gmt)), by="days")
## Probabilty.of.migration
All.Days.extended<-seq(min(All.Days)-86400, max(All.Days)+86400, by="days")
# ver 0.4.2
# correct twilights for polar periods..
start_no_polar<-start
start_no_polar[2]<-min(abs(start_no_polar[2]), 65.73) * sign(start_no_polar[2])
# these are all potential twilights that we would have for the start point...
Potential.twilights<-sort(c(suntools::sunriset(matrix(start_no_polar, nrow=1), All.Days.extended, direction="sunrise", POSIXct.out=TRUE)[,2], suntools::sunriset(matrix(start_no_polar, nrow=1), All.Days.extended, direction="sunset", POSIXct.out=TRUE)[,2]))
# now we need to add dask or dawn to it..
Sunrise<-suntools::sunriset(matrix(start_no_polar, nrow=1), Potential.twilights[1], direction="sunrise", POSIXct.out=TRUE)[,2]
Sunrises<-c(Sunrise-(3600*24), Sunrise, Sunrise+(3600*24))
Sunset<-suntools::sunriset(matrix(start_no_polar, nrow=1), Potential.twilights[1], direction="sunset", POSIXct.out=TRUE)[,2]
Sunsets<-c(Sunset-(3600*24), Sunset, Sunset+(3600*24))
First.twilight<-ifelse(which.min(c(min(abs(difftime(Sunrises,Potential.twilights[1], units="mins"))), min(abs(difftime(Sunsets,Potential.twilights[1], units="mins")))))==1, "dawn", "dusk")
Index.tab<-data.frame(Date=Potential.twilights)
if (First.twilight=="dusk") {
Index.tab$Dusk<-rep(c(TRUE,FALSE), times=length(All.Days.extended))
} else {
Index.tab$Dusk<-rep(c(FALSE,TRUE), times=length(All.Days.extended))
}
Index.tab$Curr.mat<-NA # this will be NA if no data and row number if there are..
Index.tab$Real.time<- as.POSIXct(NA, tz="GMT")
Index.tab$time<- as.POSIXct(NA, tz="GMT")
######################################
# !!!! this will not work if bird will move for over 12 time zones!!!
# Dusk
for (i in 1:length(processed.light$Final.dusk$Data$gmt)) {
Row2write<-which.min(abs(difftime(processed.light$Final.dusk$Data$gmt[i], Index.tab$Date[Index.tab$Dusk==TRUE], units="mins")))
Index.tab$time[Index.tab$Dusk==TRUE][Row2write]<-processed.light$Final.dusk$Data$gmt[i]
Index.tab$Real.time[Index.tab$Dusk==TRUE][Row2write]<-processed.light$Final.dusk$Data$gmt.adj[i]
Index.tab$Curr.mat[Index.tab$Dusk==TRUE][Row2write]<-i
}
# Dawn
for (i in 1:length(processed.light$Final.dawn$Data$gmt)) {
Row2write<-which.min(abs(difftime(processed.light$Final.dawn$Data$gmt[i], Index.tab$Date[Index.tab$Dusk==FALSE], units="mins")))
Index.tab$time[Index.tab$Dusk==FALSE][Row2write]<-processed.light$Final.dawn$Data$gmt[i]
Index.tab$Real.time[Index.tab$Dusk==FALSE][Row2write]<-processed.light$Final.dawn$Data$gmt.adj[i]
Index.tab$Curr.mat[Index.tab$Dusk==FALSE][Row2write]<-i
}
# cutting empty ends..
while (is.na(Index.tab$Curr.mat[1])) Index.tab<-Index.tab[-1,]
while (is.na(Index.tab$Curr.mat[nrow(Index.tab)])) Index.tab<-Index.tab[-nrow(Index.tab),]
Index.tab$Point<-NA
#First.Point<-which.min(sp::spDistsN1(Grid[,1:2], start, longlat=TRUE))
Grid_sf<-Grid[,1:2] |> sf::st_multipoint() |> sf::st_sfc() |> sf::st_cast('POINT')
First.Point<-sf::st_nearest_feature(sf::st_point(start),Grid_sf, longlat = TRUE)
Index.tab$Point[1]<-First.Point
#
# I decided that Curr.mat is not needed anymore
#Index.tab$Main.Index[,-which(names(Index.tab$Main.Index)=="Curr.mat")]
# this need to double checked
Index.tab$yday<-as.POSIXlt(Index.tab$Date, tz="GMT")$yday
return(Index.tab)
}
geologger.sampler.create.arrays<-function(Index.tab, Grid, start, stop=start) {
# this function wil take in Index.table with all the proposals and will create an array for case of missed data..
# the main feature is that it can account for missing data..
Index.tab.old<-Index.tab
Index.tab$proposal.index<-NA
Index.tab$proposal.index[Index.tab$Dusk==TRUE]<-"Dusk"
Index.tab$proposal.index[Index.tab$Dusk==FALSE]<-"Dawn"
Missed.twilights<-which(is.na(Index.tab$Curr.mat))
if (length(Missed.twilights)>0) {
Deleted.Rows.Count<-0
for (i in Missed.twilights) {
i=i-Deleted.Rows.Count
Index.tab$proposal.index[i-1]<-"Comb"
Index.tab$Decision[i-1]<-(Index.tab$Decision[i-1]+Index.tab$Decision[i])/2 #
if (Index.tab$Direction[i-1] != Index.tab$Direction[i]) {
Index.tab$Direction[i-1]<-0
Index.tab$Kappa[i-1]<-0
}
Index.tab$M.mean[i-1]<-Index.tab$M.mean[i-1]+Index.tab$M.mean[i]
Index.tab$M.sd[i-1]<-sqrt((Index.tab$M.sd[i-1])^2+(Index.tab$M.sd[i])^2)
Index.tab<-Index.tab[-i,]
Deleted.Rows.Count=Deleted.Rows.Count+1
}
}
output<-list()
Matrix.Index.Table<- Index.tab[,which(names(Index.tab) %in% c("Decision", "Direction", "Kappa", "M.mean", "M.sd", "yday", "Real.time", "time", "Dusk", "Loess.se.fit", "Loess.n")) ]
# I also remove last line as bird was not flying after last twilight
Matrix.Index.Table<-Matrix.Index.Table[-nrow(Matrix.Index.Table),]
# main index will have the same amount of rows we have in twilight matrices without first..
Main.Index<-c()
Main.Index$Biol.Prev<-1:nrow(Matrix.Index.Table)
Main.Index$proposal.index<-Index.tab$proposal.index[-nrow(Index.tab)]
Main.Index<-as.data.frame(Main.Index)
Indices<-list(Matrix.Index.Table=Matrix.Index.Table, Main.Index=Main.Index)
output$Indices=Indices
# ok now we need to add to output all other parts..
output$Spatial<-list(Grid=Grid) #[,1:2]
output$Spatial$Behav.mask<-as.integer(Grid[,3])
#output$Spatial$start.point<-which.min(sp::spDistsN1(Grid[,1:2], start, longlat=TRUE))
Grid_sf<-Grid[,1:2] |> sf::st_multipoint() |> sf::st_sfc() |> sf::st_cast('POINT')
output$Spatial$start.point<-sf::st_nearest_feature(sf::st_point(start),Grid_sf, longlat = TRUE)
output$Spatial$start.location<-start
if (!is.na(stop[1]) ) {
#output$Spatial$stop.point<-which.min(sp::spDistsN1(Grid[,1:2], stop, longlat=TRUE))
output$Spatial$stop.point<- output$Spatial$start.point<-sf::st_nearest_feature(sf::st_point(start),Grid_sf, longlat = TRUE)
output$Spatial$stop.location<-stop
} else {
output$Spatial$stop.point<-NA
output$Spatial$stop.location<-NA
}
return(output)
}
Try the FLightR package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.