source: CONFIG_DEVT/ORCHIDEE_OL_TP/ENSEMBLE/convert_nc_txt.R @ 5587

Last change on this file since 5587 was 5570, checked in by aclsce, 3 years ago

Created ORCHIDEE_OL_TP configuration : temporary configuration to be used during prectical session.
File size: 5.9 KB
Line 
1#!/usr/bin/env Rscript
2args = commandArgs(trailingOnly=TRUE)
3
4require(ncdf4)
5require(reshape2)
6require(plyr)
7library(dplyr)
8library(data.table)
9library(lubridate)
10
11print(args)
12
13##############################################################################################
14# initialisation variables
15#############################################################################################
16var <- c("NPP","GPP","IND","DIAMETER_DOM","HEIGHT_DOM","LAI","MOISTRESS","PLANT_STATUS","RDI","CCTRW_001","CCTRW_002","CCTRW_003")
17age=c(160,37,37,96,80,80,150,105,47,44,30,135,150,12,30,63,40,20,105,50,50)
18monthOrder <- c("janv.","févr.","mars","avril","mai","juin","juil.","août","sept.","oct.","nov.","déc.")
19LatOrderSite <- c('RU-Cok','FI-Kaa','FI-Sod','SE-Deg','SE-Fla','IS-Gun','FI-Hyy','RU-Zot','SE-Sk2','SE-SK1','SE-Nor','UK-Gri','DK-Fou','RU-Fyo','SE-Faj','UK-ESa','UK-EBu','UK-AMo','DK-Lva','DK-Ris','DK-Sor','RU-Ha2','RU-Ha1','RU-Ha3','NL-Lut','IE-Ca1','PL-wet','NL-Loo','NL-Hor','NL-Haa','IE-Dri','NL-Ca1','NL-Lan','UK-Her','NL-Mol','UK-PL3','BE-Bra','DE-Meh','UK-Tad','UK-Ham','DE-Geb','DE-Hai','DE-Tha','DE-Gri','DE-Kli','BE-Jal','BE-Lon','DE-Wet','BE-Vie','DE-Bay','CZ-BK1','SK-Tat','CZ-wet','FR-Gri','FR-Hes','FR-Fon','HU-Mat','CH-Lae','CH-Oe2','CH-Oe1','AT-Neu','HU-Bug','IT-Ren','IT-Mal','IT-MBo','IT-Lav','FR-Lq1','FR-Lq2','IT-LMa','IT-PT1','IT-Cas','FR-LBr','IT-Non','FR-Pue','IT-SRo','IT-Lec','IT-Pia','IT-Ro1','IT-Ro2','ES-VDA','IT-Amp','IT-Col','IT-Cpz','IT-Cp2','IT-BCi','ES-LMa','ES-ES1','ES-ES2','PT-Esp','PT-Mi1','PT-Mi2')
20
21#########################################################################################
22#   extract result from netcdf data
23##########################################################################################
24Table.Summary <- NULL
25h=0
26for (l in var) {
27    h <- h+1
28    j=0
29    Table.final <- NULL
30    for (i in LatOrderSite ) {
31        j <- j+1
32        file_path <- paste0(args[1],i,'_ALL.nc')
33        if (file.exists(file_path)==FALSE) next
34        else print(paste('OK pour :',i))
35        grid.nc <- nc_open(file_path)   
36        lon <- ncvar_get(grid.nc,"lon")
37        lat <- ncvar_get(grid.nc,"lat")
38        time <- ncvar_get(grid.nc,"time_counter")
39        if (args[5]=='D')tunits <- seq(as.Date("2000/1/1"), by = "day", length.out = length(time))
40        else tunits <- seq(as.Date("2000/1/1"), by = "month", length.out = length(time))
41        PFT  <- ncvar_get(grid.nc,"PFT")
42        tmp_array <- t(ncvar_get(grid.nc,l))   
43        tmp_array <- melt(tmp_array, value.name = l)
44        tmp_array[,"Date"] <- tunits[tmp_array[,"Var1"]]
45        tmp_array[,"Var1"] <- NULL
46        tmp_array[,"Lat"] <- round(lat,4)
47        tmp_array[,"Lon"] <- lon
48        tmp_array[,"Site"] <- substr(i,1,6)
49        if(j==1) { Table.final <- as.data.frame(tmp_array) }
50        else { Table.final <- bind_rows(Table.final,as.data.frame(tmp_array)) }   
51        nc_close(grid.nc)   
52    }
53    if(h==1) { Table.Summary <- Table.final }
54    else { Table.Summary <- left_join(Table.Summary,Table.final) }
55}
56if(args[5]=='D') {
57   Table.Summary$Day <- factor(format(Table.Summary$Date, "%d"))
58   Table.Summary$DOY <- factor(strftime(Table.Summary$Date, format = "%j")) 
59}
60Table.Summary$Month <- factor(format(Table.Summary$Date, "%b"),levels=monthOrder)
61Table.Summary$Year <- as.numeric(format(Table.Summary$Date , "%Y"))
62print(summary(Table.Summary))
63fwrite(Table.Summary,paste0(args[2],'Table.Summary_',args[7],'.txt'), sep=' ',quote=FALSE)
64############################################################################################
65
66monthOrder <- c("janv.","févr.","mars","avril","mai","juin","juil.","août","sept.","oct.","nov.","déc.")
67reloc <- time_length(interval(start = ymd("1970-01-01"),end = ymd("2000-01-01") ),unit = "seconds")
68climfile2 <- c("FR-Hes_1997-2006.nc","FI-Sod_2000-2006.nc","IT-Non_2001-2006.nc","DE-Bay_1996-1999.nc","DE-Wet_2002-2006.nc","SE-Nor_1996-2005.nc","SE-Fla_1996-2002.nc")
69var.f=c('sum','mean')
70var.o=c("GPPalt")
71var.oo=c("GPP")
72var.conv=c(1,1,1800)
73###########################################################################################
74# extract climate and GPP obs from netcdf data
75############################################################################################
76climate.final <- NULL
77for(l in c("FR","IT","DE","DK","FI","SE","NL","BE","CH","RU","CZ","OLD")){
78    if(l=='OLD') climfile=climfile2
79    else climfile=list.files('/home/surface5/nraoult/FLUXNET2015/done/',pattern = l)
80    for(i in climfile){
81    k=0;s=1
82    Site <- substr(i,1,6)
83    clim.t=data.frame(DOY=1:366)
84    if(l=='OLD') obs.nc <- nc_open(paste0('/home/orchideeshare/igcmg/IGCM/BC/OOL/OL2/FLUXNET/',i))
85    else obs.nc <- nc_open(paste0('/home/surface5/nraoult/FLUXNET2015/done/',i))
86    if(l=='OLD') var.final=var.oo
87    else var.final=var.o
88    for (j in var.final){
89        k=k+1
90        func=var.f[k]
91        conv=var.conv[k]
92        var.nc <- ncvar_get(obs.nc, j)
93        var.nc <- ifelse(var.nc>10e30,NA,var.nc)
94        lat <- ncvar_get(obs.nc,"lat")
95        lon <- ncvar_get(obs.nc,"lon")
96        time <- ncvar_get(obs.nc, "tstep")+ reloc
97        class(time) = c('POSIXt','POSIXct')
98        df <- data.frame(date=time,obs=var.nc*conv)
99        df$Month <- factor(format(df$date, "%b"),levels = monthOrder)
100        df$Year <- as.numeric(format(df$date , "%Y"))
101        df$DOY <- factor(strftime(df$date, format = "%j"))
102        ob <- tapply(df[,'obs'],list(df[,'DOY'],df[,'Year']),func)
103        clim.t <- data.frame(clim.t,apply(ob,1,max,na.rm=TRUE),
104            apply(ob,1,min,na.rm=TRUE),apply(ob,1,median,na.rm=TRUE))
105        names(clim.t)[c(s+1,s+2,s+3)]=c(paste0('GPP.max'),paste0('GPP.min'),paste0('GPP.med'))
106        s=s+3
107    }
108    clim.t$Site <- Site
109    clim.t$Lat <- round(lat,4)
110    clim.t$Lon <- round(lon,4)
111    climate.final <- rbind(climate.final,clim.t)
112    nc_close(obs.nc)     
113    }
114}
115#fwrite(climate.final,'climate_final.txt', sep=' ',quote=FALSE, na.string=NA)
116##########################################################################################
117
Note: See TracBrowser for help on using the repository browser.