Abstract

This file contains all the code for reproducing the figures and tables for modeling temperature data from roaming ocean Argo floats as discussed in Section 8.6.


## Note the start time
start.time<-proc.time()[3]
# figurepath <- "YourPath"   # Folder to save or read the figures from
# filepath <- "YourFilePath" # Folder to save or read the tables from 
figurepath <- "figures/chap8figures/argo_floats_2003/"
filepath <- "txttables/"
longrun <- FALSE  # Should we run the model to check the results. 
# If FALSE we use the results, e.g. tables, figures and model fits from previous runs

library(bmstdr)
## Warning in .recacheSubclasses(def@className, def, env): undefined subclass
## "numericVector" of class "Mnumeric"; definition not updated
library(ggplot2)
library(akima)
library(tidyr)
library(RColorBrewer)
library(ggpubr)
library(GGally)

1 Code to reproduce Figure 8.20


deep <- argo_floats_atlantic_2003[argo_floats_atlantic_2003$depth==3, ]
deep$month <- factor(deep$month)

 ## temperature
 ptmp <- ggplot(deep, aes(x=lon, y=temp, color=month), size=1) +
   geom_point(shape=16, size=1) +
   stat_smooth(method = "lm", col = "red") + 
   labs( x = "Longitude",  y= "Temperature") +
   theme(legend.position = "right", legend.key.size = unit(0.2, "cm"))
 
 ptmp
## `geom_smooth()` using formula 'y ~ x'

 
 ggsave(paste0(figurepath, "lon_v_temp_scatter.png"),  width=4, height=3, dpi=300)
## `geom_smooth()` using formula 'y ~ x'

 ## max temperature
 plat <- ggplot(deep, aes(x=lat, y=temp, color=month), size=1) +
   geom_point(shape=16, size=1) +
   stat_smooth(method = "lm", col = "red") + 
   labs( x = "Latitude",  y= "Temperature") +
   theme(legend.position = "right", legend.key.size = unit(0.2, "cm"))
 
 plat
## `geom_smooth()` using formula 'y ~ x'

 ggsave(paste0(figurepath, "lat_v_temp_scatter.png"),  width=4, height=3, dpi=300)
## `geom_smooth()` using formula 'y ~ x'
 
 ggarrange(plat, ptmp, nrow=1, common.legend = TRUE, legend = "bottom")
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

 ggsave(paste0(figurepath, "lat_and_lon_v_temp_scatter.png"),  width=4, height=3, dpi=300)

2 Reproduce Figure 8.21


if (longrun) {
deep$x2inter <- deep$xinter*deep$xinter
deep$month <- factor(deep$month)

deep$lat2 <- (deep$lat)^2
deep$sin1 <- round(sin(deep$time*2*pi/365), 4)
deep$cos1 <- round(cos(deep$time*2*pi/365), 4)
deep$sin2 <- round(sin(deep$time*4*pi/365), 4)
deep$cos2 <- round(cos(deep$time*4*pi/365), 4)

head(deep)
## scaling variables 
deep[, c( "xlat2", "xsin1", "xcos1", "xsin2", "xcos2")] <- 
  scale(deep[,c("lat2", "sin1", "cos1", "sin2", "cos2")])

f2 <- temp ~ xlon + xlat + xlat2+ xinter + x2inter 
M2atl <- Bmoving_sptime(formula=f2, data = deep, coordtype="lonlat", coords = 1:2,
                N=1100, burn.in=100, validrows =NULL, mchoice =T)
summary(M2atl)
# plot(M2atl)
save(M2atl, deep, file="argo/argo.RData")
names(M2atl)

table8.8 <- M2atl$params
dput(table8.8, file=paste0(filepath, "table8.8.txt"))
} else {
  load(file="argo/argo.RData")
}
listofdraws <- rstan::extract(M2atl$fit)
names(listofdraws)
## [1] "beta"     "phi"      "sigma_sq" "tau_sq"   "xbmodel"  "bigS"     "lp__"
dim(listofdraws$xbmodel)
## [1] 1000 2378
dim(listofdraws$bigS)
## [1]  1000 19666
dat <- M2atl$datatostan
names(dat)
##  [1] "n"               "tn"              "m2"              "p"              
##  [5] "missing"         "ntmiss"          "ntobs"           "data_miss_idx"  
##  [9] "data_obs_idx"    "time"            "nots"            "ots"            
## [13] "nts"             "start_row"       "fin_row"         "n_misst"        
## [17] "misst"           "Cdist"           "dmat"            "yobs"           
## [21] "X"               "sigma2_prior"    "tau2_prior"      "phidist"        
## [25] "prior_phi_param"

