Commit 7260e55e authored by Mathieu DORAY's avatar Mathieu DORAY
Browse files
parents 0e559564 87fe8adc
......@@ -50,18 +50,26 @@ lat=seq(36,60,length.out = 2401);length(lat)
# Default WGACEGG grid
xa1<-(-10.2); xa2<-(-1); ya1<-35.8; ya2<-48
# WGACEGG grid compatible with PELGAS grid
xa1<-(-12); xa2<-(14); ya1<-35; ya2<-61
# Define cell dimensions
ax=0.25;ay=0.25
define.grid.poly.mask(xa1,xa2,ya1,ya2,ax,ay,shelf=TRUE)
xg # grid cell centers' longitudes
yg # grid cell centers' latitudes
# 1.2. Define extended WGACEGG grid raster ----
lxae=seq(-10.7,13.3,0.25)
lyae=seq(35.8,60.3,0.25)
rae=raster(nrow=length(lyae),ncol=length(lxae),xmn=min(lxae),xmx=max(lxae),ymn=min(lyae),ymx=max(lyae))
lxae=seq(-12,14,0.25)
lyae=seq(35,61,0.25)
rae=raster(nrow=length(lyae),ncol=length(lxae),xmn=min(lxae),xmx=max(lxae),
ymn=min(lyae),ymx=max(lyae))
res(rae)=c(0.25,0.25)
rae.df=as.data.frame(rae,xy=TRUE)
sort(unique(rae.df$x))
sort(unique(rae.df$y))
# 1.3. Eventually, area of interest selection -------
x1=-12; x2=13; y1=36; y2=60
x1=-12; x2=14; y1=35; y2=61
ilon=(1:length(lon))[lon>(x1) & lon<(x2)]
jlat=(1:length(lat))[lat>y1 & lat<y2]
asp=1.2
......@@ -153,7 +161,7 @@ for (i in 1:length(year)) {
rasteri=flip(rasteri, 2)
names(rasteri)=paste(year[i],kkar,sep='-')
#plot(rasteri,main=paste('SSChla',paste(year[i],kkar,sep='-')))
# Resample to gridmaps resolution
# 2.2.1. Resample to gridmaps resolution ----
rasteri2=resample(rasteri,rae, fun='bilinear')
rasteri2m=mask(rasteri2,coastlsp,inverse=TRUE)
# par(mfrow=c(1,2))
......@@ -261,7 +269,7 @@ for (i in 1:length(year)) {
rasteri=flip(rasteri, 2)
names(rasteri)=paste(year[i],kkar,sep='-')
#plot(rasteri,main=paste('SST',paste(year[i],kkar,sep='-')))
# Resample to gridmaps resolution
# 3.1.1. Resample to gridmaps resolution ----
rasteri2=resample(rasteri,rae, fun='bilinear')
rasteri2m=mask(rasteri2,coastlsp,inverse=TRUE)
names(rasteri2m)
......
......@@ -297,11 +297,30 @@ lines(ts.fit.sschla.ES2017$fitted.values,col=2)
components.ES2017.ts = decompose(msschla.ES.ts2017)
plot(components.ES2017.ts)
# -> no liear trend
# -> no linear trend
# Seawifs precision: +/- 35% ...
mean(msschla.ES.ts2010)*0.35
# -> non significant linear trend
# 4.1.4. PELGAS series -------
msschla.pelgas=cellStats(mchla.rb.PELGAS,mean,na.rm=TRUE)
msschla.pelgas.ts = ts(msschla.pelgas, start = c(2000,1), frequency = 1)
length(msschla.pelgas)
plot(msschla.pelgas.ts)
year1=grep('PELGAS2008',names(msschla.pelgas))[1]
msschla.pelgas.ts2008 = ts(msschla.pelgas[year1:length(msschla.pelgas)],
start = c(2008,1), frequency = 1)
plot(msschla.pelgas.ts2008,ylab='Mean surface Chl-a concentration')
lmChlaSat.PELGAS=tslm(msschla.pelgas.ts~trend)
summary(lmChlaSat.PELGAS)
lmChlaSat.PELGAS.ts2008=tslm(msschla.pelgas.ts2008~trend)
summary(lmChlaSat.PELGAS.ts2008)
# Seawifs precision: +/- 35% ...
mean(msschla.pelgas.ts2008)*0.35
# -> non significant linear trend
# 4.2. SST -----
#*********************************
# Select different areas
......@@ -426,9 +445,9 @@ h + geom_smooth(aes(group = 1), size = 2, method = "auto", se = TRUE)
# 4.2.3.2. SST time series -----
msssst.GoL=cellStats(msst.rb.yearmonth.GoL,mean,na.rm=TRUE)
msssst.GoL.ts = ts(msssst.GoL, start = c(2000,1), frequency = 12)
plot(msssst.GoL.ts,ylab='Mean surface Chl-a concentration')
plot(msssst.GoL.ts,ylab='Mean SST')
msssst.GoL.ts2010 = ts(msssst.GoL[121:228], start = c(2010,1), frequency = 12)
plot(msssst.GoL.ts2010,ylab='Mean surface Chl-a concentration')
plot(msssst.GoL.ts2010,ylab='Mean SST')
# Time series decomposition
components.GoL.ts = decompose(msssst.GoL.ts)
plot(components.GoL.ts)
......@@ -483,4 +502,22 @@ lines(ts.fit.sssst.defipel2010$fitted.values,col=2)
# -> +0.005844 °C / month = 0.070128°C / year since 2010 in ES area: +0.6°C total
# -> significant linear trend
# 4.1.5. PELGAS series -------
msst.pelgas=cellStats(msst.rb.PELGAS,mean,na.rm=TRUE)
msst.pelgas.ts = ts(msst.pelgas, start = c(2000,1), frequency = 1)
length(msst.pelgas)
plot(msst.pelgas.ts)
year1=grep('PELGAS2010',names(msst.pelgas))[1]
msst.pelgas.ts2010 = ts(msst.pelgas[year1:length(msst.pelgas)],
start = c(2010,1), frequency = 1)
plot(msst.pelgas.ts2010,ylab='Mean SST')
lmsstSat.PELGAS=tslm(msst.pelgas.ts~trend)
summary(lmsstSat.PELGAS)
lmsstSat.PELGAS.ts2010=tslm(msst.pelgas.ts2010~trend)
summary(lmsstSat.PELGAS.ts2010)
# Seawifs precision: +/- 35% ...
mean(msst.pelgas.ts2010)*0.35
# -> non significant linear trend
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment