Commit 4a06822c authored by DORAY's avatar DORAY
Browse files

sat script clean up

parent ea536ffa
......@@ -128,9 +128,12 @@ if (exists('sschla.rasterStack')){
}
# Recreate stack from scratch if problem in export
#rm('sschla.rasterStack')
#lidik=NULL
#year=seq(1998,2021,1)
#
rm('sschla.rasterStack')
#
lidik=NULL
#
year=seq(1998,2021,1)
jj=seq(from=1,to=366,by=3)
year=seq(min(lyears.Chla),max(lyears.Chla))
......@@ -332,8 +335,9 @@ library(chron)
#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='')
year=1998:2021
idi=paste('SSChla',paste(min(year),max(year),sep='-'),sep = "_")
years.chla=1998:2021
idi=paste('SSChla',paste(min(years.chla),max(years.chla),
sep='-'),sep = "_")
filei2 = paste(path.export.grid.sat,'Chla/raster/',
paste('rasterStack',idi,sep='_'),'.grd',sep='')
# read in a single raster layer
......@@ -457,7 +461,7 @@ PELGAS.voyages.timings$enddatef=strptime(
PELGAS.voyages.timingsf=sort(strptime(union(PELGAS.voyages.timings$startdatef,
PELGAS.voyages.timings$enddatef),
format="%Y-%m-%d",tz="GMT"))
length(PELGAS.voyages.timingsf);length(PELGAS.voyages.timingsf.labels)
length(PELGAS.voyages.timingsf)
# Define PELGAS intervals
chla.rs.times$inPELGAS=cut(chla.rs.times$time,breaks=PELGAS.voyages.timingsf,
......@@ -550,9 +554,12 @@ writeRaster(q90chla.rb.PELGASs,filei8,overwrite=TRUE)
#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='')
year=2000:2021
idi=paste('SST',paste(min(year),max(year),sep='-'),sep = "_")
filei2 = paste(path.export.grid.sat,'SST/raster/',paste('rasterStack',idi,sep='_'),'.grd',sep='')
years.sst=2000:2021
idi=paste('SST',paste(min(years.sst),
max(years.sst),sep='-'),sep = "_")
filei2 = paste(path.export.grid.sat,'SST/raster/',
paste('rasterStack',idi,sep='_'),'.grd',
sep='')
# read in a single raster layer
#chla.r=raster(filei2)
#plot(chla.r)
......@@ -567,24 +574,39 @@ plot(sst.rb)
# 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))
mdy=month.day.year(sst.rs.times$jday,origin=c(month = 1, day = 1, year = sst.rs.times$year))
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))
mdy=month.day.year(
sst.rs.times$jday,origin=c(month = 1, day = 1,
year = sst.rs.times$year))
sst.rs.times$month=mdy$month
sst.rs.times$day=mdy$day
head(sst.rs.times)
sst.rs.times$smonth=as.character(sst.rs.times$month)
sst.rs.times[nchar(sst.rs.times$smonth)==1,'smonth']=paste('0',
sst.rs.times[nchar(sst.rs.times$smonth)==1,'smonth'],
sep='')
sst.rs.times$myid=paste(sst.rs.times$year,sst.rs.times$smonth,sep='-')
sst.rs.times$date=as.Date(paste(sst.rs.times$year,sst.rs.times$month,sst.rs.times$day,sep='-'))
sst.rs.times$time=strptime(paste(sst.rs.times$year,sst.rs.times$month,sst.rs.times$day,sep='-'),
format="%Y-%m-%d",tz="GMT")
sst.rs.times$group6d=rep(seq(length(sst.rs.times$time)/6),each=6)
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='-')
sst.rs.times[nchar(sst.rs.times$smonth)==1,'smonth']=
paste('0',sst.rs.times[nchar(sst.rs.times$smonth)==1,
'smonth'],sep='')
sst.rs.times$myid=paste(sst.rs.times$year,
sst.rs.times$smonth,sep='-')
sst.rs.times$date=as.Date(
paste(sst.rs.times$year,sst.rs.times$month,
sst.rs.times$day,sep='-'))
sst.rs.times$time=strptime(
paste(sst.rs.times$year,sst.rs.times$month,
sst.rs.times$day,sep='-'),
format="%Y-%m-%d",tz="GMT")
sst.rs.times$group6d=rep(seq(length(sst.rs.times$time)/6),
each=6)
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='-')
# 4.2.2. Weekly averages ------------
sst.rb.6d=aggregate(sst.rb,fact=c(1,1,6))
......@@ -596,9 +618,13 @@ names(sst.rb.6d)=substr(format(sst.rs.aDates6d[,2]),1,10)
# 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)
msst.rb.yearmonth=stackApply(
sst.rb,indices=as.numeric(factor(sst.rs.times$myid)),
fun=mean,na.rm=TRUE)
names(msst.rb.yearmonth)=lmyid
sdsst.rb.yearmonth=stackApply(sst.rb,indices=as.numeric(factor(sst.rs.times$myid)),fun=sd,na.rm=TRUE)
sdsst.rb.yearmonth=stackApply(
sst.rb,indices=as.numeric(factor(sst.rs.times$myid)),
fun=sd,na.rm=TRUE)
names(sdsst.rb.yearmonth)=lmyid
plot(msst.rb.yearmonth,col=gradient.colors)
......@@ -610,67 +636,83 @@ sst.rb.monthyear=aggregate(sst.rb,fact=c(1,1,30))
# 4.2.4. Average during PELGAS surveys ---------
# *********************************************
# Define PELGAS intervals
sst.rs.times$inPELGAS=cut(sst.rs.times$time,breaks=PELGAS.voyages.timingsf,
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)),
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)),
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)),
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)[
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='')
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)[
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='')
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)[
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='')
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)
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(years.sst),max(years.sst),sep='-'),
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='')
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='-'),
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='')
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='')
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='')
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)
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