Section 8.5: Assessing annual trends in ocean chlorophyll levels

Sujit K. Sahu

University of Southampton

2022-03-30

Abstract

This file contains all the code for reproducing the figures and tables for assessing trends in ocean chlorophyll as discussed in Section 8.5.

## Loading required package: Rcpp
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
## Warning in .recacheSubclasses(def@className, def, env): undefined subclass
## "numericVector" of class "Mnumeric"; definition not updated
## 
## ## spTimer version: 3.3.1
## Loading required package: sp
## rgdal: version: 1.5-23, (SVN revision 1121)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: /usr/share/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: /usr/share/proj
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.

figurepath <- "figures/chap8figures/ocean_chl/"
filepath <- "txttables/"
# The  folder below contains the results of MCMC runs 
outputpath <- "mcmc_output_files/ocean_chl/"
mappath <- "mapfiles"
dpath <- "datafiles/"
savename <- paste0(dpath, "C8_chl_datasets.Rdata")
load(savename)

1 Code to run the models.

This takes more than 24 hours.


# "datasets"  "knots"     "lats"      "Longhurst" "lons"      "path"      "savename" 
# Loaded  the above  
## Go to model run 

area <- unique(datasets$area) # 23 unique areas  
a <- table(datasets$area)  ## Number of data rows in each area 
b <- table(knots$area)  ## Number of knots in each area
cbind(a, b) ## See the numbers of data rows and knots 
# options(contrasts=c("contr.treatment", "contr.poly")) ## Month is a factor
options(contrasts=c("contr.sum", "contr.poly")) ## Month is a factor
datasets$month <- as.factor(datasets$month)
datasets$chl <- as.numeric(datasets$chl)


newknots <- F
N <- 2000
burn.in <- 1000
islm <- F # Should it fit linear models 
# Run this code twice once with islm <- T and another time with islm =F
pb <- txtProgressBar(min = 0, max = length(area), style = 3)   # set progress bar
narea <- length(area)

for(k in 1:narea){
# for(k in 22:23)
  #  j <- area[1]     
  #  k <- 8 
  # k <- 1
  
  j <- area[k]
  
  dataj <- datasets[datasets$area==j, ]
  dim(dataj) 
  sites <- unique(dataj$site)
  
  setTxtProgressBar(pb, k)
if (newknots) { 
    ### Set up new knots
    maxlon <- max(dataj$Longitude)
    minlon <- min(dataj$Longitude)
    maxlat <- max(dataj$Latitude)
    minlat <- min(dataj$Latitude)
    
    knotgrid <-spT.grid.coords(Longitude=c(maxlon-0.1, minlon+0.1),
                               Latitude=c(maxlat-0.1, minlat+0.1), 
                               by=round(c(maxlon-minlon, maxlat-minlat)/4)) #4 degree grid spacing
    ###selcting only those knots that are within  the required region.
    #drawing polygon outline region
    longuse[which(is.na(longuse))] <- 0 
    cont <- contourLines(lons,lats,longuse,levels=1)
    temp2 <- 0
    for(i in 1:length(cont)){
      temp <- cont[[i]]
      located <- point.in.polygon(knotgrid[,1],knotgrid[,2],temp$x,temp$y)
      located <- located + temp2
      temp2 <- located
    }
    ind <- located!=0
    knotgrid <- knotgrid[ind,]
} else { 
    knotgrid <- as.matrix(knots[knots$area==j, 1:2])
} # if new knots 
  
  cat("k=", k,  "j=", j, "n=", nrow(dataj), "n locations=", length(sites), "nknots=", nrow(knotgrid), "\n")
  
  #time series setup
  # time.data <- spT.time(t.series=196, segments=1) 

  sst <- dataj$SST
  sst <- scale(sst, center=mean(sst), scale=sd(sst))
  dataj$sst < sst 
  f1 <- chl ~ time+sst + month
  if (islm) { 
     lmodel <- Bsptime(model="lm", data=dataj, formula =f1, N=N)
     savename <- paste0(outputpath, j,"_BGC_model_lm.Rdata")
     save(lmodel, file=savename)
  } else { 
   Model <- Bsptime(package="spTimer", data=dataj, formula =f1, model="GPP", 
                    coordtype="lonlat", coords=2:3, 
                    N=N, burn.in=burn.in, knots.coords=knotgrid, n.report=10, 
                    prior.phi ="FIXED", prior.phi.param=0.006)
   savename <- paste0(outputpath, j,"_spatial_model_nm.Rdata")
   save(Model, file=savename)
  }
 
}
close(pb)

    
round(Model$params[2, ]*100*12, 4)
## Model Run complete 

