An illustrated script to using analyzepsat for bathymetric correction and utilization distribution (UD) generation

Here is some R code to get started -

### Load libraries
### ------------------------------------------------------------------------
### Install library if needed
# devtools::install_github('galuardi/analyzepsat', dep = T, ref = 'v4.0')
library(analyzepsat)
library(kftrack)

### Temporary fix
### ------------------------------------------------------------------------
 if (!any('devtools' == installed.packages()[,"Package"])){
  install.packages('devtools', repos="http://cran.rstudio.com/")}
require(devtools)
source_url("https://raw.githubusercontent.com/positioning/kalmanfilter/master/support/analyzepsat-hotfix.r")

### Load data
### ------------------------------------------------------------------------
data(big.241)
### Add dummy columns for sea surface temperature (10 C) and fish diving depth (-10 m)
big.241$sst = 10
big.241$maxz=-10

### Get bathymetry from NOAA server. Longitude goes from -180 to 180.
### -----------------------------------------------------------------
bath = get.bath.data(-180,-150,-10,50, res=1)

### Run the Kalman filter routine
### Using only datestamp, longitude and latitude from manufacturer outputs
### ------------------------------------------------------------------------
fit1 = kftrack(big.241[,1:5])

### Prepare for bathymetric correction
### ----------------------------------
ftrack1 = prepb(fit1, big.241)

### Make sure longitude is -180 to 180. We here use the summy column for maxz
ftrack1$Lon_E = ftrack1$Lon_E-360

### Do the bathymetric correction
btrack1 = make.btrack(ftrack1, bath)

### Plotting
### ----------------------------------
data(gmt3, package="trackit")

### Plot 1 - Show entire track, corrected with bathymetry
par(mar=c(4,4,6,4)) 
plot(gmt3, xlim = c(-160, -145)+360, ylim = c(10,25), typ='l', axes=F)
rect(par("usr")[1],par("usr")[3],par("usr")[2],par("usr")[4],col = "lightblue")
plot.btrack(btrack1, offset=360, add=T, ci=T)
polygon(gmt3, col = 'black')
degAxis(1); degAxis(2)
.add.month.scale()

### Zoom into the area with points on land
par(mfrow=c(1,2), mar=c(4,4,6,4))
### Plot 2 - Original solution from Kalman filter
plot(gmt3, xlim = c(-159, -153)+360, ylim = c(20,22), typ='l', axes=F)
rect(par("usr")[1],par("usr")[3],par("usr")[2],par("usr")[4],col = "lightblue")
plot.btrack(ftrack1, offset=360, add=T)
polygon(gmt3, col = 'black')
plot.btrack(ftrack1, offset=360, add=T)
degAxis(1); degAxis(2)
.add.month.scale()
### Plot 3 - Bathymetric correction applied
plot(gmt3, xlim = c(-159, -153)+360, ylim = c(20,22), typ='l', axes=F)
rect(par("usr")[1],par("usr")[3],par("usr")[2],par("usr")[4],col = "lightblue")
plot.btrack(btrack1, offset=360, add=T, ci=T)
polygon(gmt3, col = 'black')
plot.btrack(btrack1, offset=360, add=T)
degAxis(1); degAxis(2)

### Get the utilization distribution
### ----------------------------------
btrack1[,8] = btrack1[,8]+360  # add 360 back to longitude so the gmt polygons match
kd1 = track2KD(btrack1,  xsize= .1,  ysize= .1,  range.x=c(-160, -140)+360, range.y= c(0,30))
ud1 = kern2UD(kd1)

### Plotting
### ----------------------------------
 if (!any('adehabitat' == installed.packages()[,"Package"])){
  install.packages('adehabitat', repos="http://cran.rstudio.com/")}
require(adehabitat)
data(gmt3, package="trackit")

x11()
ilabs = c(.1, .3, .5,.75, .95)*100
axis.args=list( at=ilabs, labels=paste(round(ilabs), '%') )
### Plot UD
image(ud1, zlim = c(0,.99),  xlim = c(-160, -140)+360, ylim = c(10,25),
 col = gray.colors(100,start = 0, end = 1), xlab='', ylab='',add=F, axes=F)
rect(par("usr")[1],par("usr")[3],par("usr")[2],par("usr")[4],col = "lightblue")
image(ud1, zlim = c(0,.99),  xlim = c(-160, -140)+360, ylim = c(10,25),
 col = gray.colors(100,start = 0, end = 1), xlab='', ylab='',add=T, axes=F)
	degAxis(1); 
	degAxis(2);
	polygon(gmt3, col = 'green')
	box()
### Add UD legend	
image.plot(matrix(c(0,100)), col =  gray.colors(100,start = 0, end = 1),
 smallplot = c(.65,.70,.25,.85) , legend.only=T, axis.args=axis.args)

# add the track
data(ATL)
data(myramps)
plot.btrack(btrack1, offset=0, add=T, ci=F)

# now use the fixed kernel density
fkern = kernelUD(btrack1[,8:9])
ukern = getvolumeUD(fkern)

Side-by-side comprison shows that fixed kernel methods overestimate the utilization

Summary

The track2KD and kern2UD functions in analyzepsat utilize the error structure estimated through the KF methods on these pages to directly estimate the utilization (Galuardi et. al, in prep).