Commit 6f9b5bb2 authored by md0276b's avatar md0276b
Browse files

SST&Chla average per PELGAS survey

parent 6ceea289
......@@ -14,6 +14,7 @@ library(rgdal)
library(raster)
library(rasterVis)
library(maptools)
library(sp)
data(gradient_colors)
paletteGC=rasterTheme(region=gradient.colors)
......@@ -33,8 +34,8 @@ pref=echosonde
# Rasters export path
path.export.grid.sat=paste(pref,'Satellite/Atlantic/',sep='')
# 1. Define grid and polygon selection
# to resample sat images on gridmap grid (0.25°)----------
# 1. Define grid and polygon selection ------
# to resample sat images on gridmap grid (0.25°)
#*********************************
#patnam='c:/ptg/work/Atlantic/'
# westernmost_longitude: -12.000
......@@ -90,7 +91,7 @@ legend('topleft',c('Discarded','Selected'),pch=c(1,16))
coastlm=map(database="worldHires",xlim=c(x1,x2),ylim=c(y1,y2),col=gray(.8),fill=T,bg="white",add=T)
coastlsp <- map2SpatialPolygons(coastlm, IDs=coastlm$names,proj4string=CRS("+proj=longlat"))
### 2. Create Chl-a raster stack --------
#2. Create Chl-a raster bricks --------
# **********************************
# 2.1. Define path
patnam=paste(pref,"Satellite/Atlantic/Chla/Rdata/",sep='')
......@@ -214,14 +215,14 @@ cat("Raster stack saved as:", filei2, "\n")
# #r=stack(r,mr)
# r=stack(r,mr,stdr)
# 3. SST rasters (3-days composites) ----
# 3. Create SST rasters bricks ----
#******************************************
# 2.1. Define paths
# 3.1. Define paths
patnam2=paste(pref,"Satellite/Atlantic/SST/Rdata/",sep='')
ficnam2=list.files(path=patnam2) # lire les noms de fichier
ficnam2=ficnam2[substr(ficnam2,17,20)=="data"] #SST
# Grid per species and time-steps (1-day)
# 3.1. Bind raster layers per 1-day time-step --------
#******************************************
if (exists('sst.rasterStack')){
lidik2=names(sst.rasterStack)
......@@ -312,13 +313,13 @@ filei2 = paste(path.export.grid.sat,'SST/raster/',paste('rasterStack',idi,sep='_
# read in a raster stack
sst.rasterStack=stack(filei2)
# 3. Aggregate rasterStacks -------------
# 4. Average rasters over time -------------
# ************************************************
library(raster)
library(chron)
# 3.1. Chla ----------
# 4.1. Chla ----------
#**********************
# 3.1.1. Eventually, import SSChla raster stack ------
# 4.1.1. Eventually, import SSChla raster stack ------
#echosonde='/run/user/1000/gvfs/smb-share:domain=ifr,server=echosonde,share=data,user=mdoray/'
#donnees2='/run/user/1000/gvfs/smb-share:domain=ifr,server=nantes,share=donnees2,user=mdoray/'
#path.export.grid.sat=paste(pref,'Satellite/Data/Atlantic/',sep='')
......@@ -340,7 +341,7 @@ plot(chla.rb)
# Set negative values to NA
chla.rb[chla.rb<=0]=NA
# 3.1.2. Aggregate Chla raster stacks ------------
# Average Chla raster bricks over time
chla.rs.times=data.frame(rb.names=names(chla.rb))
chla.rs.times$year=as.numeric(substr(chla.rs.times$rb.names,2,5))
chla.rs.times$jday=as.numeric(substr(chla.rs.times$rb.names,7,9))
......@@ -382,16 +383,15 @@ chla.rs.aDates30d$myid=paste(chla.rs.aDates30d$Group.2,
#chla.rs
#chla.rs.6d
# weekly averages
chla.rb.6d=aggregate(chla.rb,fact=c(1,1,2))
# much faster
# 4.1.2. Weekly averages ----------
chla.rb.6d=aggregate(chla.rb,fact=c(1,1,2)) # much faster
chla.rb
chla.rb.6d
dim(chla.rs.aDates6d)
# correct layer names
names(chla.rb.6d)=substr(format(chla.rs.aDates6d[,2]),1,10)
# monthly averages:
# 4.1.3. Monthly averages ----------
# with stackApply:
quantile90=function(x,...){quantile(x,probs=.9,na.rm=TRUE)}
lmyid=unique(chla.rs.times$myid)
......@@ -413,7 +413,10 @@ plot(mchla.rb.yearmonth)
plot(sdchla.rb.yearmonth)
plot(q90chla.rb.yearmonth)
# 3.1.2.1. Mean and SD Chla during PELGAS surveys ---------
# or, with 30 days bins
chla.rb.monthyear=aggregate(chla.rb,fact=c(1,1,30))
# 4.1.4. Average during PELGAS surveys ---------
# *********************************************
# Import PELGAS dates
# database login to check reference tables
......@@ -447,33 +450,94 @@ PELGAS.voyages.timingsf=sort(strptime(union(PELGAS.voyages.timings$startdatef,
format="%Y-%m-%d",tz="GMT"))
length(PELGAS.voyages.timingsf);length(PELGAS.voyages.timingsf.labels)
# Define PELGAS intervals
chla.rs.times$inPELGAS=cut(chla.rs.times$time,breaks=PELGAS.voyages.timingsf,
labels=FALSE)
mchla.rb.PELGAS=stackApply(chla.rb,
indices=as.numeric(factor(chla.rs.times$myid)),
fun=mean,na.rm=TRUE)
names(mchla.rb.yearmonth)=lmyid
# 3.1.3. Save mean Chla raster stacks to disk ---------
unique(chla.rs.times$inPELGAS)
# Average chla bricks according to PELGAS intervals
mchla.rb.PELGAScut=stackApply(chla.rb,
indices=as.numeric(factor(chla.rs.times$inPELGAS)),
fun=mean,na.rm=TRUE)
sdchla.rb.PELGAScut=stackApply(chla.rb,
indices=as.numeric(factor(chla.rs.times$inPELGAS)),
fun=sd,na.rm=TRUE)
q90chla.rb.PELGAScut=stackApply(chla.rb,
indices=as.numeric(factor(chla.rs.times$inPELGAS)),
fun=quantile90,na.rm=TRUE)
# Select mean chla layers corresponding to PELGAS intervals
mchla.rb.PELGAS=subset(mchla.rb.PELGAScut,names(mchla.rb.PELGAScut)[
substr(names(mchla.rb.PELGAScut),7,8)%in%
seq(1,max(chla.rs.times$inPELGAS,na.rm = TRUE),2)])
names(mchla.rb.PELGAS)=paste('PELGAS',seq(2000,2019,1),sep='')
plot(mchla.rb.PELGAS)
sdchla.rb.PELGAS=subset(sdchla.rb.PELGAScut,names(sdchla.rb.PELGAScut)[
substr(names(sdchla.rb.PELGAScut),7,8)%in%
seq(1,max(chla.rs.times$inPELGAS,na.rm = TRUE),2)])
names(sdchla.rb.PELGAS)=paste('PELGAS',seq(2000,2019,1),sep='')
plot(sdchla.rb.PELGAS)
q90chla.rb.PELGAS=subset(q90chla.rb.PELGAScut,names(q90chla.rb.PELGAScut)[
substr(names(q90chla.rb.PELGAScut),7,8)%in%
seq(1,max(chla.rs.times$inPELGAS,na.rm = TRUE),2)])
names(q90chla.rb.PELGAS)=paste('PELGAS',seq(2000,2019,1),sep='')
plot(q90chla.rb.PELGAS)
# 4.1.4.1. Select cells in PELGAS polygon -----
data("PELGASpolygon")
PELGASpolygon.sp=Polygon(PELGASpolygon)
PELGASpolygon.sps = Polygons(list(PELGASpolygon.sp), "s1")
PELGASpolygon.sppol = SpatialPolygons(list(PELGASpolygon.sps))
plot(PELGASpolygon.sppol)
mchla.rb.PELGAS0 <- crop(mchla.rb.PELGAS, PELGASpolygon.sppol)
mchla.rb.PELGASs <- mask(mchla.rb.PELGAS0,
PELGASpolygon.sppol)
plot(mchla.rb.PELGASs)
sdchla.rb.PELGAS0 <- crop(sdchla.rb.PELGAS, PELGASpolygon.sppol)
sdchla.rb.PELGASs <- mask(sdchla.rb.PELGAS0,
PELGASpolygon.sppol)
q90chla.rb.PELGAS0 <- crop(q90chla.rb.PELGAS, PELGASpolygon.sppol)
q90chla.rb.PELGASs <- mask(q90chla.rb.PELGAS0,
PELGASpolygon.sppol)
# 4.1.5. Save mean Chla raster stacks to disk ---------
#************************************************
year=1998:2021
idi=paste('mSSChlaYearMonth',paste(min(year),max(year),sep='-'),sep = "_")
filei3 = paste(path.export.grid.sat,'Chla/raster/',paste('rasterStack',idi,sep='_'),'.grd',sep='')
idi=paste('mSSChlaYearMonth',paste(min(years.chla),max(years.chla),sep='-'),
sep = "_")
filei3 = paste(path.export.grid.sat,'Chla/raster/',paste('rasterStack',idi,
sep='_'),'.grd',sep='')
writeRaster(mchla.rb.yearmonth,filei3,overwrite=TRUE)
idi=paste('sdSSChlaYearMonth',paste(min(year),max(year),sep='-'),sep = "_")
filei4 = paste(path.export.grid.sat,'Chla/raster/',paste('rasterStack',idi,sep='_'),'.grd',sep='')
idi=paste('sdSSChlaYearMonth',paste(min(years.chla),max(years.chla),sep='-'),
sep = "_")
filei4 = paste(path.export.grid.sat,'Chla/raster/',paste('rasterStack',idi,
sep='_'),'.grd',sep='')
writeRaster(sdchla.rb.yearmonth,filei4,overwrite=TRUE)
idi=paste('q90SSChlaYearMonth',paste(min(year),max(year),sep='-'),sep = "_")
filei5 = paste(path.export.grid.sat,'Chla/raster/',paste('rasterStack',idi,sep='_'),'.grd',sep='')
idi=paste('q90SSChlaYearMonth',paste(min(years.chla),max(years.chla),sep='-'),
sep = "_")
filei5 = paste(path.export.grid.sat,'Chla/raster/',paste('rasterStack',idi,
sep='_'),'.grd',sep='')
writeRaster(q90chla.rb.yearmonth,filei5,overwrite=TRUE)
# or, with 30 days bins
chla.rb.monthyear=aggregate(chla.rb,fact=c(1,1,30))
# 3.2. Create SST raster stack----------
year1PELGAS=2000;yearNPELGAS=2019
idi=paste('mSSChlaPELGAS',paste(year1PELGAS,yearNPELGAS,sep='-'),
sep = "_")
filei6 = paste(path.export.grid.sat,'Chla/raster/',paste('rasterStack',idi,
sep='_'),'.grd',sep='')
writeRaster(mchla.rb.PELGASs,filei6,overwrite=TRUE)
idi=paste('sdSSChlaPELGAS',paste(year1PELGAS,yearNPELGAS,sep='-'),sep = "_")
filei7 = paste(path.export.grid.sat,'Chla/raster/',paste('rasterStack',idi,
sep='_'),'.grd',sep='')
writeRaster(sdchla.rb.PELGASs,filei7,overwrite=TRUE)
idi=paste('q90SSChlaPELGAS',paste(year1PELGAS,yearNPELGAS,sep='-'),sep = "_")
filei8 = paste(path.export.grid.sat,'Chla/raster/',paste('rasterStack',idi,
sep='_'),'.grd',sep='')
writeRaster(q90chla.rb.PELGASs,filei8,overwrite=TRUE)
# 4.2. SST ----------
#**********************
# 3.2.1. Eventually, import SSChla raste stack ------
# 4.2.1. Eventually, import SSChla raste stack ------
#echosonde='/run/user/1000/gvfs/smb-share:domain=ifr,server=echosonde,share=data,user=mdoray/'
#donnees2='/run/user/1000/gvfs/smb-share:domain=ifr,server=nantes,share=donnees2,user=mdoray/'
#path.export.grid.sat=paste(echosonde,'Satellite/Atlantic/',sep='')
......@@ -492,7 +556,7 @@ plot(sst.rb)
# sst.rb=brick(sst.rasterStack): very long
# 3.2.2. Aggregate Chla raster stacks ------------
# Aggregate SST raster stacks
sst.rs.times=data.frame(rb.names=names(sst.rb))
sst.rs.times$year=as.numeric(substr(sst.rs.times$rb.names,2,5))
sst.rs.times$jday=as.numeric(substr(sst.rs.times$rb.names,7,9))
......@@ -513,14 +577,14 @@ sst.rs.aDates6d=aggregate(sst.rs.times$time,by=list(sst.rs.times$group6d),mean)
sst.rs.aDates30d=aggregate(sst.rs.times$time,by=list(sst.rs.times$month,sst.rs.times$year),mean,na.rm=TRUE)
sst.rs.aDates30d$myid=paste(sst.rs.aDates30d$Group.2,sst.rs.aDates30d$Group.1,sep='-')
# weekly averages: much faster
# 4.2.2. Weekly averages ------------
sst.rb.6d=aggregate(sst.rb,fact=c(1,1,6))
sst.rb
sst.rb.6d
# correct layer names
names(sst.rb.6d)=substr(format(sst.rs.aDates6d[,2]),1,10)
# monthly averages:
# 4.2.3. Monthly averages ------------
lmyid=unique(sst.rs.times$myid)
# or easier, with stackApply:
msst.rb.yearmonth=stackApply(sst.rb,indices=as.numeric(factor(sst.rs.times$myid)),fun=mean,na.rm=TRUE)
......@@ -534,11 +598,70 @@ plot(sdsst.rb.yearmonth,col=gradient.colors)
# or, with 30 days bins
sst.rb.monthyear=aggregate(sst.rb,fact=c(1,1,30))
# 3.2.3. Save mean sst raster stacks to disk ---------
# 4.2.4. Average during PELGAS surveys ---------
# *********************************************
# Define PELGAS intervals
sst.rs.times$inPELGAS=cut(sst.rs.times$time,breaks=PELGAS.voyages.timingsf,
labels=FALSE)
unique(sst.rs.times$inPELGAS)
# Average chla bricks according to PELGAS intervals
msst.rb.PELGAScut=stackApply(sst.rb,
indices=as.numeric(factor(sst.rs.times$inPELGAS)),
fun=mean,na.rm=TRUE)
sdsst.rb.PELGAScut=stackApply(sst.rb,
indices=as.numeric(factor(sst.rs.times$inPELGAS)),
fun=sd,na.rm=TRUE)
q90sst.rb.PELGAScut=stackApply(sst.rb,
indices=as.numeric(factor(sst.rs.times$inPELGAS)),
fun=quantile90,na.rm=TRUE)
# Select mean chla layers corresponding to PELGAS intervals
msst.rb.PELGAS=subset(msst.rb.PELGAScut,names(msst.rb.PELGAScut)[
substr(names(msst.rb.PELGAScut),7,8)%in%
seq(1,max(sst.rs.times$inPELGAS,na.rm = TRUE),2)])
names(msst.rb.PELGAS)=paste('PELGAS',seq(2000,2019,1),sep='')
plot(msst.rb.PELGAS)
sdsst.rb.PELGAS=subset(sdsst.rb.PELGAScut,names(sdsst.rb.PELGAScut)[
substr(names(sdsst.rb.PELGAScut),7,8)%in%
seq(1,max(sst.rs.times$inPELGAS,na.rm = TRUE),2)])
names(sdsst.rb.PELGAS)=paste('PELGAS',seq(2000,2019,1),sep='')
plot(sdsst.rb.PELGAS)
q90sst.rb.PELGAS=subset(q90sst.rb.PELGAScut,names(q90sst.rb.PELGAScut)[
substr(names(q90sst.rb.PELGAScut),7,8)%in%
seq(1,max(sst.rs.times$inPELGAS,na.rm = TRUE),2)])
names(q90sst.rb.PELGAS)=paste('PELGAS',seq(2000,2019,1),sep='')
plot(q90sst.rb.PELGAS)
# Select cells in PELGAS polygon
msst.rb.PELGAS0 <- crop(msst.rb.PELGAS, PELGASpolygon.sppol)
msst.rb.PELGASs <- mask(msst.rb.PELGAS0,PELGASpolygon.sppol)
plot(msst.rb.PELGASs)
sdsst.rb.PELGAS0 <- crop(sdsst.rb.PELGAS, PELGASpolygon.sppol)
sdsst.rb.PELGASs <- mask(sdsst.rb.PELGAS0,PELGASpolygon.sppol)
# 4.2.5. Save mean SST raster stacks to disk ---------
#************************************************
idi=paste('mSSTYearMonth',paste(min(year),max(year),sep='-'),sep = "_")
filei6 = paste(path.export.grid.sat,'SST/raster/',paste('rasterStack',idi,sep='_'),'.grd',sep='')
writeRaster(msst.rb.yearmonth,filei6,overwrite=TRUE)
idi=paste('sdSSTYearMonth',paste(min(year),max(year),sep='-'),sep = "_")
filei7 = paste(path.export.grid.sat,'SST/raster/',paste('rasterStack',idi,sep='_'),'.grd',sep='')
writeRaster(sdsst.rb.yearmonth,filei7,overwrite=TRUE)
idi=paste('mSSTYearMonth',paste(min(years.sst),max(years.sst),sep='-'),
sep = "_")
filei9 = paste(path.export.grid.sat,'SST/raster/',paste('rasterStack',idi,
sep='_'),'.grd',sep='')
writeRaster(msst.rb.yearmonth,filei9,overwrite=TRUE)
idi=paste('sdSSTYearMonth',paste(min(years.sst),max(years.sst),sep='-'),
sep = "_")
filei10 = paste(path.export.grid.sat,'SST/raster/',paste('rasterStack',idi,
sep='_'),'.grd',sep='')
writeRaster(sdsst.rb.yearmonth,filei10,overwrite=TRUE)
year1PELGAS=2000;yearNPELGAS=2019
idi=paste('mSSTPELGAS',paste(year1PELGAS,yearNPELGAS,sep='-'),
sep = "_")
filei11 = paste(path.export.grid.sat,'SST/raster/',paste('rasterStack',idi,
sep='_'),'.grd',sep='')
writeRaster(msst.rb.PELGASs,filei11,overwrite=TRUE)
idi=paste('sdSSTPELGAS',paste(year1PELGAS,yearNPELGAS,sep='-'),sep = "_")
filei12 = paste(path.export.grid.sat,'SST/raster/',paste('rasterStack',idi,
sep='_'),'.grd',sep='')
writeRaster(sdsst.rb.PELGASs,filei12,overwrite=TRUE)
......@@ -12,18 +12,19 @@ echosonde='/run/user/1000/gvfs/smb-share:domain=ifr,server=echosonde,share=data,
donnees2='/run/user/1000/gvfs/smb-share:domain=ifr,server=nantes,share=donnees2,user=mdoray/'
# path to raster maps on echosonde drive
prefix='Y:/'
prefix='/media/mathieu/IfremerMData/'
path.export.grid.sat=paste(prefix,'Satellite/Atlantic/',sep='')
year=1998:2021
years.chla=1998:2021
# import Chl-a
idi=paste('mSSChlaYearMonth',paste(min(year),max(year),sep='-'),sep = "_")
idi=paste('mSSChlaYearMonth',paste(min(years.chla),max(years.chla),sep='-'),sep = "_")
filei3 = paste(path.export.grid.sat,'Chla/raster/',paste('rasterStack',idi,sep='_'),'.grd',sep='')
idi=paste('sdSSChlaYearMonth',paste(min(year),max(year),sep='-'),sep = "_")
idi=paste('sdSSChlaYearMonth',paste(min(years.chla),max(years.chla),sep='-'),sep = "_")
filei4 = paste(path.export.grid.sat,'Chla/raster/',paste('rasterStack',idi,sep='_'),'.grd',sep='')
idi=paste('q90SSChlaYearMonth',paste(min(year),max(year),sep='-'),sep = "_")
idi=paste('q90SSChlaYearMonth',paste(min(years.chla),max(years.chla),sep='-'),sep = "_")
filei5 = paste(path.export.grid.sat,'Chla/raster/',paste('rasterStack',idi,sep='_'),'.grd',sep='')
# read in raster bricks
# Read in raster bricks per month and year
mchla.rb.yearmonth=brick(filei3)
sdchla.rb.yearmonth=brick(filei4)
q90chla.rb.yearmonth=brick(filei5)
......@@ -35,16 +36,44 @@ lmyid=unique(chla.rs.times$myid)
q98mchla.rb.yearmonth=cellStats(mchla.rb.yearmonth,stat=quantile98,na.rm=TRUE)
summary(q98mchla.rb.yearmonth)
# Read in raster bricks per PELGAS survey
year1PELGAS=2000;yearNPELGAS=2019
idi=paste('mSSChlaPELGAS',paste(year1PELGAS,yearNPELGAS,sep='-'),
sep = "_")
filei6 = paste(path.export.grid.sat,'Chla/raster/',paste('rasterStack',idi,
sep='_'),'.grd',sep='')
mchla.rb.PELGAS=brick(filei6)
idi=paste('sdSSChlaPELGAS',paste(year1PELGAS,yearNPELGAS,sep='-'),sep = "_")
filei7 = paste(path.export.grid.sat,'Chla/raster/',paste('rasterStack',idi,
sep='_'),'.grd',sep='')
sdchla.rb.PELGAS=brick(filei7)
idi=paste('q90SSChlaPELGAS',paste(year1PELGAS,yearNPELGAS,sep='-'),sep = "_")
filei8 = paste(path.export.grid.sat,'Chla/raster/',paste('rasterStack',idi,
sep='_'),'.grd',sep='')
q90chla.rb.PELGAS=brick(filei8)
# import SST
idi=paste('mSSTYearMonth',paste(min(year),max(year),sep='-'),sep = "_")
years.sst=2000:2021
idi=paste('mSSTYearMonth',paste(min(years.sst),max(years.sst),sep='-'),sep = "_")
filei6 = paste(path.export.grid.sat,'SST/raster/',paste('rasterStack',idi,sep='_'),'.grd',sep='')
idi=paste('sdSSTYearMonth',paste(min(year),max(year),sep='-'),sep = "_")
idi=paste('sdSSTYearMonth',paste(min(years.sst),max(years.sst),sep='-'),sep = "_")
filei7 = paste(path.export.grid.sat,'SST/raster/',paste('rasterStack',idi,sep='_'),'.grd',sep='')
# read in raster bricks
msst.rb.yearmonth=brick(filei6)
sdsst.rb.yearmonth=brick(filei7)
plot(msst.rb.yearmonth,col=gradient.colors)
year1PELGAS=2000;yearNPELGAS=2019
idi=paste('mSSTPELGAS',paste(year1PELGAS,yearNPELGAS,sep='-'),
sep = "_")
filei11 = paste(path.export.grid.sat,'SST/raster/',paste('rasterStack',idi,
sep='_'),'.grd',sep='')
msst.rb.PELGAS=brick(filei11)
idi=paste('sdSSTPELGAS',paste(year1PELGAS,yearNPELGAS,sep='-'),sep = "_")
filei12 = paste(path.export.grid.sat,'SST/raster/',paste('rasterStack',idi,
sep='_'),'.grd',sep='')
sdsst.rb.PELGAS=brick(filei12)
# 4. Time series analysis ----
#****************************
# 4.1. Chla ----
......@@ -54,7 +83,6 @@ plot(msst.rb.yearmonth,col=gradient.colors)
# Rectangular areas
mchla.rb.yearmonth.dp <- crop(mchla.rb.yearmonth, extent(-10.7,8,35.8,51))
mchla.rb.yearmonth.BoB <- crop(mchla.rb.yearmonth, extent(-9,0,43,49))
mchla.rb.yearmonth.BoB <- crop(mchla.rb.yearmonth, extent(-9,0,43,49))
mchla.rb.yearmonth.GoL <- crop(mchla.rb.yearmonth, extent(3,8,41,45))
mchla.rb.yearmonth.atl <- crop(mchla.rb.yearmonth, extent(-10.7,0,35.8,52))
mchla.rb.yearmonth.ES <- crop(mchla.rb.yearmonth, extent(-3,-2.7,47,47.3))
......@@ -69,23 +97,30 @@ mchla.rb.yearmonth.PELGAS0 <- crop(mchla.rb.yearmonth, PELGASpolygon.sppol)
mchla.rb.yearmonth.PELGAS <- mask(mchla.rb.yearmonth.PELGAS0,
PELGASpolygon.sppol)
# Mean SSChl-a map, 2000-2018
# 4.1.0. Mean SSChl-a map, > 1998 ----
plot(mean(mchla.rb.yearmonth,na.rm=TRUE),
main='Mean surface Chl-a concentration, 2000-2018');coast()
main=paste('Mean surface Chl-a concentration,',
paste(min(years.chla),max(years.chla),sep='-')));coast()
plot(mean(mchla.rb.yearmonth.dp,na.rm=TRUE),
main='Mean surface Chl-a concentration, 2000-2018');coast()
main=paste('Mean surface Chl-a concentration,',
paste(min(years.chla),max(years.chla),sep='-')));coast()
plot(mean(mchla.rb.yearmonth.BoB,na.rm=TRUE),
main='Mean surface Chl-a concentration, 2000-2018');coast()
main=paste('Mean surface Chl-a concentration,',
paste(min(years.chla),max(years.chla),sep='-')));coast()
plot(mean(mchla.rb.yearmonth.GoL,na.rm=TRUE),
main='Mean surface Chl-a concentration, 2000-2018');coast()
main=paste('Mean surface Chl-a concentration,',
paste(min(years.chla),max(years.chla),sep='-')));coast()
plot(mean(mchla.rb.yearmonth.atl,na.rm=TRUE),
main='Mean surface Chl-a concentration, 2000-2018');coast()
main=paste('Mean surface Chl-a concentration,',
paste(min(years.chla),max(years.chla),sep='-')));coast()
plot(mean(mchla.rb.yearmonth.ES,na.rm=TRUE),
main='Mean surface Chl-a concentration, 2000-2018');coast()
main=paste('Mean surface Chl-a concentration,',
paste(min(years.chla),max(years.chla),sep='-')));coast()
plot(mean(mchla.rb.yearmonth.PELGAS,na.rm=TRUE),
main='Mean surface Chl-a concentration, 2000-2018');coast()
main=paste('Mean surface Chl-a concentration,',
paste(min(years.chla),max(years.chla),sep='-')));coast()
# Mean SSChl-a time series, 2000-2018
# Mean SSChl-a time series, >1998
# 4.1.1. Defipel area ----
# 4.1.1.1. Chl-a seasonal cycle -------
mchla.yearmonth.dp=data.frame(yearMonth=names(mchla.rb.yearmonth.dp),
......@@ -102,17 +137,18 @@ h + geom_smooth(aes(group = 1), size = 2, method = "auto", se = TRUE)
# 4.1.1.1. Chl-a time series -------
msschla.defipel=cellStats(mchla.rb.yearmonth.dp,mean,na.rm=TRUE)
msschla.defipel.ts = ts(msschla.defipel, start = c(2000,1), frequency = 12)
msschla.defipel.ts = ts(msschla.defipel, start = c(1998,1), frequency = 12)
length(msschla.defipel)
grep('X2010.1',names(msschla.defipel))
msschla.defipel.ts2010 = ts(msschla.defipel[121:228], start = c(2010,1), frequency = 12)
year1=grep('X2010.1',names(msschla.defipel))[1]
msschla.defipel.ts2010 = ts(msschla.defipel[year1:length(msschla.defipel)],
start = c(2010,1), frequency = 12)
plot(msschla.defipel.ts2010,ylab='Mean surface Chl-a concentration')
# Time series decomposition
components.ts.defipel = decompose(msschla.defipel.ts)
plot(components.ts.defipel)
# Trend analysis
library(forecast)
components.ts.defipel.trend=ts(c(components.ts.defipel$trend), start = c(2000,1), frequency = 12)
components.ts.defipel.trend=ts(c(components.ts.defipel$trend), start = c(1998,1), frequency = 12)
plot(components.ts.defipel.trend)
ts.fit.sschla.defipel <- tslm(msschla.defipel.ts ~ trend + season)
summary(ts.fit.sschla.defipel)
......@@ -124,7 +160,10 @@ resy=c(residuals(ts.fit.sschla.defipel))
qqnorm(resy);qqline(resy, col = 2)
plot(msschla.defipel.ts2010)
lines(ts.fit.sschla.defipel2010$fitted.values,col=2)
# -> -0.0112032 mg.L-1 Chla / year since 2010 in Defipel area: -0.1344384 total
# -> -0.0007804 mg.L-1 Chla / year since 2010 in Defipel area
-0.0007804*(2021-2010) # total decrease
# Seawifs precision: +/- 35% ...
# Half cnfidence interval range
mean(msschla.defipel.ts2010)*0.35
# -> non significant linear trend
......@@ -132,18 +171,21 @@ mean(msschla.defipel.ts2010)*0.35
# 4.1.2.1. Chl-a seasonal cycle -------
mchla.yearmonth.BoB=data.frame(yearMonth=names(mchla.rb.yearmonth.BoB),
year=substr(names(mchla.rb.yearmonth.BoB),2,5),
month=substr(names(mchla.rb.yearmonth.BoB),7,8),
month=as.numeric(substr(names(mchla.rb.yearmonth.BoB),7,8)),
mChla=cellStats(mchla.rb.yearmonth.BoB,stat=mean))
mchla.yearmonth.BoB$month=factor(mchla.yearmonth.BoB$month,levels=seq(12),ordered=TRUE)
mchla.yearmonth.BoB$month=factor(mchla.yearmonth.BoB$month,
levels=seq(12),ordered=TRUE)
h=ggplot(mchla.yearmonth.BoB, aes(month,mChla,colour=year))+
geom_line(aes(group= year))+ggtitle("Bay of Biscay area")
h + geom_smooth(aes(group = 1), size = 2, method = "auto", se = TRUE)
# 4.1.1.1. Chl-a time series -------
msschla.BoB=cellStats(mchla.rb.yearmonth.BoB,mean,na.rm=TRUE)
msschla.BoB.ts = ts(msschla.BoB, start = c(2000,1), frequency = 12)
msschla.BoB.ts = ts(msschla.BoB, start = c(1998,1), frequency = 12)
plot(msschla.BoB.ts,ylab='Mean surface Chl-a concentration')
msschla.BoB.ts2010 = ts(msschla.BoB[121:228], start = c(2010,1), frequency = 12)
year1=grep('X2010.1',names(msschla.BoB))[1]
msschla.BoB.ts2010 = ts(msschla.BoB[year1:length(msschla.BoB)],
start = c(2010,1), frequency = 12)
plot(msschla.BoB.ts2010,ylab='Mean surface Chl-a concentration')
# Time series decomposition
components.BoB.ts = decompose(msschla.BoB.ts)
......@@ -161,7 +203,8 @@ hist(resy)
qqnorm(resy);qqline(resy, col = 2)
plot(msschla.BoB.ts2010)
lines(ts.fit.sschla.defipel2010$fitted.values,col=2)
# -> -0.0103884 mg.L-1 Chla / year since 2010 in BoB area
# -> -0.0006458 mg.L-1 Chla / year since 2010 in BoB area
# Seawifs precision: +/- 35% ...
mean(msschla.BoB.ts2010)*0.35
# -> non significant linear trend
......@@ -169,9 +212,10 @@ mean(msschla.BoB.ts2010)*0.35
# 4.1.3.1. Chl-a seasonal cycle -------
mchla.yearmonth.GoL=data.frame(yearMonth=names(mchla.rb.yearmonth.GoL),
year=substr(names(mchla.rb.yearmonth.GoL),2,5),
month=substr(names(mchla.rb.yearmonth.GoL),7,8),
month=as.numeric(substr(names(mchla.rb.yearmonth.GoL),7,8)),
mChla=cellStats(mchla.rb.yearmonth.GoL,stat=mean))
mchla.yearmonth.GoL$month=factor(mchla.yearmonth.GoL$month,levels=seq(12),ordered=TRUE)
mchla.yearmonth.GoL$month=factor(mchla.yearmonth.GoL$month,levels=seq(12),
ordered=TRUE)
h=ggplot(mchla.yearmonth.GoL, aes(month,mChla,colour=year))+
geom_line(aes(group= year))+ggtitle("Gulf of Lion area")
h + geom_smooth(aes(group = 1), size = 2, method = "auto", se = TRUE)
......@@ -179,7 +223,9 @@ h + geom_smooth(aes(group = 1), size = 2, method = "auto", se = TRUE)
msschla.GoL=cellStats(mchla.rb.yearmonth.GoL,mean,na.rm=TRUE)
msschla.GoL.ts = ts(msschla.GoL, start = c(2000,1), frequency = 12)
plot(msschla.GoL.ts,ylab='Mean surface Chl-a concentration')
msschla.GoL.ts2010 = ts(msschla.GoL[121:228], start = c(2010,1), frequency = 12)
year1=grep('X2010.1',names(msschla.GoL))[1]
msschla.GoL.ts2010 = ts(msschla.GoL[year1:length(msschla.GoL)],
start = c(2010,1), frequency = 12)
plot(msschla.GoL.ts2010,ylab='Mean surface Chl-a concentration')
# Time series decomposition
components.GoL.ts = decompose(msschla.GoL.ts)
......@@ -197,7 +243,7 @@ hist(resy)
qqnorm(resy);qqline(resy, col = 2)
plot(msschla.GoL.ts2010)
lines(ts.fit.sschla.defipel2010$fitted.values,col=2)
# -> -0.0116244 mg.L-1 Chla / year since 2010 in GoL area: -0.0929952 total
# -> -0.0005090 mg.L-1 Chla / year since 2010 in GoL area: -0.0929952 total
# Seawifs precision: +/- 35% ...
mean(msschla.GoL.ts2010)*0.35
# -> non significant linear trend
......@@ -205,7 +251,7 @@ mean(msschla.GoL.ts2010)*0.35
# 4.1.3. EchoSonde area -----
mchla.yearmonth.ES=data.frame(yearMonth=names(mchla.rb.yearmonth.ES),
year=substr(names(mchla.rb.yearmonth.ES),2,5),
month=substr(names(mchla.rb.yearmonth.ES),7,8),
month=as.numeric(substr(names(mchla.rb.yearmonth.ES),7,8)),
mChla=cellStats(mchla.rb.yearmonth.ES,stat=mean))
mchla.yearmonth.ES$month=factor(mchla.yearmonth.ES$month,levels=seq(12),ordered=TRUE)
h=ggplot(mchla.yearmonth.ES, aes(month,mChla,colour=year))+
......@@ -213,10 +259,16 @@ h=ggplot(mchla.yearmonth.ES, aes(month,mChla,colour=year))+
h + geom_smooth(aes(group = 1), size = 2, method = "auto", se = TRUE)
msschla.ES=cellStats(mchla.rb.yearmonth.ES,mean,na.rm=TRUE)
msschla.ES.ts = ts(msschla.ES, start = c(2000,1), frequency = 12)
msschla.ES.ts = ts(msschla.ES, start = c(1998,1), frequency = 12)
plot(msschla.ES.ts,ylab='Mean surface Chl-a concentration')
msschla.ES.ts2010 = ts(msschla.ES[121:228], start = c(2010,1), frequency = 12)
year1=grep('X2010.1',names(msschla.ES))[1]
msschla.ES.ts2010 = ts(msschla.ES[year1:length(msschla.ES)],
start = c(2010,1), frequency = 12)
plot(msschla.ES.ts2010,ylab='Mean surface Chl-a concentration')
year1=grep('X2017.1',names(msschla.ES))[1]
msschla.ES.ts2017 = ts(msschla.ES[year1