## Check convergence 
## 

j <- 51
savename <- paste0(outputpath, j,"_spatial_model_nm.Rdata")
load(savename)
names(Model)
plot(Model$fit)


j <- 46
savename <- paste0(outputpath,  j,"_spatial_model_nm.Rdata")
load(savename)
plot(Model$fit)


j <- 33
savename <- paste0(outputpath, j,"_spatial_model_nm.Rdata")
load(savename)
plot(Model$fit)
  


j <- 51
savename <- paste0(outputpath, j,"_spatial_model_nm.Rdata")
load(savename)
plot(Model$fit)


j <- 33
dataj <- datasets[datasets$area==j, ]
dim(dataj)
head(dataj)
summary(dataj)
hist(dataj$chl)
boxplot(dataj$chl)

2 Gather the results from the model runs.


area <- unique(datasets$area) # 23 unique areas  
narea <- length(area)
a <- table(datasets$area)  ## Number of data rows in each area 
b <- table(knots$area)  ## Number of knots in each area
cbind(a, b) ## See the numbers of data rows and knots 
# options(contrasts=c("contr.treatment", "contr.poly")) ## Month is a factor
options(contrasts=c("contr.sum", "contr.poly")) ## Month is a factor
datasets$month <- as.factor(datasets$month)
datasets$chl <- as.numeric(datasets$chl)



pb <- txtProgressBar(min = 0, max = narea, style = 3)   # set progress bar
if (islm) {  
  results <- matrix(NA, nrow=23, ncol=8)
  oparms <- matrix(NA, nrow=23, ncol=6)
} else { 
  results <- matrix(NA, nrow=23, ncol=14)
  oparms <- matrix(NA, nrow=23, ncol=9)
}

for(k in 1:narea) {
  
  j <- area[k]
  cat("k=", k, "j=", j, "\n")
  setTxtProgressBar(pb, k)
  
  dataj <- datasets[datasets$area==j, ]
  dim(dataj)
  
  knotgrid <- as.matrix(knots[knots$area==j, 1:2])
  coords <- unique(dataj[, 2:3])
  dim(coords)
  if (islm)  { 
    savename <- paste0(outputpath, j,"_BGC_model_lm.Rdata")
    load(savename)
    betas <- lmodel$params[2, c(1, 3, 4)] * 100* 12
    print(betas)
    # mchoice <- lmodel$mchoice
    pms <- c( lmodel$params[1, c(1:2)], lmodel$params[15, c(1, 3, 4)])
  } else { 
    savename <- paste0(outputpath, j,"_spatial_model_nm.Rdata")
    load(savename)
    names(Model)
    head(Model$params)
    betas <- Model$params[2, c(1, 3, 4)] * 100* 12 
    pms <- c( Model$params[1, c(1, 3)], Model$params[16, c(1, 3, 4)], Model$params[17, c(1, 3, 4)])
    # mchoice <- Model$mchoice
  }
  signif <- 1
  if ( (betas[1] < 0) & (betas[3]>0) ) signif <- 0
  nt <- nrow(na.omit(dataj))
  if (islm) { 
    a <- as.numeric(c(j, nrow(dataj), nrow(coords), nrow(knotgrid), betas, signif))
  } else { 
    rho <-  as.numeric(Model$params[15, c(1, 3, 4)])
    bsst <- as.numeric(Model$params[3, c(1, 3, 4)])
    a <- as.numeric(c(j, nrow(dataj), nrow(coords), nrow(knotgrid), betas, signif, rho, bsst)) 
  }
  oparms[k, ] <- as.numeric(c(j, pms))
  results[k, ] <- a
  print(round(a, 2))
}

close(pb)



results <- data.frame(results)
# colnames(results) <-  c("area","nT", "n", "knot", "trend","lower","upper","significant", "gof", "penalty", "pmcc")
if (islm) { 
  colnames(results) <-  c("area","nT", "n", "knot", "trend","lower","upper","significant") 
  colnames(oparms) <-  c("area","intercept", "sd", "sige", "lower","upper")
} else { 
  colnames(results) <-  c("area","nT", "n", "knot", "trend","lower","upper","significant", "rho", "rholow", "rhoup", "sst", "lsst", "rsst")
  colnames(oparms) <-  c("area","intercept", "sd", "sige", "lower","upper","sigw", "low", "up")
}
results
oparms