## Extract the diagonal elements of all St 
v <- numeric()
x <- numeric()
m <- length(dat$nts)
for (i in 1:m) {
  k <- dat$nts[i]
  a1 <- 1:k^2
  a2 <- seq(from=1, to =k^2, length=k)
  b1 <- rep(0, k^2)
  b1[which(a1==a2)] <- 1
  cbind(a1, b1)
  v <- c(v, a1) 
  x <- c(x, b1) ## indicates if the corresponding index is a diagonal 
}


varsamps <- listofdraws$bigS[, which(x>0)]/365
dim(varsamps)
## [1] 1000 2378
        xbeta <- listofdraws$xbmodel
dim(xbeta)
## [1] 1000 2378
sdsamps <- sqrt(varsamps)

ntot <- nrow(xbeta) * ncol(xbeta)

zsamp <- matrix(rnorm(ntot), nrow=nrow(xbeta), ncol=ncol(xbeta))
dim(zsamp)
## [1] 1000 2378
zsamp <- xbeta + zsamp * sdsamps
ansummary <- get_parameter_estimates(zsamp)
head(ansummary)
##       mean         sd      low       up
## 1 3.215707 0.09122472 3.050348 3.405181
## 2 4.768590 0.06088208 4.650099 4.881378
## 3 6.939425 0.07800954 6.780590 7.090811
## 4 5.528516 0.06879682 5.394417 5.668185
## 5 6.921571 0.05874630 6.813468 7.035608
## 6 4.086425 0.07112330 3.944540 4.228669
summary(ansummary$mean)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.026   5.117   6.496   6.273   7.705  11.016
summary(ansummary$sd)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.05326 0.06040 0.06593 0.07008 0.07669 0.12929
summary(ansummary$low)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.847   4.993   6.379   6.137   7.572  10.816
summary(ansummary$up)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.208   5.236   6.615   6.411   7.846  11.223


pdata <- cbind(deep, ansummary)


atlmap <- map_data("world", xlim=c(-70, 10), ylim=c(15, 65))
head(atlmap)
##        long      lat group order   region subregion
## 1 -63.00122 18.22178     1     1 Anguilla      <NA>
## 2 -63.16001 18.17139     1     2 Anguilla      <NA>
## 3 -63.15332 18.20029     1     3 Anguilla      <NA>
## 4 -63.02603 18.26973     1     4 Anguilla      <NA>
## 5 -62.97959 18.26480     1     5 Anguilla      <NA>
## 6 -63.00122 18.22178     1     6 Anguilla      <NA>
atlmap <- atlmap[atlmap$long < 5, ]
atlmap <- atlmap[atlmap$long > -70, ]
atlmap <- atlmap[atlmap$lat < 65, ]
atlmap <- atlmap[atlmap$lat > 10, ]

xo <- seq(from=-70, to = 5,  length=200)
yo <- seq(from= 10, to = 65, length=200)

## Temperature 
surf <- interp(pdata$lon, pdata$lat, pdata$mean,  xo=xo, yo=yo)
names(surf)
## [1] "x" "y" "z"
surf <- list(x=surf$x, y=surf$y, z=surf$z)
interp1 <- data.frame(long = surf$x, surf$z )
names(interp1)[1:length(surf$y)+1] <- surf$y
interp1 <- gather(interp1,key = lat,value =Temp,-long,convert = TRUE)

