Skip to content
Benjamin Merkel edited this page Dec 29, 2016 · 17 revisions

How to use the probGLS package with different logger brands

made by Benjamin Merkel, August 2016

---------------------- download environmental data -------

download yearly NetCDF files for (replace YEAR with appropriate number):

  • daily mean SST -> 'sst.day.mean.YEAR.v2.nc'
  • daily SST error -> 'sst.day.err.YEAR.v2.nc'
  • daily mean sea ice concentration -> 'icec.day.mean.YEAR.v2.nc'

from: http://www.esrl.noaa.gov/psd/data/gridded/data.noaa.oisst.v2.highres.html and place all into the same folder

Also, download the land mask file: 'lsmask.oisst.v2.nc' from the same directory and place it in the same folder as all the other NetCDF files

---------------------- load package ----------------------

this package depends on the packages: GeoLight, raster, sp, geosphere, ncdf4 and maps

library(probGLS) 

---------------------- BRAND Lotek -----------------------

load twilight events

lotek <- read.csv('daylogfile')
trn   <- lotek_to_dataframe(date    = as.Date(strptime(lotek$TimeS,"%d/%m/%y")),
                            sunrise = as.character(lotek$Sunrise),
                            sunset  = as.character(lotek$Sunset),
                            TRLat   = lotek$TRLat,
                            TRLon   = lotek$TRLon)

load activity file

act           <- read.csv('basiclogfile')
act$dtime     <- as.POSIXct(strptime(act[,1], format = "%H:%M:%S %d/%m/%y"),tz='UTC') # datetime object needs to be POSIXct in UTC and must be called 'dtime'
act$wetdry    <- act$WetDryState # wet dry data column must be called 'wetdry'
act           <- act[!is.na(act$dtime),]

load temperature file

td            <- read.csv('basiclogfile')
td$dtime      <- as.POSIXct(strptime(td[,1], format = "%H:%M:%S %d/%m/%y"),tz='UTC') # datetime object needs to be POSIXct in UTC and must be called 'dtime'

This step depends on the set up of the tag. Here we assume that the tag logged internal temperature every 10 min regardless of the wetdry state

td$WetDryState.before <- c(NA,head(td$WetDryState,-1))
td                    <- td[td$WetDryState==0 & td$WetDryState.before==0,]

sst deduction function

sen           <- sst_deduction(datetime = td$dtime, temp = td$IntTemp, temp.range = c(-2,30))
# use sen <- sen[sen$SST.remove==F,] if you want to use the build in removal of unlikely SST values

---------------------- BRAND BAS/biotrack ----------------

load twilight events

trn           <- trn_to_dataframe('trnfile') 

if no trn file available load ligfile and use twilightCalc function from the GeoLight library

load activity file

act           <- read.csv("actfile",sep=",", header = T,row.names = NULL)
act$dtime     <- as.POSIXct(strptime(act[,2], format = "%d/%m/%y %H:%M:%S"),tz='UTC')
act$wetdry    <- act[,4] # wet dry data column must be called 'wetdry'
act           <- act[!is.na(act$dtime),]

load temperature file

td            <- read.csv("temfile",sep=",", header = T,row.names = NULL)
td$dtime      <- as.POSIXct(strptime(td[,2], format = "%d/%m/%y %H:%M:%S"),tz='UTC')

sst deduction function

sen           <- sst_deduction(datetime = td$dtime, temp = td[,4], temp.range = c(-2,30))
# use sen <- sen[sen$SST.remove==F,] if you want to use the build in removal of unlikely SST values

---------------------- BRAND Migrate technology ----------

load twilight events

trn           <- trn_to_dataframe('trnfile') 

if no trn file available load luxfile and use twilightCalc function from the GeoLight library

load activity file

act           <- read.table('degfile',sep="\t",skip=19,header=T) 
act$dtime     <- as.POSIXct(strptime(act[,1], format = "%d.%m.%Y %H:%M:%S"),tz='UTC') # datetime object needs to be POSIXct in UTC and must be called 'dtime'
act$wetdry    <- act$wets0.20 # wet dry data column must be called 'wetdry'
act           <- act[!is.na(act$dtime),]

load temperature file

td            <- read.table('sstfile',sep="\t",skip=19,header=T) 
td$dtime      <- as.POSIXct(strptime(td[,1], format = "%d/%m/%Y %H:%M:%S"),tz='UTC')

sst deduction function

sen           <- sst_deduction(datetime = td$dtime, temp = td[,4], temp.range = c(-2,60))
# use sen <- sen[sen$SST.remove==F,] if you want to use the build in removal of unlikely SST values

---------------------- run prob_algorithm ---------------

# twilight error distribution estimation 
tw            <- twilight_error_estimation()

# see ?prob_algorithm for details
pr <- prob_algorithm(particle.number             = 1000,
                     iteration.number            = 50,
                     trn                         = trn, 
                     sensor                      = sen,  
                     act                         = act, 
                     loess.quartile              = NULL,
                     tagging.location            = c(15.15, 78.17),          # change to colony lon and lat
                     tagging.date                = as.POSIXct("2008-06-01"), # change to tagging date
                     retrieval.date              = as.POSIXct("2015-07-05"), # change to retrieval date
                     wetdry.resolution           = 3,                        # change to sampling resolution of conductivity switch in sec
                     sunrise.sd                  = tw,
                     sunset.sd                   = tw,
                     range.solar                 = c(-7,-1),
                     speed.wet                   = c(1 , 0.4, 2.5),
                     speed.dry                   = c(20, 0.2, 25), 
                     sst.sd                      = 0.5,       
                     max.sst.diff                = 3,          
                     days.around.spring.equinox  = c(10,10),   
                     days.around.fall.equinox    = c(10,10), 
                     ice.conc.cutoff             = 1,                        # ranging from 0 to 1
                     boundary.box                = c(-180,180,-90,90),
                     land.mask                   = T,
                     med.sea                     = F,        
                     black.sea                   = F,        
                     baltic.sea                  = F,      
                     caspian.sea                 = F,    
                     east.west.comp              = T,
                     backward                    = F,
                     NOAA.OI.location            = 'NetCDF file location')

---------------------- plotting --------------------------

plot lat, lon, SST vs time

plot_timeline(pr,degElevation = -3)

plot lon vs lat map

plot_map(pr)
Clone this wiki locally