if (islm) { 
  lmresults <- results
  lmoparms <- oparms
  save(lmresults,  lmoparms, file=paste0(outputpath, "C8_chl_lm_results.Rdata"))
  save(lmresults,  lmoparms, file=paste0(dpath, "C8_chl_lm_results.Rdata"))
} else { 
  spresults <- results
  spoparms <- oparms
}
# round(results, 3)
# library(xtable)

u <- results[, c(3, 5, 6, 7, 12, 13, 14, 9, 10, 11)]
# round(u, 2)
# xtable(u, digits = 3)
# xtable(u, digits = 2)
# head(results)

# head(spoparms)
# round(spoparms, 4)
# xtable(spoparms[, -1], digits=4)
save(spresults,  spoparms, file=paste0(outputpath, "C8_chl_sp_results.Rdata"))
save(spresults,  spoparms, file=paste0(dpath, "C8_chl_sp_results.Rdata"))

3 Tables 8.7 and 8.7 and Figure 8.19


long<- readOGR(dsn=path.expand(mappath), layer="Longhurst_world_v4_2010")
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/sks/Dropbox/sks/bookproject/bmstdr_Sahu/chapter_code_files/mapfiles", layer: "Longhurst_world_v4_2010"
## with 54 features
## It has 2 fields
long_df <- tidy(long)
## Regions defined for each Polygons
wmap<- readOGR(dsn=path.expand(mappath), layer="ne_110m_land")
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/sks/Dropbox/sks/bookproject/bmstdr_Sahu/chapter_code_files/mapfiles", layer: "ne_110m_land"
## with 127 features
## It has 2 fields
wmap_df <- tidy(wmap)
## Regions defined for each Polygons
wmap_df <-wmap_df[-(which(wmap_df$group==112.2)),]#remove caspian sea



#trend import and set trends where   0 with in credibility interval to equal 0 
load(paste0(dpath, "C8_chl_sp_results.Rdata"))

