Commit b0d0ebf0 authored by md0276b's avatar md0276b
Browse files

polish

parent 13ec9e8e
......@@ -25,14 +25,15 @@ 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/'
# on external drive
pref='G:/'
pref='/media/mathieu/IfremerMData/'
# on network
pref=echosonde
# Rasters export path
path.export.grid.sat=paste(pref,'Satellite/Data/Atlantic/',sep='')
# 1. Define grid and polygon selection
# to resample sat images to gridmap resolution (0.25°)----------
# to resample sat images on gridmap grid (0.25°)----------
#*********************************
#patnam='c:/ptg/work/Atlantic/'
# westernmost_longitude: -12.000
......@@ -88,7 +89,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. Chl-a rasters --------
### 2. Create Chl-a raster stack --------
# **********************************
# 2.1. Define path
patnam=paste(pref,"Satellite/Data/Atlantic/Chla/Rdata/",sep='')
......@@ -97,6 +98,10 @@ 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/'
# Years range
lyears.Chla=unique(substr(ficnam,1,4))
# path.export = path to raster maps output folder
#pref='/media/mathieu/IfremerMData/'
# 2.1.1. Check file list in output folder ----------
......@@ -106,7 +111,7 @@ lofiles.chla=list.files(paste(path.export.grid.sat,'/Chla/raster/',sep=''))
# 2.2. Convert to SSChla raster stacks (3-days composites) -----
#******************************************
jj=seq(from=1,to=366,by=3)
year=2000:2021
year=seq(min(lyears.Chla),max(lyears.Chla))
tabmax=1000 # chla max value
for (i in 1:length(year)) {
for (k in 1:length(jj)) {
......@@ -143,49 +148,51 @@ for (i in 1:length(year)) {
# par(mfrow=c(1,1))
if (i==1&k==1){
#chla.df=tablo.long
sschla.rasterBrick=rasteri2m
sschla.rasterStack=rasteri2m
}else{
#chla.df=rbind(chla.df,tablo.long)
sschla.rasterStack=stack(sschla.rasterBrick,rasteri2m)
sschla.rasterStack=stack(sschla.rasterStack,rasteri2m)
}
}
}
}
# concert sschla.rasterStack to brick
sschla.rasterBrick=brick(sschla.rasterStack)
# sschla.rasterBrick=brick(sschla.rasterStack)
# 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
}
names(stdr)=paste('SD.',paste(min(lyears),max(lyears),'+',round(omr),sep='-'),sep='')
#r=stack(r,mr)
r=stack(r,mr,stdr)
paste('rasterStack',idi,sep='_'),'.grd',sep='')
writeRaster(sschla.rasterStack,filei2,overwrite=TRUE)
cat("Raster stack saved as:", filei2, "\n")
rm('sschla.rasterStack')
# # 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
# }
# names(stdr)=paste('SD.',paste(min(lyears),max(lyears),'+',round(omr),sep='-'),sep='')
# #r=stack(r,mr)
# r=stack(r,mr,stdr)
# 3. SST rasters (3-days composites) ----
#******************************************
# 2.1. Define paths
patnam2=paste(echosonde,"Satellite/Atlantic/SST/data/",sep='')
patnam2=paste(pref,"Satellite/Data/Atlantic/SST/Rdata/",sep='')
ficnam2=list.files(path=patnam2) # lire les noms de fichier
ficnam2=ficnam2[substr(ficnam2,17,20)=="data"] #SST
......@@ -251,7 +258,8 @@ rev(names(sst.rasterBrick))
# (~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='')
filei2 = paste(path.export.grid.sat,'SST/raster/',
paste('rasterBrick',idi,sep='_'),'.grd',sep='')
writeRaster(sst.rasterBrick,filei2,overwrite=TRUE)
# import raster
......@@ -270,22 +278,22 @@ library(raster)
library(chron)
# 3.1. Chla ----------
#**********************
# 3.1.1. Eventually, import SSChla raste stack ------
# 3.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(echosonde,'Satellite/Atlantic/',sep='')
year=2000:2018
path.export.grid.sat=paste(pref,'Satellite/Data/Atlantic/',sep='')
year=1998:2019
idi=paste('SSChla',paste(min(year),max(year),sep='-'),sep = "_")
filei2 = paste(path.export.grid.sat,'Chla/raster/',
paste('rasterStack',idi,sep='_'),'.grd',sep='')
# read in a single raster layer
chla.r=raster(filei2)
plot(chla.r)
# read in a raster stack
chla.rs=stack(filei2)
plot(chla.rs)
# read in a raster brick
rm('sschla.rasterStack')
# read in a raster stack (not efficient)
# chla.rs=stack(filei2)
# plot(chla.rs)
# rm('sschla.rasterStack')
# read in a raster brick: better
chla.rb=brick(filei2)
plot(chla.rb)
......@@ -293,9 +301,9 @@ plot(chla.rb)
chla.rb[chla.rb<=0]=NA
# 3.1.2. Aggregate Chla raster stacks ------------
chla.rs.times=data.frame(rs.names=names(chla.rs))
chla.rs.times$year=as.numeric(substr(chla.rs.times$rs.names,2,5))
chla.rs.times$jday=as.numeric(substr(chla.rs.times$rs.names,7,9))
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))
mdy=month.day.year(chla.rs.times$jday,
origin=c(month = 1, day = 1, year = chla.rs.times$year))
chla.rs.times$month=mdy$month
......@@ -312,7 +320,11 @@ chla.rs.times$date=as.Date(paste(chla.rs.times$year,
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=rep(seq(length(chla.rs.times$time)/2),each=2)
chla.rs.times$group6d
if (is.integer)
group6d0=rep(seq(floor(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),
mean)
chla.rs.aDates30d=aggregate(chla.rs.times$time,
......@@ -444,7 +456,7 @@ 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. SST ----------
# 3.2. Create SST raster stack----------
#**********************
# 3.2.1. Eventually, import SSChla raste stack ------
echosonde='/run/user/1000/gvfs/smb-share:domain=ifr,server=echosonde,share=data,user=mdoray/'
......@@ -466,9 +478,9 @@ plot(sst.rb)
# sst.rb=brick(sst.rasterStack): very long
# 3.2.2. Aggregate Chla raster stacks ------------
sst.rs.times=data.frame(rs.names=names(sst.rb))
sst.rs.times$year=as.numeric(substr(sst.rs.times$rs.names,2,5))
sst.rs.times$jday=as.numeric(substr(sst.rs.times$rs.names,7,9))
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$month=mdy$month
sst.rs.times$day=mdy$day
......
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