Project

General

Profile

« Previous | Next » 

Revision 4d20f64c

Added by Alberto Guzman over 10 years ago

Minor changes to output messages and updated a few script paths after moving to larger disk.

View differences:

climate/research/world/interpolation/Database_stations_covariates_processing_function_01062014.R
114 114

  
115 115
  drv <- dbDriver("PostgreSQL")
116 116
  db <- dbConnect(drv, dbname=db.name)
117
 
117

  
118 118
  time1<-proc.time()    #Start stop watch
119 119
  list_s<-format_s(stat_reg$STAT_ID)
120 120
  data2<-dbGetQuery(db, paste("SELECT *
......
140 140
  data_d <-data_reg  #data_d: daily data containing the query without screening
141 141
  #data_reg <-subset(data_d,mflag=="0" | mflag=="S") #should be input arguments!!
142 142
  #Transform the query to be depending on the number of flags
143

  
143
  
144 144
  data_reg <-subset(data_d, mflag %in% qc_flags_stations) #screening using flags
145 145
  #data_reg2 <-subset(data_d,mflag==qc_flags_stations[1] | mflag==qc_flags_stations[2]) #screening using flags
146 146

  
......
171 171
  ###################################################################
172 172
  ### STEP 4: Extract values at stations from covariates stack of raster images
173 173
  #Eventually this step may be skipped if the covariates information is stored in the database...
174

  
175 174
  #s_raster<-stack(file.path(in_path,infile_covariates))                   #read in the data stack
176 175
  s_raster<-brick(infile_covariates)                   #read in the data stack
177 176
  names(s_raster)<-covar_names               #Assigning names to the raster layers: making sure it is included in the extraction
......
181 180
  #create a shape file and data_frame with names
182 181
  data_RST<-as.data.frame(stat_val)                                            #This creates a data frame with the values extracted
183 182
  data_RST_SDF<-cbind(data_reg,data_RST)
183

  
184 184
  coordinates(data_RST_SDF)<-coordinates(data_reg) #Transforming data_RST_SDF into a spatial point dataframe
185 185
  CRS_reg<-proj4string(data_reg)
186 186
  proj4string(data_RST_SDF)<-CRS_reg  #Need to assign coordinates...
......
216 216
  #year_start_clim: set at the start of the script
217 217
  time1<-proc.time()    #Start stop watch
218 218
  list_s<-format_s(stat_reg$STAT_ID)
219
  print(proc.time()-time1)
219 220
  data_m<-dbGetQuery(db, paste("SELECT *
220 221
                               FROM ghcn
221 222
                               WHERE element=",shQuote(var),
......
246 247
  #In Venezuela and other regions where there are not many stations...mflag==S should be added..see Durenne etal.2010.
247 248

  
248 249
  #d<-subset(data_m,mflag==qc_flags_stations[1] | mflag==qc_flags_stations[2])
249
  d<-subset(data_m,mflag %in% qc_flags_stations)
250
  d<-subset(data_m,mflag %in% qc_flags_stations) 
250 251
  
251 252
  #Add screening here ...May need some screeing??? i.e. range of temp and elevation...
252 253
  
climate/research/world/interpolation/GAM_fusion_analysis_raster_prediction_multisampling_01062014.R
178 178
  infile_daily_rds <-sub(".shp",".rds",infile_daily)
179 179
  if (file.exists(infile_daily_rds) == TRUE){
180 180
     ghcn<-readRDS(infile_daily_rds)
181
     CRS_interp<-proj4string(ghcn)
181 182
  }else{
182 183
    ghcn<-readOGR(dsn=dirname(infile_daily),layer=sub(".shp","",basename(infile_daily)))
183 184
    CRS_interp<-proj4string(ghcn)                       #Storing projection information (ellipsoid, datum,etc.)
......
186 187

  
187 188
  infile_locs_rds<-sub(".shp",".rds",infile_locs)
188 189
  if (file.exists(infile_locs_rds) == TRUE){
189
     readRDS(infile_locs_rds) 
190
    stat_loc<-readRDS(infile_locs_rds) 
190 191
  }else{
191 192
    stat_loc<-readOGR(dsn=dirname(infile_locs),layer=sub(".shp","",basename(infile_locs)))
192 193
    saveRDS(stat_loc,infile_locs_rds)
......
278 279
  }
279 280
  
280 281
  if (interpolation_method %in% c("gam_CAI","kriging_CAI", "gwr_CAI")){
282
    num_cores2 = as.integer(num_cores) + 2
281 283
    list_param_runClim_KGCAI<-list(j,s_raster,covar_names,lst_avg,list_models,dst,sampling_month_obj,var,y_var_name, out_prefix,out_path)
282 284
    names(list_param_runClim_KGCAI)<-c("list_index","covar_rast","covar_names","lst_avg","list_models","dst","sampling_month_obj","var","y_var_name","out_prefix","out_path")
283
    clim_method_mod_obj<-mclapply(1:length(sampling_month_obj$ghcn_data), list_param=list_param_runClim_KGCAI, runClim_KGCAI,mc.preschedule=FALSE,mc.cores = num_cores) #This is the end bracket from mclapply(...) statement
285
    clim_method_mod_obj<-mclapply(1:length(sampling_month_obj$ghcn_data), list_param=list_param_runClim_KGCAI, runClim_KGCAI,mc.preschedule=FALSE,mc.cores = num_cores2) #This is the end bracket from mclapply(...) statement
284 286
    #test<-runClim_KGCAI(1,list_param=list_param_runClim_KGCAI)
287

  
285 288
    save(clim_method_mod_obj,file= file.path(out_path,paste(interpolation_method,"_mod_",y_var_name,out_prefix,".RData",sep="")))
286 289
    list_tmp<-vector("list",length(clim_method_mod_obj))
287 290
    for (i in 1:length(clim_method_mod_obj)){
......
306 309
  #method_mod_obj<-mclapply(1:1, runGAMFusion,mc.preschedule=FALSE,mc.cores = 1) #This is the end bracket from mclapply(...) statement
307 310

  
308 311
  #Set raster options for daily steps
309
  rasterOptions(timer=TRUE,chunksize=1e+05)
312
  #rasterOptions(timer=TRUE,chunksize=1e+05)
313
  #This is for high station count areas
314
  rasterOptions(timer=TRUE,chunksize=1e+04)
310 315

  
311 316
  cat("Predictions at the daily scale:",file=log_fname,sep="\n", append=TRUE)
312 317
  t1<-proc.time()
......
516 521
    save(raster_prediction_obj,file= file.path(out_path,paste("raster_prediction_obj_",interpolation_method,"_", y_var_name,out_prefix,".RData",sep="")))
517 522
    
518 523
  }
519
  
524

  
520 525
  return(raster_prediction_obj)
521 526
}
522 527

  
climate/research/world/interpolation/GAM_fusion_function_multisampling_01062014.R
27 27
    mod <-in_models[[i]] #accessing GAM model ojbect "j"
28 28
    raster_name<-out_filename[[i]]
29 29
    if (inherits(mod,"gam")) {           #change to c("gam","autoKrige")
30
      raster_pred<- predict(object=r_stack,model=mod,na.rm=FALSE,block.size=1000) #Using the coeff to predict new values.
31
      raster_pred<- predict(object=r_stack,model=mod,na.rm=FALSE) #Using the coeff to predict new values.
30
      #raster_pred<- predict(object=r_stack,model=mod,na.rm=FALSE,block.size=1000) #Using the coeff to predict new values.
31
      #raster_pred<- predict(object=r_stack,model=mod,na.rm=FALSE) #Using the coeff to predict new values.
32
      raster_pred<- predict(object=r_stack,model=mod,filename=raster_name,na.rm=FALSE,overwrite=TRUE) #Using the coeff to predict new values. 
32 33
      names(raster_pred)<-"y_pred"  
33
      writeRaster(raster_pred, filename=raster_name,overwrite=TRUE)  #Writing the data in a raster file format...(IDRISI)
34
      #writeRaster(raster_pred, filename=raster_name)#,overwrite=TRUE)  #Writing the data in a raster file format...(IDRISI)
34 35
      #print(paste("Interpolation:","mod", j ,sep=" "))
35 36
      list_rast_pred[[i]]<-raster_name
36 37
    }
......
45 46
  #This functions several models and returns model objects.
46 47
  #Arguments: - list of formulas for GAM models
47 48
  #           - fitting data in a data.frame or SpatialPointDataFrame
48
  #Output: list of model objects 
49
  #Output: list of model objects
49 50
  list_fitted_models<-vector("list",length(list_formulas))
50 51
  for (k in 1:length(list_formulas)){
51 52
    formula<-list_formulas[[k]]
......
74 75
  list_formulas<-lapply(list_models,as.formula,env=.GlobalEnv) #mulitple arguments passed to lapply!!
75 76
  cname<-paste("mod",1:length(list_formulas),sep="") #change to more meaningful name?
76 77
  
77
  names(list_out_filename)<-cname  
78
  
78
  names(list_out_filename)<-cname   
79 79
  ##Now carry out prediction
80
  if(method_interp=="gam"){    
81
    
80
  if(method_interp=="gam"){ 
82 81
    #First fitting
83 82
    mod_list<-fit_models(list_formulas,data_df) #only gam at this stage
84 83
    names(mod_list)<-cname
......
164 163
  date_month <-strptime(sampling_month_dat$date[j], "%Y%m%d")   # interpolation date being processed
165 164
  month_no <-strftime(date_month, "%m")          # current month of the date being processed
166 165
  LST_month<-paste("mm_",month_no,sep="") # name of LST month to be matched
167
  LST_name <-LST_month  
166
  LST_name <-LST_month
167

  
168
  print(paste("Processing month ",j))  
168 169
  #### STEP 2: PREPARE DATA
169 170
    
170 171
  #change here...use training data...
......
234 235
  ## Select the relevant method...
235 236
  
236 237
  if (interpolation_method=="gam_CAI"){
237
    
238 238
    #First fitting
239 239
    mod_list<-fit_models(list_formulas,data_month) #only gam at this stage
240 240
    names(mod_list)<-cname
......
273 273
  mod_krtmp1<-fitclim
274 274
  model_name<-"mod_kr"
275 275
  
276
  clim_rast<-interpolate(LST,fitclim) #interpolation using function from raster package
276
  #clim_rast<-interpolate(LST,fitclim) #interpolation using function from raster package
277 277
  
278 278
  #Write out modeled layers
279 279
  data_name<-paste(var,"_clim_month_",as.integer(month_no),"_",model_name,"_",prop_month,
280 280
                   "_",run_samp,sep="")
281 281
  raster_name_clim<-file.path(out_path,paste("CAI_",data_name,out_prefix,".tif", sep=""))
282

  
283
  clim_rast<-interpolate(LST,fitclim) #interpolation using function from raster package
282 284
  writeRaster(clim_rast, filename=raster_name_clim,overwrite=TRUE)  #Writing the data in a raster file format...(IDRISI)
283
  
285
    
284 286
  #Adding to current objects
285 287
  mod_list[[model_name]]<-mod_krtmp1
286 288
  #rast_bias_list[[model_name]]<-raster_name_bias
......
749 751
      rast_clim_mod <- stack(rast_clim_list)
750 752
      names(rast_clim_mod) <- names(rast_clim_list)
751 753
      rast_clim_month <- subset(rast_clim_mod,1) #example layer to interpolate to
752
      
754
           
753 755
      daily_delta_rast<-interpolate(rast_clim_month,fitdelta) #Interpolation of the of the daily devation
754 756
      #there is only one daily devation (delta) sruface in this case
755 757
      
......
758 760
      data_name<-paste("daily_delta_",y_var_name,"_",model_name,"_",sampling_month_dat$prop[index_m],"_",sampling_month_dat$run_samp[index_m],"_",
759 761
                       sampling_dat$date[index_d],"_",sampling_dat$prop[index_d],"_",sampling_dat$run_samp[index_d],sep="")
760 762
      raster_name_delta<-file.path(out_path,paste(interpolation_method,"_",var,"_",data_name,out_prefix,".tif", sep=""))
763

  
764
      
761 765
      writeRaster(daily_delta_rast, filename=raster_name_delta,overwrite=TRUE)  #Writing the data in a raster file format...(IDRISI)
762 766
      
767
     
763 768
      list_daily_delta_rast[[1]] <- raster_name_delta
764 769
      list_mod_krtmp2[[1]] <- mod_krtmp2
765 770
    }
......
804 809
  }
805 810
  
806 811
  if(use_clim_image==TRUE){
807
    
808 812
    # User can choose to join daily and monthly station before interpolation:
809 813
    #-this ensures that the delta difference is "more" exact since its starting point is basesd on average value but there is risk to loose some stations
810 814
    #may need to change this option later!!
......
870 874
        #rast_clim_list<-rast_clim_yearlist[[index_m]]  #select relevant monthly climatology image ...
871 875
        rast_clim_month <- subset(rast_clim_mod,1) #example layer to interpolate to
872 876
        
877
        
873 878
        daily_delta_rast<-interpolate(rast_clim_month,fitdelta) #Interpolation of the bias surface...
874 879
        
875 880
        #list_daily_delta_rast[[k]] <- raster_name_delta
......
877 882
        data_name<-paste("daily_delta_",y_var_name,"_",model_name,"_",sampling_month_dat$prop[index_m],"_",sampling_month_dat$run_samp[index_m],"_",
878 883
                         sampling_dat$date[index_d],"_",sampling_dat$prop[index_d],"_",sampling_dat$run_samp[index_d],sep="")
879 884
        raster_name_delta<-file.path(out_path,paste(interpolation_method,"_",var,"_",data_name,out_prefix,".tif", sep=""))
880
        
885
              
886
        daily_delta_rast<-interpolate(rast_clim_month,fitdelta) #Interpolation of the of the daily devation
887
        #there is only one daily devation (delta) sruface in this case
881 888
        writeRaster(daily_delta_rast, filename=raster_name_delta,overwrite=TRUE)  #Writing the data in a raster file format...(IDRISI)
882
        #writeRaster(r_spat, NAflag=NA_flag_val,filename=raster_name,bylayer=TRUE,bandorder="BSQ",overwrite=TRUE)   
889
          #writeRaster(r_spat, NAflag=NA_flag_val,filename=raster_name,bylayer=TRUE,bandorder="BSQ",overwrite=TRUE)   
883 890
        
884 891
        #raster_name_delta <- list_daily_delta_rast
885 892
        #mod_krtmp2 <- list_mod_krtmp2
climate/research/world/interpolation/climatologyScripts/climatology.pyx
1 1
cimport cython
2 2
import numpy as np
3 3
cimport numpy as np
4
import scipy.interpolate as inter
5
import matplotlib.pyplot as plt
6
from scipy.interpolate import interp1d
7
from scipy.optimize import curve_fit
8
from scipy.optimize import leastsq
9
import sys
10
import pylab
11
from cpython cimport bool
4 12
@cython.boundscheck(False)
5 13
@cython.wraparound(False)
6 14

  
7
def addLST(np.ndarray [np.double_t, ndim=3] meanTotal, np.ndarray [np.double_t, ndim=2] lst, np.ndarray [np.uint8_t, ndim=2] qa,np.ndarray [np.float32_t, ndim=3] nobsTotal):
15

  
16
def addLST(np.ndarray [np.float32_t, ndim=3] meanTotal, np.ndarray [np.double_t, ndim=2] lst, np.ndarray [np.uint8_t, ndim=2] qa,np.ndarray [np.float32_t, ndim=3] nobsTotal):
8 17
   cdef int yDim = 1200
9 18
   cdef int xDim = 1200
10 19
   cdef int zDim = 4 #Debug 5
......
55 64
                       
56 65
   return meanTotal,nobsTotal        
57 66

  
58
def mergeBands(np.ndarray [np.double_t, ndim=3] meanTotal,np.ndarray [np.float32_t, ndim=3] nobsTotal,int nobsLow):
67
def mergeBands(np.ndarray [np.float32_t, ndim=3] meanTotal,np.ndarray [np.float32_t, ndim=3] nobsTotal,int nobsLow):
59 68
   cdef int yDim = 1200
60 69
   cdef int xDim = 1200
61 70
   cdef int zDim = 4
......
65 74
   cdef double lstVal = 0
66 75
   cdef double nobsVal = 0
67 76
   cdef int qaVal = 0
77
   cdef double nodata = -9999.0
68 78
     
69
   cdef np.ndarray[np.double_t, ndim=2] outLST = np.zeros([yDim,xDim], dtype=np.double) 
79
   cdef np.ndarray[np.float32_t, ndim=2] outLST = np.ones([yDim,xDim], dtype=np.float32) * -9999.0
70 80
   cdef np.ndarray[np.float32_t, ndim=2] outNobs = np.zeros([yDim,xDim], dtype=np.float32)
71 81
   cdef np.ndarray[np.int8_t, ndim=2] outQA = np.zeros([yDim,xDim], dtype=np.int8)
72 82

  
83
   #Temp
84
   nobsCheck = [20,10,5]
85

  
73 86
   for y in range(0,yDim):
74 87
     for x in range(0,xDim):
75
       
76 88
       #From high qa(1) to low(4)
77 89
       for z in range(0,4):
78
         lstVal = meanTotal[<unsigned int>z,<unsigned int>y,<unsigned int>x]
79 90
         nobsVal = nobsTotal[<unsigned int>z,<unsigned int>y,<unsigned int>x]
80 91

  
81 92
         #It's a climatology so we should have at least nobs > nobsLow
82 93
         if nobsVal > nobsLow:
83
           outLST[<unsigned int>y,<unsigned int>x] = lstVal
94
           lstVal = meanTotal[<unsigned int>z,<unsigned int>y,<unsigned int>x]
95
           outLST[<unsigned int>y,<unsigned int>x] = lstVal/nobsVal
84 96
           outNobs[<unsigned int>y,<unsigned int>x] = nobsVal
85 97
           outQA[<unsigned int>y,<unsigned int>x] = z+1
86 98
           break
99
        
100
       #Fill any that were not filled with 20 using 10 or 5?
101
       lstVal = outLST[<unsigned int>y,<unsigned int>x]
102
       if lstVal < -999:
103
         for z in range(0,4):
104
           nobsVal = nobsTotal[<unsigned int>z,<unsigned int>y,<unsigned int>x]     
105

  
106
           if nobsVal > 5:
107
             lstVal = meanTotal[<unsigned int>z,<unsigned int>y,<unsigned int>x]
87 108

  
109
             outLST[<unsigned int>y,<unsigned int>x] = lstVal/nobsVal
110
             outNobs[<unsigned int>y,<unsigned int>x] = nobsVal
111
             outQA[<unsigned int>y,<unsigned int>x] = z+1+5 #DEBUG adding 5 
112
             break 
113
      
88 114
   return outLST,outNobs,outQA
89 115
         
116
def temporalFit(np.ndarray [np.float32_t, ndim=3] lstAll,np.ndarray [np.float32_t, ndim=3] nobsAll,
117
                np.ndarray [np.int8_t, ndim=3] qaAll,fitMet="spline", int nobsLow=20,tension=3):
118
   cdef int yDim = 1200
119
   cdef int xDim = 1200
120
   cdef int zDim = 5
121
   cdef int mDim = 12
122
   cdef int y = 0
123
   cdef int x = 0
124
   cdef int z = 0
125
   cdef int m = 0
126
   cdef int i = 0
127
   cdef double lstVal = 0
128
   
129
   cdef double nobsVal = 0
130
   cdef int qaVal = 0
131
   cdef double nodata = -9999.0
132
   cdef int qaFillVal = 11 
133
   cdef int count = 0 
134
   cdef int filled = 0 
135
  
136
   cdef np.ndarray[np.float32_t, ndim=1] ySeries = np.ones([mDim], dtype=np.float32) * -9999.0
137
   cdef np.ndarray[np.int32_t, ndim=1] xSeries = np.zeros([mDim], dtype=np.int32)
138

  
139
   for x in range(0,12):
140
     xSeries[<unsigned int>x] = x
141

  
142
   for y in range(0,yDim):
143
     for x in range(0,xDim):
144
       #Get the time series for pixel
145
       for m in range(0,mDim):
146
         ySeries[<unsigned int>m] = lstAll[<unsigned int>m,<unsigned int>y,<unsigned int>x]
147
         
148
       #Get rid of no data
149
       ySub = ySeries[ySeries>-999]
150
       xSub = xSeries[ySeries>-999]
151

  
152
       #Means water pixel, there's not enough to fit
153
       #or fulltimeseries       
154
       if (len(ySub) < 7) or (len(ySub) == 12):
155
          continue
156

  
157
       count += 1
158

  
159
       if fitMet == "spline":
160
         #fitY = inter.InterpolatedUnivariateSpline(xSub, ySub,k=5)
161
         #fitY = fitY(xSeries)
162
         spline = inter.UnivariateSpline(xSub, ySub,k=tension,s=0.5)
163
         fitY = spline(xSeries)
164
  
165
       for i in range(0,12):
166
         lstVal = lstAll[<unsigned int>i,<unsigned int>y,<unsigned int>x]
167
         if lstVal < -999:
168
           filled += 1
169
           lstAll[<unsigned int>i,<unsigned int>y,<unsigned int>x] = fitY[<unsigned int>i]
170
           qaAll[<unsigned int>i,<unsigned int>y,<unsigned int>x] = 6 #6 is fit from spline
171
           #What to do about nobs? 
172
           nobsAll[<unsigned int>i,<unsigned int>y,<unsigned int>x] = nobsLow 
173

  
174
   print count,filled
175
   return lstAll,nobsAll,qaAll
176

  
177
def getAllNobs(np.ndarray [np.float32_t, ndim=4] nobsTotal):
178
  cdef int yDim = 1200
179
  cdef int xDim = 1200
180
  cdef int zDim = 5
181
  cdef int mDim = 12
182
  cdef int y = 0
183
  cdef int x = 0
184
  cdef int z = 0
185
  cdef double nobsVal0 = 0
186
  cdef double nobsVal1 = 0
187
  cdef double nobsVal2 = 0
188
  cdef double nobsVal = 0
189

  
190
  cdef np.ndarray[np.float32_t,ndim=3] allNobs = np.zeros([mDim,yDim,xDim], dtype=np.float32)
191

  
192
  for m in range(0,12):
193
    for y in range(0,yDim):
194
      for x in range(0,xDim):
195
        nobsVal = 0
196
        #for z in range(0,4):
197
          #nobsVal0 = nobsTotal[<unsigned int>m[0]-1,<unsigned int>z,<unsigned int>y,<unsigned int>x]
198
          #nobsVal1 = nobsTotal[<unsigned int>m[1]-1,<unsigned int>z,<unsigned int>y,<unsigned int>x]
199
          #nobsVal2 = nobsTotal[<unsigned int>m[2]-1,<unsigned int>z,<unsigned int>y,<unsigned int>x]
200

  
201
        #nobsVal = nobsVal + nobsTotal[<unsigned int>m,<unsigned int>z,<unsigned int>y,<unsigned int>x]
202
        allNobs[<unsigned int>m,<unsigned int>y,<unsigned int>x] = nobsTotal[<unsigned int>m,3,<unsigned int>y,<unsigned int>x]
203
        #allNobs[<unsignedint>m,<unsigned int>y,<unsigned int>x] = nobsVal
204

  
205
  return allNobs
206

  
207
def getTS(np.ndarray [np.float32_t, ndim=3] meanTotal,outDir,tile):
208
   cdef int yDim = 1200
209
   cdef int xDim = 1200
210
   cdef int mDim = 12
211
   cdef int y = 0
212
   cdef int x = 0
213
   cdef int z = 0
214

  
215
   cdef double lstVal = 0
216

  
217
   cdef np.ndarray[np.float32_t, ndim=1] ySeries = np.ones([mDim], dtype=np.float32) * -9999.0
218
   cdef np.ndarray[np.int32_t, ndim=1] xSeries = np.zeros([mDim], dtype=np.int32)
219

  
220
   fPtr = open(outDir+"timeSeriesOut_"+tile+".txt","w+")
221
   if fPtr is None:
222
      print "Error opening in file"
223
      return None
224

  
225
   for x in range(0,12):
226
     xSeries[<unsigned int>x] = x
227

  
228
   for y in range(0,yDim):
229
     for x in range(0,xDim):
230
       #Get the time series for pixel
231
       for m in range(0,mDim):
232
         ySeries[<unsigned int>m] = meanTotal[<unsigned int>m,<unsigned int>y,<unsigned int>x]
233
         
234
       #Get rid of no data
235
       ySub = ySeries[ySeries>-999]
236
       xSub = xSeries[ySeries>-999]
237

  
238
       #Most likely means water pixel or there's not enough to fit      
239
       if (len(ySub) < 5) or (len(ySub) == 12):
240
          continue
241

  
242
       stArr = str(map(str, ySeries))
243
       stArr = stArr[1:-1].replace("'","")
244
       fPtr.write(stArr)
245
       fPtr.write('\n')
90 246

  
247
   fPtr.close()
248
   return True
91 249

  
climate/research/world/interpolation/climatologyScripts/mosaicByMonth.py
84 84
           #Now mosaic all files in outFnTxt
85 85
           if l == 'mean': 
86 86
              outType = '-ot Float32'
87
              nodata = "-9999"
87 88
           elif l == 'nobs': 
88 89
              outType = '-ot Int16'
90
              nodata = "0"
89 91
           elif l == 'qa':
90
              outType = '-ot Byte'            
92
              outType = '-ot Byte'
93
              nodata = "0"            
91 94

  
92
           inputStr = ['gdal_merge.py','--config','GDAL_CACHEMAX=1500',outType,'-n',"0","-a_nodata","0","-o",outFn,"--optfile",outFnTxt]
95
           inputStr = ['gdal_merge.py','--config','GDAL_CACHEMAX=1500',outType,'-n',nodata,"-a_nodata","0","-o",outFn,"--optfile",outFnTxt]
93 96
           #inputStr = ['gdalwarp','-srcnodata',"0","-dstnodata","0","-wm","500","-overwrite","--optfile",outFnTxt,outFn]
94 97
           
95 98
           #out = 0
climate/research/world/interpolation/covariates_production_temperature_01062014.R
404 404
    lst_pat<-"LST_Day_1km" #for the time being change at later stage...
405 405
  }
406 406

  
407

  
407 408
  if (process_LST == FALSE){
408 409
    mean_m_list<-rep(NA, 11)
409 410
    nobs_m_list<-rep(NA, 11)
......
569 570
  xy <-coordinates(r1)  #get x and y projected coordinates...
570 571
  CRS_interp<-proj4string(r1)
571 572

  
572
  if (identical(CRS_interp,CRS_locs_WGS84)){  
573
  if (process_LST == FALSE){  
573 574
    xy_latlon<-xy
574 575
  }else{
575 576
    xy_latlon<-project(xy, CRS_interp, inv=TRUE) # find lat long for projected coordinats (or pixels...)
......
602 603
  
603 604
  #Screen values given valid value ranges for specified variables
604 605
  s_raster<-screening_val_r_stack_fun(list_val_range,s_raster)
605
  
606

  
606 607
  #Write out stack of number of change 
607 608
  data_name<-paste("covariates_",out_region_name,"_",sep="")
608 609
  raster_name<-paste(data_name,var,"_",out_suffix,".tif", sep="")
climate/research/world/interpolation/interpolation_method_day_function_multisampling_01062014.R
29 29
    raster_name<-out_filename[[i]]
30 30
    if (inherits(mod,"gam")) {           #change to c("gam","autoKrige")
31 31
      #raster_pred<- predict(object=r_stack,model=mod,na.rm=FALSE) #Using the coeff to predict new values. 
32
      #raster_pred<-predict(object=r_stack,filename=raster_name,model=mod,na.rm=FALSE,overwrite=TRUE) #Using the coeff to predict new values.
33
      raster_pred<-predict(object=r_stack,model=mod,na.rm=FALSE,filename=raster_name,overwrite=TRUE) 
34
      names(raster_pred)<-"y_pred"        
35

  
32
      raster_pred<-predict(object=r_stack,filename=raster_name,model=mod,na.rm=FALSE,overwrite=TRUE) #Using the coeff to predict new values. 
33
        #raster_pred<-predict(object=r_stack,model=mod,na.rm=FALSE,filename=raster_name)#,overwrite=TRUE) 
34
      names(raster_pred)<-"y_pred"
35
      
36
  
36 37
      #Doing this in one step above
37 38
      #writeRaster(raster_pred, filename=raster_name,overwrite=TRUE)  #Writing the data in a raster file format...(IDRISI)
38 39
      #print(paste("Interpolation:","mod", j ,sep=" "))
......
232 233
  list_fitted_models<-vector("list",length(list_formulas))
233 234
  for (k in 1:length(list_formulas)){
234 235
    formula<-list_formulas[[k]]
235
    mod<- try(gam(formula, data=data_training)) #change to any model!!
236
    mod<- try(gam(formula, data=data_training)) #change to any model!! 
236 237
    #mod<- try(autoKrige(formula, input_data=data_s,new_data=s_sgdf,data_variogram=data_s))
237 238
    model_name<-paste("mod",k,sep="")
238 239
    assign(model_name,mod) 
climate/research/world/interpolation/master_script_03012014.R
44 44

  
45 45
#Adding command line arguments to use mpiexec
46 46
args<-commandArgs(TRUE)
47
script_path<-"/nobackup/aguzman4/climateLayers/climateCode/environmental-layers/climate/research/oregon/interpolation"
48
dataHome<-"/nobackup/aguzman4/climateLayers/interp/testdata/"
49
script_path2<-"/nobackup/aguzman4/climateLayers/climateCode/environmental-layers/climate/research/world/interpolation"
47
script_path<-"/nobackupp4/aguzman4/climateLayers/climateCode/environmental-layers/climate/research/oregon/interpolation"
48
dataHome<-"/nobackupp4/aguzman4/climateLayers/interp/testdata/"
49
script_path2<-"/nobackupp4/aguzman4/climateLayers/climateCode/environmental-layers/climate/research/world/interpolation"
50 50

  
51 51
#CALLED FROM MASTER SCRIPT:
52 52

  
......
68 68
source(file.path(script_path,"GAM_fusion_function_multisampling_validation_metrics.R"))
69 69

  
70 70
#stages_to_run<-c(0,2,3,4,5) #MRun only raster fitting, prediction and assessemnt (providing lst averages, covar brick and met stations)
71
stages_to_run<-c(0,2,3,4,0)
71
stages_to_run<-c(0,2,0,0,0)
72
#stages_to_run<-c(0,0,3,0,0)
72 73

  
73 74
#If stage 2 is skipped then use previous covar object
74 75
covar_obj_file<-args[8]
......
81 82
out_suffix<-args[2]                             #Regional suffix
82 83
out_suffix_modis<-args[3]                       #pattern to find tiles produced previously     
83 84

  
84
interpolation_method<-c("gam_fusion") #other otpions to be added later
85
#interpolation_method<-c("gam_fusion") #other otpions to be added later
86
interpolation_method<-c("gam_CAI")
85 87

  
86 88
out_path <- args[4] #"/nobackup/aguzman4/climateLayers/output/"
87 89

  
......
162 164

  
163 165
names(list_param_covar_production)<-c("var","out_path","lc_path","infile_modis_grid","infile_elev","infile_canheight",
164 166
                                      "infile_distoc","list_tiles_modis","infile_reg_outline","CRS_interp","CRS_locs_WGS84","out_region_name",
165
                                      "buffer_dist","list_val_range","out_suffix","out_suffix_modis","ref_rast_name","hdfdir","covar_names","process_LST") 
167
                                   "buffer_dist","list_val_range","out_suffix","out_suffix_modis","ref_rast_name","hdfdir","covar_names","process_LST") 
166 168

  
167 169
## Modify to store infile_covar_brick in output folder!!!
168 170
if (stages_to_run[2]==2){
......
186 188
range_years<-c("2010","2011") #right bound not included in the range!!
187 189
range_years_clim<-c("2000","2011") #right bound not included in the range!!
188 190
infile_ghncd_data <-"/nobackupp4/aguzman4/climateLayers/interp/testdata/data_workflow/inputs/ghcn/v2.92-upd-2012052822/ghcnd-stations.txt"                              #This is the textfile of station locations from GHCND
189
#qc_flags_stations<-c("0","S")    #flags allowed for screening after the query from the GHCND??
190
qc_flags_stations<-c("0")    #flags allowed for screening after the query from the GHCND??
191
qc_flags_stations<-c("0","S")    #flags allowed for screening after the query from the GHCND??
192
#qc_flags_stations<-c("0")    #flags allowed for screening after the query from the GHCND??
191 193

  
192 194
#infile_covariates and infile_reg_outline defined in stage 2 or at the start of script...
193 195

  
......
201 203

  
202 204
if (stages_to_run[3]==3){
203 205
  list_outfiles<-database_covariates_preparation(list_param_prep)
206
}
207
if ((stages_to_run[4]==4) | (stages_to_run[5]==5)){
208
    list_outfiles <-load_obj(met_stations_outfiles_obj_file)
204 209
}else{
205
  if ((stages_to_run[4]==4) | (stages_to_run[5]==5)){
206
    list_outfiles <-load_obj(met_stations_outfiles_obj_file) 
207
  }
208
  else{
209 210
    quit()
210 211
  }
211
}
212 212

  
213 213
############### STAGE 4: RASTER PREDICTION #################
214 214

  

Also available in: Unified diff