table8.6  <- spresults[, c(3, 5, 6, 7, 12, 13, 14, 9, 10, 11)]
dput(table8.6, file=paste0(filepath, "table8.6.txt"))
table8.6
##       n        trend        lower         upper          sst         lsst
## 1   173 -0.106102100 -0.142728882 -7.037372e-02 -0.046986457 -0.050530231
## 2  1586 -0.046088147 -0.065460073 -3.103700e-02 -0.079130642 -0.079650117
## 3  1526 -0.068805369 -0.075776728 -6.215582e-02 -0.020537675 -0.020832965
## 4  1180 -0.119802936 -0.141456777 -1.031129e-01 -0.062884750 -0.063632214
## 5  1263 -0.146522742 -0.168255852 -1.279857e-01 -0.058924672 -0.059614534
## 6   412  0.005101933 -0.099497609  1.077368e-01 -0.025383000 -0.030018242
## 7   406  0.017854971 -0.039989562  6.532805e-02  0.040567496  0.038376049
## 8   307 -0.198794606 -0.320731250 -8.219274e-02 -0.219575994 -0.223843846
## 9   441 -0.052445696 -0.102055295  9.727836e-05 -0.227320521 -0.231175035
## 10  408  0.184953303  0.117369008  2.602423e-01 -0.077769789 -0.083435143
## 11  671 -0.001178155 -0.043011111  3.134055e-02 -0.046990226 -0.048382075
## 12 1311 -0.028870493 -0.035998722 -2.243106e-02  0.008542313  0.008278530
## 13  703  0.379061048  0.256525285  4.804198e-01 -0.242975063 -0.247961708
## 14  556 -0.038120110 -0.048786360 -2.858894e-02 -0.104365436 -0.105210184
## 15  116  0.011091832 -0.032688067  5.610633e-02 -0.174611044 -0.180757150
## 16  439 -0.045750066 -0.102699380  1.082080e-02 -0.071497448 -0.074697309
## 17 1794 -0.009312275 -0.029882176  1.748509e-02 -0.006384200 -0.007052242
## 18 3222  0.006924106 -0.002380876  1.267688e-02 -0.016758293 -0.017078511
## 19  294  0.108547013  0.032801734  1.791103e-01 -0.213098569 -0.220958425
## 20  767  0.038991481  0.011085737  6.618488e-02 -0.071161138 -0.072638154
## 21 1259 -0.067475723 -0.072860827 -6.123082e-02 -0.056094445 -0.056486846
## 22 3785  0.087890101  0.080860212  9.502144e-02 -0.017532131 -0.017899531
## 23 1018 -0.015568382 -0.024296274 -7.079297e-03 -0.010212372 -0.010557788
##            rsst       rho    rholow     rhoup
## 1  -0.043867672 0.6792846 0.6259729 0.7260672
## 2  -0.078587439 0.8619465 0.8540520 0.8698209
## 3  -0.020252972 0.8667368 0.8585748 0.8745530
## 4  -0.062119207 0.7819697 0.7712434 0.7927855
## 5  -0.058146090 0.8246454 0.8142565 0.8346552
## 6  -0.020343142 0.7805984 0.7613509 0.7994939
## 7   0.042573928 0.8806588 0.8653919 0.8961570
## 8  -0.214871903 0.7535259 0.7289676 0.7775122
## 9  -0.223802746 0.8085399 0.7889748 0.8263388
## 10 -0.072230167 0.4988337 0.4639981 0.5308764
## 11 -0.045535928 0.8732969 0.8618487 0.8853612
## 12  0.008776193 0.8506085 0.8417469 0.8591717
## 13 -0.238146490 0.8847860 0.8738431 0.8955735
## 14 -0.103459906 0.6987095 0.6776266 0.7198695
## 15 -0.168623835 0.5425496 0.4782246 0.6045384
## 16 -0.068024932 0.5113895 0.4839967 0.5393902
## 17 -0.005607702 0.8220268 0.8150430 0.8290603
## 18 -0.016439477 0.9305193 0.9268816 0.9344224
## 19 -0.204723789 0.6879025 0.6557002 0.7178821
## 20 -0.069442650 0.6587983 0.6406926 0.6768000
## 21 -0.055679954 0.6666128 0.6544032 0.6798979
## 22 -0.017208117 0.8967740 0.8733662 0.9200848
## 23 -0.009864346 0.9078956 0.8992100 0.9158933

table8.7 <- spoparms[, -1]
round(table8.7, 4)
##       intercept     sd   sige  lower  upper   sigw    low     up
##  [1,]    0.2903 0.2869 0.0095 0.0093 0.0096 0.0164 0.0148 0.0181
##  [2,]    0.1664 0.1646 0.0020 0.0020 0.0020 0.0053 0.0052 0.0055
##  [3,]    0.1153 0.1145 0.0005 0.0005 0.0005 0.0020 0.0019 0.0020
##  [4,]    0.1846 0.1828 0.0022 0.0022 0.0023 0.0060 0.0059 0.0062
##  [5,]    0.1805 0.1787 0.0035 0.0035 0.0035 0.0076 0.0074 0.0078
##  [6,]    0.2813 0.2703 0.0363 0.0360 0.0366 0.0634 0.0601 0.0669
##  [7,]    0.2817 0.2774 0.0161 0.0159 0.0162 0.0201 0.0190 0.0214
##  [8,]    0.2900 0.2788 0.0643 0.0636 0.0651 0.0776 0.0729 0.0826
##  [9,]    0.2594 0.2551 0.0128 0.0127 0.0130 0.0219 0.0209 0.0230
## [10,]    0.3666 0.3594 0.0352 0.0349 0.0355 0.0896 0.0847 0.0948
## [11,]    0.1584 0.1549 0.0122 0.0121 0.0123 0.0200 0.0192 0.0209
## [12,]    0.0876 0.0869 0.0003 0.0003 0.0003 0.0014 0.0014 0.0015
## [13,]    0.3117 0.3039 0.0669 0.0664 0.0674 0.0732 0.0702 0.0764
## [14,]    0.1041 0.1032 0.0008 0.0008 0.0008 0.0044 0.0043 0.0046
## [15,]    0.3015 0.2965 0.0111 0.0109 0.0113 0.0368 0.0334 0.0405
## [16,]    0.4320 0.4248 0.0225 0.0223 0.0228 0.1078 0.1025 0.1130
## [17,]    0.1987 0.1969 0.0046 0.0046 0.0047 0.0178 0.0174 0.0182
## [18,]    0.0910 0.0903 0.0007 0.0007 0.0007 0.0015 0.0014 0.0015
## [19,]    0.4002 0.3926 0.0749 0.0741 0.0759 0.3637 0.3405 0.3875
## [20,]    0.2352 0.2323 0.0055 0.0055 0.0055 0.0170 0.0165 0.0177
## [21,]    0.0774 0.0767 0.0003 0.0003 0.0003 0.0021 0.0020 0.0021
## [22,]    0.2147 0.2140 0.0150 0.0149 0.0150 0.6205 0.5002 0.7377
## [23,]    0.0837 0.0828 0.0005 0.0005 0.0006 0.0014 0.0014 0.0015
dput(table8.7, file=paste0(filepath, "table8.7.txt"))

