## 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