head(interp1)
##        long lat Temp
## 1 -70.00000  10   NA
## 2 -69.62312  10   NA
## 3 -69.24623  10   NA
## 4 -68.86935  10   NA
## 5 -68.49246  10   NA
## 6 -68.11558  10   NA
summary(interp1$Temp)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   2.028   5.223   6.850   6.569   7.884  11.012   26248
pcolors <- topo.colors(5)
P <- ggplot() +  
  geom_raster(data=interp1, aes(x = long, y = lat,fill = Temp)) + 
  # scale_fill_gradientn(colours=c("#0000FFFF","#FFFFFFFF","#FF0000FF")) + 
  scale_fill_gradientn(colours=pcolors) +  
  geom_contour(data=interp1, aes(x = long, y = lat,z = Temp)) + 
  geom_polygon(data=atlmap, aes(x=long, y=lat, group=group),
               color="black", size = 0.6, fill=NA) + 
    geom_point(data =deep, aes(x=lon, y=lat, colour=month), size=1) +
  labs( title= "Annual temperature in deep ocean in 2003", x="Longitude", y = "Latitude") +
  ggsn::scalebar(data =atlmap, dist = 1000, location = "bottomleft", transform=T, dist_unit = "km",
                 st.dist = .05, st.size =5, height = .05, st.bottom=T, model="WGS84") +
  ggsn::north(data=atlmap, location="topright", symbol=12) 

P 
## Warning: Removed 26248 rows containing non-finite values (stat_contour).


ggsave(paste0(figurepath, "temp_deep.png"),  width=8.27, height=5, dpi=300)
## Warning: Removed 26248 rows containing non-finite values (stat_contour).

## sd of temperature 

surf <- interp(pdata$lon, pdata$lat, pdata$sd,  xo=xo, yo=yo)
names(surf)
## [1] "x" "y" "z"
surf <- list(x=surf$x, y=surf$y, z=surf$z)
interp1 <- data.frame(long = surf$x, surf$z )
names(interp1)[1:length(surf$y)+1] <- surf$y
interp1 <- gather(interp1,key = lat,value =sd,-long,convert = TRUE)

head(interp1)
##        long lat sd
## 1 -70.00000  10 NA
## 2 -69.62312  10 NA
## 3 -69.24623  10 NA
## 4 -68.86935  10 NA
## 5 -68.49246  10 NA
## 6 -68.11558  10 NA
summary(interp1$sd)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.054   0.062   0.069   0.073   0.081   0.128   26248

psd <- ggplot() +  
  geom_raster(data=interp1, aes(x = long, y = lat,fill = sd)) + 
  # scale_fill_gradientn(colours=c("#0000FFFF","#FFFFFFFF","#FF0000FF")) + 
  scale_fill_gradientn(colours=pcolors) + 
  geom_contour(data=interp1, aes(x = long, y = lat,z = sd)) + 
  geom_polygon(data=atlmap, aes(x=long, y=lat, group=group),
               color="black", size = 0.6, fill=NA) + 
  geom_point(data =deep, aes(x=lon, y=lat, colour=month), size=1) +
  labs( title= "sd of annual temperature in deep ocean in 2003", x="Longitude", y = "Latitude") +
  ggsn::scalebar(data = atlmap, dist = 1000, location = "bottomleft", transform=T, dist_unit = "km",
                 st.dist = .05, st.size =5, height = .05, st.bottom=T, model="WGS84") +
  ggsn::north(data=atlmap, location="topright", symbol=12) 

psd 
## Warning: Removed 26248 rows containing non-finite values (stat_contour).


ggsave(paste0(figurepath, "sd_temp_deep.png"),  width=8.27, height=5, dpi=300)
## Warning: Removed 26248 rows containing non-finite values (stat_contour).

# All done 
end.time <- proc.time()[3]
comp.time<- (end.time-start.time)/60
# comp.time<-fancy.time(comp.time)
print(comp.time)
##   elapsed 
## 0.2494167

Section 8.6: Modeling temperature data from roaming ocean Argo floats

Sujit K. Sahu

2022-03-30

1 Code to reproduce Figure 8.20

2 Reproduce Figure 8.21