islm <- F  ## Do for the spatial model 
if (islm)  { trend <- lmresults$trend
} else trend <-  spresults$trend
long_df$area <- as.numeric(long_df$id)
long_df$trend <- NA
#sorting differing region orders(stephs longhurst --> marineplan longhurst order-->our provnice order)
provsorder <- c(NA,NA,13,11,15,4,10,1,8,NA,NA,NA,NA,NA,NA,NA,14,NA,NA,NA,2,3,NA,NA,NA,NA,NA,NA,NA,18,19,12,16,17,23,20,6,5,7,9,NA,NA,NA,NA,NA,NA,NA,NA,NA,21,22,NA,NA,NA)
for(i in 1:54){
  long_df$trend[which(long_df$id==i)] <- trend[provsorder[i]]
}


#centrepoints for region labelling
long_df$lonc <- NA
long_df$latc <- NA
centroids_df <- as.data.frame(coordinates(long))

centroids_df[5,] <- c(-55.3,40.7) #small number of regions centrepoints not optimal
centroids_df[33,] <- c(149.3,39) #small number of regions centrepoints not optimal
centroids_df[36,] <- c(160.6,-36.3)
centroids_df[51,] <- c(75,-36.2)
for(i in c(4,5,6,7,8,9,10,18,22,23,31,32,33,34,35,36,37,38,39,40,41,51,52)){
  long_df$lonc[which(long_df$area==i)] <- centroids_df$V1[i]
  long_df$latc[which(long_df$area==i)] <- centroids_df$V2[i]
}
#correct labels
long_df$truid <- NA
provsorder <- c(NA,NA,NA,13,11,15,4,10,1,8,NA,NA,NA,NA,NA,NA,NA,14,NA,NA,NA,2,3,NA,NA,NA,NA,NA,NA,NA,18,19,12,16,17,23,20,6,5,7,9,NA,NA,NA,NA,NA,NA,NA,NA,NA,21,22,NA,NA)#slightly different due to null region

for(i in 1:54){
  long_df$truid[which(long_df$id==i)] <- provsorder[i]
}
# range(trend)



#long_df$truid <- long_df$area
###create plot
com <- rev(c("firebrick4","firebrick2","white","dodgerblue2","dodgerblue4"))#colour palette
chlmap <-  ggplot( data=long_df, aes(x=long, y=lat, group = group,fill=trend)) +
  scale_fill_gradientn(colours=com,na.value="black",limits=c(-0.2, 0.4))+
  geom_polygon(colour='black',size=0.25)+ #longhurst outlines with half width size
  geom_text(data=na.omit(long_df),aes(label = truid, x = lonc, y = latc,family="Times"),size=4) +
  geom_polygon(data=wmap_df, aes(x=long, y=lat, group = group),colour="black", fill="gainsboro",size=0.25 )+
  coord_equal() + guides(fill=guide_colorbar(title=expression(paste("Trend  ", "(%yr"^" -1",")"))))+
  ylab("Latitude")+xlab("Longitude")+
  coord_cartesian(ylim = c(-75, 75),expand=F)+
  scale_y_continuous(breaks=seq(-60,60,30))  +
  scale_x_continuous(breaks=seq(-180,180,60))  +
  annotate("text",x=-133,y=-38.1,label="21",size=4,colour='black',family="Times")+ #extra label for region 21 to show that it wraps
  theme_bw()+theme(text=element_text(family="Times"))
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
chlmap

ggsave(paste0(figurepath, "space_trend_longhurst.png"),  width=8.27, height=3.44, dpi=300)
# 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.41385