Skip to content
Tim Lam edited this page Apr 4, 2018 · 4 revisions

Obtaining confidence intervals and regions

Assuming you have obtained a converged "fit", here's how you can get the CIs:

    lat <- fit$most.prob.track[, 2]
    std <- sqrt(fit$var.most.prob.track[, 4])
    low <- lat - 2 * std
    hig <- lat + 2 * std

    lon <- fit$most.prob.track[, 1]
    std <- sqrt(fit$var.most.prob.track[, 1])
    low <- lon - 2 * std
    hig <- lon + 2 * std

To export the confidence regions into a shapefile, use the CI2shp in Analyzepsat:

if (!any('devtools' == installed.packages()[,"Package"])){
	install.packages('devtools', repos="http://cran.rstudio.com/")}
require(devtools)

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

library(analyzepsat)
### Apply temporary fix
source_url("https://raw.githubusercontent.com/positioning/kalmanfilter/master/support/analyzepsat-hotfix.r")

### Get an example run of ukfsst
### ------------------------------------------------------------------------
library(ukfsst)
example(blue.shark)
track <- blue.shark
print(fit)

### Obtain data from the model fit
### ------------------------------------------------------------------------
fmat <- data.frame(track[,c(3,2,1)],
                   fit$var.most.prob.track,fit$most.prob.track,0,0)
names(fmat) <- c("Year", "Month", "Day", "V11", "V12", "V21", 
                       "V22", "Lon_E", "Lat_N", "max_depth", "SST")

### Generate shapefile
### ------------------------------------------------------------------------
### Shapefile will be saved in working directory
### Use function getwd() to see where
CI2shp(fmat, fname = "testshp", level = 0.95, npoints = 100, 
       proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs "))