00__creation_TestsData_v6.R 11 KB
Newer Older
unknown's avatar
Test  
unknown committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168

  DIR_WD <- "C:/Data/Hydrologie/Codes/airGR_dev/airGR_dev_0.7/";

  DIR_EXPORT_RDA  <- paste(DIR_WD,"files_data/",sep="");
  DIR_EXPORT_TXT  <- paste(DIR_WD,"files_misc/airGR_AdvancedExample/Test_DataInput/",sep="");

  ### DIR_DATA_INPUT_FR  <- "C:/Data/Hydrologie/BD/BD_FR_IRSTEA2010/BD_BASSINS/";  
  DIR_DATA_INPUT_FR  <- paste(DIR_WD,"tmp_data_input_DO_NOT_CIRCULATE/",sep="");
  ### DIR_DATA_INPUT_US  <- "C:/Data/Hydrologie/BD/BD_US_MOPEX/BD_MOPEX_430_DLY/";
  DIR_DATA_INPUT_US  <- paste(DIR_WD,"tmp_data_input_DO_NOT_CIRCULATE/",sep="");

  ListBasinCode <- c("J4902010","13337000"); ListBasinDB <- c("FR","US"); LIST_TMP <- list();
  ### ListBasinCode <- c("J4902010","12358500"); ListBasinDB <- c("FR","US"); LIST_TMP <- list();
  ### ListBasinCode <- c("M3340910","12358500"); ListBasinDB <- c("FR","US"); LIST_TMP <- list();
  for(iBasin in 1:length(ListBasinCode)){

    ##__IMPORT_________________________________________________________________________________####
      BasinCode <- NULL; BasinName <- NULL; BasinArea_km2_NEW <- NULL; HypsoData <- NULL;
      TabDatesT <- NULL; TabDatesR <- NULL; TabObsQls <- NULL; TabObsQmm <- NULL; TabObsP <- NULL; TabObsT <- NULL; TabObsE <- NULL; 
      if(ListBasinDB[iBasin]=="FR"){
        BasinCode <- ListBasinCode[iBasin];
        DIR_DATA_INPUT   <- DIR_DATA_INPUT_FR;  
        BasinCharactFile <- paste(DIR_DATA_INPUT,"00_Liste_BV.txt",sep=""); 
        HypsoDataFile    <- paste(DIR_DATA_INPUT,"_HypsoData_Pierre.txt",sep="");
        ##File_check
        DataImportFile <- paste(DIR_DATA_INPUT,BasinCode,"_BV.txt",sep="");
        FileExists <- file.exists(DataImportFile); if(FileExists==FALSE){ print("Error: DataImportFile not found",quote=FALSE); stop("EXECUTION STOPPED",call.=FALSE); }
        ##_Basin_Characteristics
        if(file.exists(BasinCharactFile)){ 
          TAB_CHARACT <- read.csv2(file=BasinCharactFile,header=TRUE,skip=0,stringsAsFactors=FALSE);
          IndBasin <- which(TAB_CHARACT$Code  == BasinCode);
          if(length(IndBasin)==1){ 
          BasinName         <- TAB_CHARACT$Nom[IndBasin]; 
          BasinArea_km2     <- as.numeric(TAB_CHARACT$S_CEM[IndBasin]);
          BasinArea_km2_NEW <- 10*round(BasinArea_km2/10);
          rm(TAB_CHARACT); } ### memory clear
        }
        if(file.exists(HypsoDataFile)){
          TAB_HYPSO <- read.csv2(file=HypsoDataFile,header=TRUE,skip=0,stringsAsFactors=FALSE);
          iR <- which(formatC(TAB_HYPSO$CODE,format="d",width=8,flag="0")==BasinCode);
          iC <- which(colnames(TAB_HYPSO)=="Zmin"):which(colnames(TAB_HYPSO)=="Zmax");
          HypsoData <- as.numeric(TAB_HYPSO[iR,iC]);  ### min, q01, q02, ... , q98, q99, max
          rm(TAB_HYPSO); ### memory clear
        }
        ##DataSeries
        Format <- c("A8","X1","I8","X1","A5","X1","A7","X1","F5.0","X1","F5.0","X1","F5.0","X1","F5.0","X1","F5.0");
        TAB_DATA <- read.fortran(file=DataImportFile,skip=51,header=FALSE,Format);
        TabDatesT  <- as.character(TAB_DATA[,1]);
        TabDatesR  <- as.POSIXlt(strptime(TAB_DATA[,1],format="%Y%m%d",tz="UTC"));
        TabObsQls  <- TAB_DATA[,2]; TabObsQls[TabObsQls<0] <- NA;   ### observed runoff (in l/s)
        TabObsQm3s <- TabObsQls/1000; ### observed runoff (in m3/s)
        TabObsQmm  <- TabObsQm3s*86.4/BasinArea_km2_NEW;  ### observed runoff (in mm/d)
        TabObsP    <- TAB_DATA[,5];   ### precipitation (catchment average in mm)
        TabObsT    <- TAB_DATA[,7];   ### air temp (catchment average in degre C)
        TabObsE    <- TAB_DATA[,8];   ### potential evap (catchment average in mm/d)
        rm(TAB_DATA); ### memory clear
      }
      if(ListBasinDB[iBasin]=="US"){
        BasinCode <- ListBasinCode[iBasin];
        DIR_DATA_INPUT   <- DIR_DATA_INPUT_US;
        BasinCharactFile <- paste(DIR_DATA_INPUT,"_List_US_416_LAURENT.txt",sep=""); 
        HypsoDataFile    <- paste(DIR_DATA_INPUT,"_HypsoData_431.txt",sep="");
        ##File_check
        DataImportFile <- paste(DIR_DATA_INPUT,BasinCode,".dly",sep="");
        FileExists <- file.exists(DataImportFile); if(FileExists==FALSE){ print("Error: DataImportFile not found",quote=FALSE); stop("EXECUTION STOPPED",call.=FALSE); }
        ##_Basin_Characteristics
        if(file.exists(BasinCharactFile)){ 
          TAB_CHARACT <- read.csv2(file=BasinCharactFile,header=TRUE,skip=0,stringsAsFactors=FALSE);
          IndBasin <- which(formatC(TAB_CHARACT[,1],format="d",width=8,flag="0")  == BasinCode);
          if(length(IndBasin)==1){ 
          BasinName         <- TAB_CHARACT[IndBasin,11];
          BasinArea_km2     <- as.numeric(TAB_CHARACT[IndBasin,10]);
          BasinArea_km2_NEW <- 10*round(BasinArea_km2/10);
          rm(TAB_CHARACT); } ### memory clear
        }
        if(file.exists(HypsoDataFile)){
          TAB_HYPSO <- read.csv2(file=HypsoDataFile,header=TRUE,skip=0,stringsAsFactors=FALSE);
          iR <- which(formatC(TAB_HYPSO$CODE,format="d",width=8,flag="0")==BasinCode);
          iC <- which(colnames(TAB_HYPSO)=="Zmin"):which(colnames(TAB_HYPSO)=="Zmax");
          HypsoData <- as.numeric(TAB_HYPSO[iR,iC]);  ### min, q01, q02, ... , q98, q99, max
          rm(TAB_HYPSO); ### memory clear
        }
        ##DataSeries
        Format <- c("A8","F10.0","F10.0","F10.0","F10.0","F10.0");
        TAB_DATA <- read.fortran(file=DataImportFile,skip=0,header=FALSE,Format);
        TabDatesT   <- TAB_DATA[,1]; TabDatesT <- gsub(pattern=" ",replacement="0",TabDatesT);
        TabDatesR   <- as.POSIXlt(strptime(TabDatesT,format="%Y%m%d",tz="UTC"));
        TabObsQmm   <- as.numeric(TAB_DATA[,4]); TabObsQmm[TabObsQmm<0] <- NA;   ### observed runoff (in mm)
        TabObsQls   <- TabObsQmm*BasinArea_km2_NEW/86.4*1000;   ### observed runoff (in l/s)
        TabObsP     <- as.numeric(TAB_DATA[,2]);   ### precipitation (catchment average in mm)
        TabObsTmin  <- as.numeric(TAB_DATA[,6]);   ### air temp min (catchment average in degre C)
        TabObsTmax  <- as.numeric(TAB_DATA[,5]);   ### air temp max (catchment average in degre C)
        TabObsT     <- (TabObsTmin+TabObsTmax)/2;  ### air temp (catchment average in degre C)
        TabObsE     <- as.numeric(TAB_DATA[,3]);   ### potential evap (catchment average in mm/d)
        rm(TAB_DATA); ### memory clear
      }
    

    ##__DATA_MODIF___and__BASINDATA_CREATION___________________________________________________####
      ##Name_changes_and_temporal_delay_to_mask_data_origin
      if(iBasin==1){
        BasinCode_NEW  <- "L0123001";
        BasinName_NEW  <- "Banjo River at Paterson Creek";
        Select_OLD     <- which(TabDatesR==as.POSIXlt(strptime("01/01/1972",format="%d/%m/%Y",tz="UTC"))):which(TabDatesR==as.POSIXlt(strptime("31/12/2000",format="%d/%m/%Y",tz="UTC")));
      }
      if(iBasin==2){
        BasinCode_NEW  <- "L0123002";
        BasinName_NEW  <- "Snowy River at Orroral Valley Homestead";
        Select_OLD     <- which(TabDatesR==as.POSIXlt(strptime("01/01/1972",format="%d/%m/%Y",tz="UTC"))):which(TabDatesR==as.POSIXlt(strptime("31/12/2000",format="%d/%m/%Y",tz="UTC")));
      }
      TabDatesR_NEW <- as.POSIXlt(seq(from=as.POSIXlt(strptime("01/01/1984",format="%d/%m/%Y",tz="UTC")),to=as.POSIXlt(strptime("31/12/2012",format="%d/%m/%Y",tz="UTC")),by="days"));
      if(identical(TabDatesR$mday[Select_OLD],TabDatesR_NEW$mday)==FALSE){ stop("STOP"); }
      ##BasinData
      BasinData               <- list();
      BasinData$BasinCode     <- BasinCode_NEW;
      BasinData$BasinName     <- BasinName_NEW; 
      BasinData$BasinArea_km2 <- BasinArea_km2_NEW;
      BasinData$HypsoData     <- HypsoData;
      BasinData$TabDatesT     <- format(TabDatesR_NEW,"%Y%m%d");
      BasinData$TabDatesR     <- TabDatesR_NEW;
      BasinData$TabObsQls     <- TabObsQls[Select_OLD];
      BasinData$TabObsQmm     <- TabObsQmm[Select_OLD];
      BasinData$TabObsP       <- TabObsP[Select_OLD];
      BasinData$TabObsT       <- TabObsT[Select_OLD];
      BasinData$TabObsE       <- TabObsE[Select_OLD];

      
    ##__SAVE___________________________________________________________________________________####
      LIST_TMP[[iBasin]] <- BasinData;
  }


    ##__EXPORT_________________________________________________________________________________####
      for(iBasin in 1:length(ListBasinCode)){
      BasinData <- LIST_TMP[[iBasin]];

      ##format_RData
      FileExport_R <- paste(DIR_EXPORT_RDA,BasinData$BasinCode,".rda",sep="");
      BasinInfo    <- list(BasinCode=BasinData$BasinCode,BasinName=BasinData$BasinName,BasinArea=BasinData$BasinArea_km2,HypsoCurve=BasinData$HypsoData);
      BasinObs     <- data.frame(DatesR=BasinData$TabDatesR,
                                 P=BasinData$TabObsP,T=BasinData$TabObsT,E=BasinData$TabObsE,
                                 Qls=BasinData$TabObsQls,Qmm=BasinData$TabObsQmm)   
      save(BasinInfo,BasinObs,file=FileExport_R,compress="xz")



      ##format_CSV
      FileExport_T <- paste(DIR_EXPORT_TXT,BasinData$BasinCode,".txt"  ,sep="");
      TabObsQls_TMP <- BasinData$TabObsQls; TabObsQls_TMP[is.na(TabObsQls_TMP)] <- (-9); TabObsQls_TMP <- round(TabObsQls_TMP); 
      TabObsQmm_TMP <- BasinData$TabObsQmm; TabObsQmm_TMP[is.na(TabObsQmm_TMP)] <- (-9);
      MYSERIES <- data.frame(BasinData$TabDatesT,formatC(TabObsQls_TMP,format="d",width=8),formatC(cbind(TabObsQmm_TMP,BasinData$TabObsP,BasinData$TabObsT,BasinData$TabObsE),format="f",width=6,digits=2))
      write.table(file=FileExport_T,"# -------------------------------------------------------",append=FALSE,col.names=FALSE,row.names=FALSE,quote=FALSE)
      write.table(file=FileExport_T,paste("# Code       ; ",BasinData$BasinCode   ,sep="")     ,append=TRUE ,col.names=FALSE,row.names=FALSE,quote=FALSE)
      write.table(file=FileExport_T,paste("# Name       ; ",BasinData$BasinName   ,sep="")     ,append=TRUE ,col.names=FALSE,row.names=FALSE,quote=FALSE)
      write.table(file=FileExport_T,paste("# Area (km2) ; ",BasinData$BasinArea_km2,sep="")    ,append=TRUE ,col.names=FALSE,row.names=FALSE,quote=FALSE)
      write.table(file=FileExport_T,"# Units      ; YYYYMMDD ; l/s ; mm/d ; mm/d ; degC ; mm/d",append=TRUE ,col.names=FALSE,row.names=FALSE,quote=FALSE)
      write.table(file=FileExport_T,"# -------------------------------------------------------",append=TRUE ,col.names=FALSE,row.names=FALSE,quote=FALSE)
      write.table(file=FileExport_T,"    Date;       Q;   Qmm;  Ptot;  Temp;    PE"            ,append=TRUE ,col.names=FALSE,row.names=FALSE,quote=FALSE)
      write.table(file=FileExport_T,MYSERIES                                      ,sep=";"     ,append=TRUE ,col.names=FALSE,row.names=FALSE,quote=FALSE)
      }
      FileExport_T <- paste(DIR_EXPORT_TXT,"List_HypsoData",".txt"  ,sep="");
      write.table(file=FileExport_T,paste("    CODE Zmin",paste("  Z",formatC(1:99,width=2,flag="0"),collapse="",sep="")," Zmax",sep="")              ,append=FALSE,col.names=FALSE,row.names=FALSE,quote=FALSE)
      for(iBasin in 1:length(ListBasinCode)){
      write.table(file=FileExport_T,paste(LIST_TMP[[iBasin]]$BasinCode,paste(formatC(LIST_TMP[[iBasin]]$HypsoData,format="f",width=5,digits=0),collapse="",sep=""),sep=""),append=TRUE ,col.names=FALSE,row.names=FALSE,quote=FALSE); }