Commit 772bd9a9 authored by md0276b's avatar md0276b
Browse files

initial sat scripts commit

parent d7e4d3e3
# Animation cartes par an
# GG
# Auteur: Pierre Petitgas
library(fields)
library(maps)
library(mapdata)
library(EchoR)
#load("gradient_colors.rda")
#source('isobath.r')
# lon lat
patnam='c:/ptg/work/Atlantic/'
# westernmost_longitude: -12.000
# easternmost_longitude: 13.000
# southernmost_latitude: 36.000
# northernmost_latitude: 60.000
# spatial_resolution: 0.015 degrees
# load(paste(patnam,"lat_lon.data",sep=""))
lon=seq(-12,13,0.015);length(lon)
lat=seq(36,60,length.out = 2401);length(lat)
# zone GG
x1=-6; x2=-1; y1=43; y2=48
x1=-12; x2=13; y1=36; y2=52
ilon=(1:length(lon))[lon>(x1) & lon<(x2)]
jlat=(1:length(lat))[lat>y1 & lat<y2]
asp=1.2
###
### chloro
###
echosonde='/run/user/1000/gvfs/smb-share:domain=ifr,server=echosonde,share=data,user=mdoray/'
patnam="c:/ptg/work/Atlantic/Chloro/"
patnam=paste(echosonde,"Satellite/Atlantic/Chla/",sep='')
patout="c:/ptg/work/Atlantic/fig_chloro_gg/"
patout='/users/mathieu/Documents/temp/'
ficnam=list.files(path=patnam) # lire les noms de fichier
ficnam=ficnam[substr(ficnam,14,17)=="data"] #chloro
# figures
jj=seq(from=1,to=366,by=3)
year=1999:2015
tabmax=20
for (i in 1:length(year)) {
cat(paste(year[i]),"\n")
for (k in 1:length(jj)) {
selfic= substr(ficnam,1,4)==paste(year[i]) & as.numeric(substr(ficnam,10,12))==jj[k] # chloro
if (sum(selfic)!=0) {
load(paste(patnam,ficnam[selfic],sep=""))
dim(tablo)
length(lon);length(lat)
tablo=tablo[ilon,jlat] #chloro
tablo[tablo>tabmax]=tabmax
if (jj[k]<10) {kkar=paste("00",jj[k],sep="")}
if (jj[k]>=10 & k<100) {kkar=paste("0",jj[k],sep="")}
if (jj[k]>=100) {kkar=paste(jj[k])}
image.plot(tablo)
image.plot(lon,lat,tablo,xlab="",ylab="",asp=asp,
col=gradient.colors,zlim=c(0,tabmax))
image.plot(lon[ilon],lat[jlat],tablo,xlab="",ylab="",asp=asp,
col=gradient.colors,zlim=c(0,tabmax))
title(paste(year[i],kkar,sep="-"))
map(database="worldHires",col=gray(.8),fill=T,bg="white",add=T)
box() #trace un cadre autour de la carte
map.axes() #trace des axes en bas et ? gauche
if (k<10) {nn=paste("00",k,sep="")}
if (k>=10 & k<100) {nn=paste("0",k,sep="")}
if (k>=100) {nn=paste(k)}
dev.print(png, file=paste(patout,"chloro_",year[i],"_",nn,".png",sep=""),width=500, height=500)
}
}
}
###
### SST
###
patnam='c:/ptg/work/Atlantic/SST/'
patout="c:/ptg/work/Atlantic/fig_sst_gg/"
ficnam=list.files(path=patnam) # lire les noms de fichier
ficnam=ficnam[substr(ficnam,17,20)=="data"] #sst
# figures
jj=seq(from=1,to=366,by=3)
year=1999:2015
tabmax=24; tabmin=8
for (i in 1:length(year)) {
cat(paste(year[i]),"\n")
for (k in 1:length(jj)) {
selfic= substr(ficnam,4,7)==paste(year[i]) & as.numeric(substr(ficnam,13,15))==jj[k] # sst
if (sum(selfic)!=0) {
load(paste(patnam,ficnam[selfic],sep=""))
tablo=sst[ilon,jlat] #sst
tablo[tablo<0]=NA; tablo[tablo>tabmax]=tabmax
if (jj[k]<10) {kkar=paste("00",jj[k],sep="")}
if (jj[k]>=10 & k<100) {kkar=paste("0",jj[k],sep="")}
if (jj[k]>=100) {kkar=paste(jj[k])}
image.plot(lon[ilon],lat[jlat],tablo,xlab="",ylab="",asp=asp,
col=gradient.colors,zlim=c(tabmin,tabmax))
title(paste(year[i],kkar,sep="-"))
map(database="worldHires",col=gray(.8),fill=T,bg="white",add=T)
box() #trace un cadre autour de la carte
map.axes() #trace des axes en bas et ? gauche
if (k<10) {nn=paste("00",k,sep="")}
if (k>=10 & k<100) {nn=paste("0",k,sep="")}
if (k>=100) {nn=paste(k)}
dev.print(png, file=paste(patout,"sst_",year[i],"_",nn,".png",sep=""),
width=500, height=500)
}
}
}
library(R.utils)
library(tiff)
library(fields)
conv.mj2jj.f=function(mm,jc,year) {
# conversion (jour calendaire, mois) -> jour julien
jm=c(31,28,31,30,31,30,31,31,30,31,30,31)
if ( is.element(year,seq(2000,2016,by=4)) ) {
jm=c(31,29,31,30,31,30,31,31,30,31,30,31)
}
if (mm>1) { jj=jc+cumsum(jm)[(mm-1)] }
else {jj=jc}
return(jj)
}
conv.n2char3.f=function(nn) {
if (nchar(nn)==1) {cn=paste("00",nn,sep="")}
else if (nchar(nn)==2) {cn=paste("0",nn,sep="")}
else {cn=paste(nn)}
return(cn)
}
############## read1_SST ######
### B: c'est le disque "/home8/begmeil/gwel_public/atlantic/
### diff?rents produits selon les ann?es: Pathfinder, Ostia, Odyssea
### depuis 04/2017 B: "/iota1/satcoast-public/atlantic/"
iota1="/run/user/1000/gvfs/smb-share:domain=ifr,server=iota1,share=satcoast-public,user=mdoray/"
#
path.atlantic=paste(iota1,'atlantic/',sep='')
#path.atlantic=iota1
# dossier o? je vais stocker les r?sultats
patout="C:/ptg/WORK/Atlantic/SST/"
patout="/users/mathieu/Documents/Data/SST/"
# 1999-0101 ? 2009-1231: Pathfinder
for(annee in 1999:2009) {
patin=paste(path.atlantic,"pathfinder/",annee,"/SST/",sep="")
fic=dir(path=patin,pattern="tif.gz")
fic2=dir(path=patin,pattern="png")
if(length(fic)>=1) {
for (k in 1:length(fic)) {
fichier=paste(patin,fic[k],sep="")
fichier=paste(patin,"'",fic[k],"'",sep="")
fichier2=paste(patin,fic2[k],sep="")
sortie=paste(patout,gsub("[.]gz$", "", fic[k]),sep="")
#isGzipped(fichier)
#gunzip(filename=fichier, destname=sortie, overwrite=T, remove=F, BFR.SIZE=1e+07)
#gunzip(filename=fichier, destname=sortie,overwrite=T)
system(paste('cd',patin,'&&','cp',fic[k],patout))
try(system(paste('cp ',paste(fichier,sep=''),patout)))
system(paste('cp ',fichier2,patout))
try(system(paste('gunzip',paste(patout,fic[k],sep=''))))
# Je modifie le fichier pour qu'il ait la m?me structure que celui de chloro
mat=readTIFF(source=sortie); mat[mat<0]=NA
mat=t(mat)
nl=dim(mat)[1]; nc=dim(mat)[2]
sst=matrix(NA,nrow=nl,ncol=nc)
for (j in 1:nc) { sst[,j]=mat[,(nc+1)-j] }
# Je sauve le fichier dans un .data et d?truit le .tif
ye=as.numeric(substr(fic[k],4,7))
mo=as.numeric(substr(fic[k],8,9)); jo=as.numeric(substr(fic[k],10,11))
joju=conv.mj2jj.f(mo,jo,ye)
joju=conv.n2char3.f(joju)
ficout=paste(patout,substr(fic[k],1,11),"-",joju,".data",sep="")
save(sst,file=ficout)
file.remove(sortie)
}
}
}
# 2010-0101 ? 2010-0831 : Ostia
annee=2010
patout="C:/ptg/WORK/Atlantic/SST/ostia/"
patin=paste("B:/ostia/",annee,"/SST/",sep="")
fic=dir(path=patin,pattern="tif.gz")
if(length(fic)>=1) {
for (k in 1:length(fic)) {
fichier=paste(patin,fic[k],sep="")
sortie=paste(patout,gsub("[.]gz$", "", fic[k]),sep="")
gunzip(filename=fichier, destname=sortie, overwrite=T, remove=F, BFR.SIZE=1e+07)
# Je modifie le fichier pour qu'il ait la m?me structure que celui de chloro
mat=readTIFF(source=sortie); mat[mat<0]=NA
mat=t(mat)
nl=dim(mat)[1]; nc=dim(mat)[2]
sst=matrix(NA,nrow=nl,ncol=nc)
for (j in 1:nc) { sst[,j]=mat[,(nc+1)-j] }
# Je sauve le fichier dans un .data et d?truit le .tif
ye=as.numeric(substr(fic[k],4,7))
mo=as.numeric(substr(fic[k],8,9)); jo=as.numeric(substr(fic[k],10,11))
joju=conv.mj2jj.f(mo,jo,ye)
joju=conv.n2char3.f(joju)
ficout=paste(patout,substr(fic[k],1,11),"-",joju,".data",sep="")
save(sst,file=ficout)
file.remove(sortie)
}
}
# 2010-0901 ? 2013-1231 : Odyssea
patout="C:/ptg/WORK/Atlantic/SST/odyssea/"
for(annee in 2010:2013) {
patin=paste("B:/odyssea/",annee,"/SST/",sep="")
fic=dir(path=patin,pattern="tif.gz")
if(length(fic)>=1) {
for (k in 1:length(fic)) {
fichier=paste(patin,fic[k],sep="")
sortie=paste(patout,gsub("[.]gz$", "", fic[k]),sep="")
gunzip(filename=fichier, destname=sortie, overwrite=T, remove=F, BFR.SIZE=1e+07)
# Je modifie le fichier pour qu'il ait la m?me structure que celui de chloro
mat=readTIFF(source=sortie); mat[mat<0]=NA
mat=t(mat)
nl=dim(mat)[1]; nc=dim(mat)[2]
sst=matrix(NA,nrow=nl,ncol=nc)
for (j in 1:nc) { sst[,j]=mat[,(nc+1)-j] }
# Je sauve le fichier dans un .data et d?truit le .tif
ye=as.numeric(substr(fic[k],4,7))
mo=as.numeric(substr(fic[k],8,9)); jo=as.numeric(substr(fic[k],10,11))
joju=conv.mj2jj.f(mo,jo,ye)
joju=conv.n2char3.f(joju)
ficout=paste(patout,substr(fic[k],1,11),"-",joju,".data",sep="")
save(sst,file=ficout)
file.remove(sortie)
}
}
}
# 2014-0101 ? 2015-1231 : Odyssea
patout="C:/ptg/WORK/Atlantic/SST/"
for(annee in 2014:2015) {
patin=paste("B:/odyssea/",annee,"/SST/",sep="")
fic=dir(path=patin,pattern="tif.gz")
if(length(fic)>=1) {
for (k in 1:length(fic)) {
fichier=paste(patin,fic[k],sep="")
sortie=paste(patout,gsub("[.]gz$", "", fic[k]),sep="")
gunzip(filename=fichier, destname=sortie, overwrite=T, remove=F, BFR.SIZE=1e+07)
# Je modifie le fichier pour qu'il ait la m?me structure que celui de chloro
mat=readTIFF(source=sortie); mat[mat<0]=NA
mat=t(mat)
nl=dim(mat)[1]; nc=dim(mat)[2]
sst=matrix(NA,nrow=nl,ncol=nc)
for (j in 1:nc) { sst[,j]=mat[,(nc+1)-j] }
# Je sauve le fichier dans un .data et d?truit le .tif
ye=as.numeric(substr(fic[k],4,7))
mo=as.numeric(substr(fic[k],8,9)); jo=as.numeric(substr(fic[k],10,11))
joju=conv.mj2jj.f(mo,jo,ye)
joju=conv.n2char3.f(joju)
ficout=paste(patout,substr(fic[k],1,11),"-",joju,".data",sep="")
save(sst,file=ficout)
file.remove(sortie)
}
}
}
# 2016-0101 ? 2016-1231 : Odyssea
patout="C:/ptg/WORK/Atlantic/SST/"
patout="/users/mathieu/Documents/Data/SST/"
iota1='/users/mathieu/.gvfs/smb-share:domain=ifr,server=iota1,share=satcoast-public,user=mdoray/'
donnees22="/users/mathieu/.gvfs/smb-share:domain=ifr,server=nantes,share=donnees2,user=mdoray/"
iota1=paste()
#read.table(paste(donnees22,'Campagnes/PELGAS/Data/Acoustique/PELGAS0012_sAbiom_esdus.txt',sep=''),sep=';')
for(annee in 2016:2018) {
patin=paste(iota1,"atlantic/odyssea/",annee,"/SST/",sep="")
fic=dir(path=patin,pattern="tif.gz")
if(length(fic)>=1) {
for (k in 1:length(fic)) {
fichier=paste(patin,fic[k],sep="")
sortie=paste(patout,gsub("[.]gz$", "", fic[k]),sep="")
gunzip(filename=fichier, destname=sortie, overwrite=T, remove=F, BFR.SIZE=1e+07)
# Je modifie le fichier pour qu'il ait la m?me structure que celui de chloro
mat=readTIFF(source=sortie); mat[mat<0]=NA
mat=t(mat)
nl=dim(mat)[1]; nc=dim(mat)[2]
sst=matrix(NA,nrow=nl,ncol=nc)
for (j in 1:nc) { sst[,j]=mat[,(nc+1)-j] }
# Je sauve le fichier dans un .data et d?truit le .tif
ye=as.numeric(substr(fic[k],4,7))
mo=as.numeric(substr(fic[k],8,9)); jo=as.numeric(substr(fic[k],10,11))
joju=conv.mj2jj.f(mo,jo,ye)
joju=conv.n2char3.f(joju)
ficout=paste(patout,substr(fic[k],1,11),"-",joju,".data",sep="")
save(sst,file=ficout)
file.remove(sortie)
}
}
}
library(R.utils)
library(ncdf4)
# Dossier o? sont les donn?es source
# setwd("/home8/begmeil/gwel/public/atlantic/analysis/cdf/")
# setwd("//begmeil/gwel_public/atlantic/analysis/cdf/")
iota1='/run/user/1000/gvfs/smb-share:domain=ifr,server=iota1,share=satcoast-public,user=mdoray/'
setwd("//iota1/satcoast-public/atlantic/analysis/cdf/")
setwd(paste(iota1,"atlantic/analysis/cdf/",sep=''))
# dossier o? je vais stocker les r?sultats
dossier_out="C:/ptg/WORK/Atlantic/Chloro/"
dossier_out='/users/mathieu/Documents/temp/'
# je vais de 1998 ?? 2012
for(annee in 1998:2012)
{
# Je cherche tous les sous-dossiers qui commencent par 0, 1, 2 ou 3
for(jour in dir(path=as.character(annee), pattern="^[0-3]"))
{
print(paste("Annee : ",annee," Jour : ", jour))
# je cherche le nom du fichier compress?
fic=dir(file.path(annee, jour),pattern=".bz2")
# S'il y en a un...
if(length(fic==1))
{
# Je construis son chemin complet
fichier=file.path(annee,jour, fic)
# Cr?ation du nom du fichier d?compress?
sortie=file.path(dossier_out,gsub("[.]bz2$", "", fic))
# Je d?compresse
bunzip2(fichier, destname=sortie, overwrite=T, remove=FALSE, BFR.SIZE=1e+07)
tab1 <- nc_open(sortie)
# Je lis la variable analysed_chl_a
tablo=ncvar_get(tab1,"analysed_chl_a")
nc_close(tab1)
# et je sauve la variable que j'ai lue
save(tablo,file=file.path(dossier_out,paste(substr(fic,1,9), jour,".data",sep="")))
file.remove(sortie)
}
}
}
# je vais de 2013 ?? 2015
for(annee in 2013:2015)
{
# Je cherche tous les sous-dossiers qui commencent par 0, 1, 2 ou 3
for(jour in dir(path=as.character(annee), pattern="^[0-3]"))
{
print(paste("Annee : ",annee," Jour : ", jour))
# je cherche le nom du fichier compress?
fic=dir(file.path(annee, jour),pattern=".bz2")
# S'il y en a un...
if(length(fic==1))
{
# Je construis son chemin complet
fichier=file.path(annee,jour, fic)
# Cr?ation du nom du fichier d?compress?
sortie=file.path(dossier_out,gsub("[.]bz2$", "", fic))
# Je d?compresse
bunzip2(fichier, destname=sortie, overwrite=T, remove=FALSE, BFR.SIZE=1e+07)
tab1 <- nc_open(sortie)
# Je lis la variable analysed_chl_a
tablo=ncvar_get(tab1,"analysed_chl_a")
nc_close(tab1)
# et je sauve la variable que j'ai lue
save(tablo,file=file.path(dossier_out,paste(substr(fic,1,9), jour,".data",sep="")))
file.remove(sortie)
}
}
}
# je vais de 2016 ?? 2016
for(annee in 2016:2016)
{
# Je cherche tous les sous-dossiers qui commencent par 0, 1, 2 ou 3
for(jour in dir(path=as.character(annee), pattern="^[0-3]"))
{
print(paste("Annee : ",annee," Jour : ", jour))
# je cherche le nom du fichier compress?
fic=dir(file.path(annee, jour),pattern=".bz2")
# S'il y en a un...
if(length(fic==1))
{
# Je construis son chemin complet
fichier=file.path(annee,jour, fic)
# Cr?ation du nom du fichier d?compress?
sortie=file.path(dossier_out,gsub("[.]bz2$", "", fic))
# Je d?compresse
bunzip2(fichier, destname=sortie, overwrite=T, remove=FALSE, BFR.SIZE=1e+07)
tab1 <- nc_open(sortie)
# Je lis la variable analysed_chl_a
tablo=ncvar_get(tab1,"analysed_chl_a")
nc_close(tab1)
# et je sauve la variable que j'ai lue
save(tablo,file=file.path(dossier_out,paste(substr(fic,1,9), jour,".data",sep="")))
file.remove(sortie)
}
}
}
# Convert satellite maps to raster maps
# Auteur: Pierre Petitgas, Mathieu Doray
# ******************************************
library(fields)
library(maps)
library(mapdata)
library(EchoR)
library(chron)
library(ggplot2)
library(rgdal)
library(raster)
library(rasterVis)
library(maptools)
data(gradient_colors)
paletteGC=rasterTheme(region=gradient.colors)
# 1. Define grid and polygon selection ----------
#*********************************
#patnam='c:/ptg/work/Atlantic/'
# westernmost_longitude: -12.000
# easternmost_longitude: 13.000
# southernmost_latitude: 36.000
# northernmost_latitude: 60.000
# spatial_resolution: 0.015 degrees
# load(paste(patnam,"lat_lon.data",sep=""))
lon=seq(-12,13,0.015);length(lon)
lat=seq(36,60,length.out = 2401);length(lat)
# Default WGACEGG grid
xa1<-(-10.2); xa2<-(-1); ya1<-35.8; ya2<-48
# Define cell dimensions
ax=0.25;ay=0.25
define.grid.poly.mask(xa1,xa2,ya1,ya2,ax,ay,shelf=TRUE)
# 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))
res(rae)=c(0.25,0.25)
# 1.3. Eventually, area of interest selection -------
x1=-12; x2=13; y1=36; y2=60
ilon=(1:length(lon))[lon>(x1) & lon<(x2)]
jlat=(1:length(lat))[lat>y1 & lat<y2]
asp=1.2
# 1.4. Import Defipel polygon --------
DefipelPolygon.shp <- readOGR(dsn = paste(donnees2,"Satellite",sep=''), layer = "Defipel_Polygon")
plot(DefipelPolygon.shp)
DefipelPolygon.shpll=spTransform(DefipelPolygon.shp,"+proj=longlat +datum=WGS84")
DefipelPolygon.df=fortify(DefipelPolygon.shpll)
head(DefipelPolygon.df)
xpol=DefipelPolygon.df$long
ypol=DefipelPolygon.df$lat
plot(matxg,matyg,main='Grid and polygon',
xlab='',ylab='',type='p',col='grey50')
coast()
#add polygon
sel=inout(pts=cbind(as.vector(matxg),as.vector(matyg)),
poly=cbind(xpol,ypol),bound=T,quiet=T)
lines(xpol,ypol)
coast()
points(matxg[sel],matyg[sel],pch=16)
legend('topleft',c('Discarded','Selected'),pch=c(1,16))
# 1.5. Import coastline as polygons -------
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. Chl-a rasters --------
# **********************************
# 2.1. Define paths -----
# echosonde = define path to 'data' folder on 'echosonde' drive
echosonde='/run/user/1000/gvfs/smb-share:domain=ifr,server=echosonde,share=data,user=mdoray/'
patnam=paste(echosonde,"Satellite/Atlantic/Chla/data/",sep='')
# patout = path to temporary folder
patout='/users/mathieu/Documents/temp/'
ficnam=list.files(path=patnam) # lire les noms de fichier
ficnam=ficnam[substr(ficnam,14,17)=="data"] #chloro
path.grids1='/users/mathieu/Documents/temp/'
# donnees2 = path to 'donnees2' drive
donnees2='/run/user/1000/gvfs/smb-share:domain=ifr,server=nantes,share=donnees2,user=mdoray/'
# path.export = path to raster maps output folder
path.export.grid.sat=paste(echosonde,'Satellite/Atlantic/',sep='')
# 2.2. Convert to SSChla raster stacks (3-days composites) -----
#******************************************
jj=seq(from=1,to=366,by=3)
year=2000:2018
tabmax=1000 # chla max value
for (i in 1:length(year)) {
for (k in 1:length(jj)) {
cat(paste(year[i],jj[k]),"\n")
selfic= substr(ficnam,1,4)==paste(year[i]) & as.numeric(substr(ficnam,10,12))==jj[k] # chloro
if (sum(selfic)!=0) {
if (jj[k]<10) {kkar=paste("00",jj[k],sep="")}
if (jj[k]>=10 & k<100) {kkar=paste("0",jj[k],sep="")}
if (jj[k]>=100) {kkar=paste(jj[k])}
mdy=month.day.year(as.numeric(kkar),origin=c(month = 1, day = 1, year = year[i]))
if (mdy$month<10){monthi=paste('0',mdy$month,sep='')}
if (mdy$day<10){dayi=paste('0',mdy$day,sep='')}
load(paste(patnam,ficnam[selfic],sep=""))
dim(tablo)
length(lon);length(lat)
tablo=tablo[ilon,jlat] #chloro
tablo[tablo>tabmax]=tabmax
# Convert to raster
rasteri=raster(t(tablo))
xmin(rasteri)=min(lon);xmax(rasteri)=max(lon)
ymin(rasteri)=min(lat);ymax(rasteri)=max(lat)
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
rasteri2=resample(rasteri,rae, fun='bilinear')
rasteri2m=mask(rasteri2,coastlsp,inverse=TRUE)
# par(mfrow=c(1,2))
# plot(rasteri,main='Max resolution')
# coast()
# plot(rasteri2m,main='15km resolution')
# coast()
# par(mfrow=c(1,1))
if (i==1&k==1){
#chla.df=tablo.long
sschla.rasterBrick=rasteri2m
}else{
#chla.df=rbind(chla.df,tablo.long)
sschla.rasterBrick=brick(sschla.rasterBrick,rasteri2m)
}
}
}
}
# 2.3. Save Chla raster stack to disk ---------
# (~200 Mo, 18 year series, 15 km resolution, 3 days composite)
idi=paste('SSChla',paste(min(year),max(year),sep='-'),sep = "_")
filei2 = paste(path.export.grid.sat,'Chla/raster/',paste('rasterBrick',idi,sep='_'),'.grd',sep='')
writeRaster(sschla.rasterBrick,filei2,overwrite=TRUE)
cat("Raster brick saved as:", filei2, "\n")
# Resample by week
plot(sschla.rasterStack)
rasteri2=aggregate(sschla.rasterStack)
# Add mean an stdev maps to stack
mr=calc(r,mean,na.rm=TRUE)
names(mr)=paste('Avg',paste(min(lyears),max(lyears),sep='-'))
# or avgMap2=calc(r, fun=mean)
omr=cellStats(mr,'mean')
if (dim(r)[3]>1){
stdr=calc(r, fun=sd,na.rm=TRUE)+omr
}else{
stdr=mr
values(stdr)=NA