Revision 4d20f64c
Added by Alberto Guzman over 10 years ago
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
Minor changes to output messages and updated a few script paths after moving to larger disk.