Commit 9ef930b9 authored by DORAY's avatar DORAY
Browse files

faster raster update

parent b0d0ebf0
......@@ -28,9 +28,9 @@ pref='G:/'
pref='/media/mathieu/IfremerMData/'
# on network
pref=echosonde
pref='X:/'
# Rasters export path
path.export.grid.sat=paste(pref,'Satellite/Data/Atlantic/',sep='')
path.export.grid.sat=paste(pref,'Satellite/Atlantic/',sep='')
# 1. Define grid and polygon selection
# to resample sat images on gridmap grid (0.25°)----------
......@@ -92,12 +92,12 @@ coastlsp <- map2SpatialPolygons(coastlm, IDs=coastlm$names,proj4string=CRS("+pro
### 2. Create Chl-a raster stack --------
# **********************************
# 2.1. Define path
patnam=paste(pref,"Satellite/Data/Atlantic/Chla/Rdata/",sep='')
patnam=paste(pref,"Satellite/Atlantic/Chla/Rdata/",sep='')
# patout = path to temporary folder
patout='/users/mathieu/Documents/temp/'
patout='/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/'
path.grids1='/Documents/temp/'
# Years range
lyears.Chla=unique(substr(ficnam,1,4))
......@@ -110,6 +110,12 @@ lofiles.chla=list.files(paste(path.export.grid.sat,'/Chla/raster/',sep=''))
# 2.2. Convert to SSChla raster stacks (3-days composites) -----
#******************************************
if (exists('sschla.rasterStack')){
lidik=names(sschla.rasterStack)
}else{
lidik=NULL
}
jj=seq(from=1,to=366,by=3)
year=seq(min(lyears.Chla),max(lyears.Chla))
tabmax=1000 # chla max value
......@@ -122,37 +128,45 @@ for (i in 1:length(year)) {
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.rasterStack=rasteri2m
idik=paste('X',year[i],'.',kkar,sep='')
if (!idik%in%lidik){
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.rasterStack=rasteri2m
}else{
#chla.df=rbind(chla.df,tablo.long)
sschla.rasterStack=stack(sschla.rasterStack,rasteri2m)
}
cat('Chla image imported','\n')
}else{
#chla.df=rbind(chla.df,tablo.long)
sschla.rasterStack=stack(sschla.rasterStack,rasteri2m)
cat('Chla image already imported','\n')
}
}else{
cat('Chla image does not exists','\n')
}
}
}
......@@ -168,7 +182,9 @@ filei2 = paste(path.export.grid.sat,'Chla/raster/',
writeRaster(sschla.rasterStack,filei2,overwrite=TRUE)
cat("Raster stack saved as:", filei2, "\n")
rm('sschla.rasterStack')
#sschla.rasterStack=stack(filei2)
#rm('sschla.rasterStack')
# # Resample by week
# plot(sschla.rasterStack)
......@@ -192,15 +208,21 @@ rm('sschla.rasterStack')
# 3. SST rasters (3-days composites) ----
#******************************************
# 2.1. Define paths
patnam2=paste(pref,"Satellite/Data/Atlantic/SST/Rdata/",sep='')
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)
#******************************************
if (exists('sst.rasterStack')){
lidik2=names(sst.rasterStack)
}else{
lidik2=NULL
}
jj=seq(from=1,to=366,by=1)
jj=seq(from=131,to=366,by=1)
year=2000:2018
year=2000:2021
#year=c(2017,2017)
tabmax=1000 # SST max value
for (i in 1:length(year)) {
......@@ -211,39 +233,47 @@ for (i in 1:length(year)) {
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(patnam2,ficnam2[selfic],sep=""))
dim(sst)
sst=sst[ilon,jlat] #SST
sst[sst>tabmax]=tabmax
# Convert to raster
rasteri=raster(t(sst))
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('SST',paste(year[i],kkar,sep='-')))
# Resample to gridmaps resolution
rasteri2=resample(rasteri,rae, fun='bilinear')
rasteri2m=mask(rasteri2,coastlsp,inverse=TRUE)
names(rasteri2m)
# 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
sst.rasterBrick=rasteri2m
idik2=paste('X',year[i],'.',kkar,sep='')
if (!idik2%in%lidik2){
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(patnam2,ficnam2[selfic],sep=""))
dim(sst)
sst=sst[ilon,jlat] #SST
sst[sst>tabmax]=tabmax
# Convert to raster
rasteri=raster(t(sst))
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('SST',paste(year[i],kkar,sep='-')))
# Resample to gridmaps resolution
rasteri2=resample(rasteri,rae, fun='bilinear')
rasteri2m=mask(rasteri2,coastlsp,inverse=TRUE)
names(rasteri2m)
# 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
sst.rasterStack=rasteri2m
}else{
#chla.df=rbind(chla.df,tablo.long)
sst.rasterStack=stack(sst.rasterStack,rasteri2m)
}
cat('SST image imported','\n')
}else{
#chla.df=rbind(chla.df,tablo.long)
sst.rasterBrick=brick(sst.rasterBrick,rasteri2m)
cat('SST image already imported','\n')
}
}else{
cat('SST image does not exists','\n')
}
}
}
......@@ -252,15 +282,16 @@ for (i in 1:length(year)) {
#sst.rasterStack2017.2=sst.rasterStack
# bind 2 rasterstacks
#sst.rasterStack=stack(sst.rasterStack2017.2,sst.rasterStack2018)
rev(names(sst.rasterBrick))
rev(names(sst.rasterStack))
#rm(sst.rasterStack)
# 2.3. Save SST raster stack to disk ---------
# (~200 Mo, 18 year series, 15 km resolution, 3 days composite)
idi=paste('SST',paste(min(year),max(year),sep='-'),sep = "_")
#idi=paste('SST',paste(min(year),2017,sep='-'),sep = "_")
filei2 = paste(path.export.grid.sat,'SST/raster/',
paste('rasterBrick',idi,sep='_'),'.grd',sep='')
writeRaster(sst.rasterBrick,filei2,overwrite=TRUE)
paste('rasterStack',idi,sep='_'),'.grd',sep='')
writeRaster(sst.rasterStack,filei2,overwrite=TRUE)
# import raster
echosonde='/run/user/1000/gvfs/smb-share:domain=ifr,server=echosonde,share=data,user=mdoray/'
......@@ -282,7 +313,7 @@ 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:2019
year=1998:2021
idi=paste('SSChla',paste(min(year),max(year),sep='-'),sep = "_")
filei2 = paste(path.export.grid.sat,'Chla/raster/',
paste('rasterStack',idi,sep='_'),'.grd',sep='')
......@@ -321,11 +352,14 @@ chla.rs.times$time=strptime(paste(chla.rs.times$year,
chla.rs.times$month,chla.rs.times$day,
sep='-'),format="%Y-%m-%d",tz="GMT")
chla.rs.times$group6d
if (is.integer)
group6d0=rep(seq(floor(length(chla.rs.times$time)/2)),each=2)
group6d0=rep(seq(ceiling(length(chla.rs.times$time)/2)),each=2)
length(group6d0)
length(chla.rs.times$time)]
chla.rs.aDates6d=aggregate(chla.rs.times$time,by=list(chla.rs.times$group6d),
length(chla.rs.times$time)
chla.rs.times$group6d=group6d0[1:length(chla.rs.times$time)]
head(chla.rs.times)
chla.rs.aDates6d=aggregate(chla.rs.times$time,
by=list(chla.rs.times$group6d),
mean)
chla.rs.aDates30d=aggregate(chla.rs.times$time,
by=list(chla.rs.times$month,chla.rs.times$year),
......@@ -442,7 +476,7 @@ names(mchla.rb.yearmonth)=lmyid
# 3.1.3. Save mean Chla raster stacks to disk ---------
#************************************************
year=2000:2018
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='')
writeRaster(mchla.rb.yearmonth,filei3,overwrite=TRUE)
......@@ -462,7 +496,7 @@ chla.rb.monthyear=aggregate(chla.rb,fact=c(1,1,30))
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:2018
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='')
# read in a single raster layer
......@@ -498,17 +532,15 @@ 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
# weekly averages: much faster
sst.rb.6d=aggregate(sst.rb,fact=c(1,1,6))
# much faster
sst.rb
sst.rb.6d
dim(sst.rs.aDates)
# correct layer names
names(sst.rb.6d)=substr(format(sst.rs.aDates6d[,2]),1,10)
# monthly averages: to be done by hand in a loop...
# lmyid=unique(sst.rs.times$myid)
#
# for (i in 1:length(lmyid)){
# print(lmyid[i])
# m=sst.rs.times$myid==lmyid[i]
......@@ -529,6 +561,7 @@ names(sst.rb.6d)=substr(format(sst.rs.aDates6d[,2]),1,10)
# names(msst.rb.yearmonth)=lmyid
# names(sdsst.rb.yearmonth)=lmyid
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)
names(msst.rb.yearmonth)=lmyid
......
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