Revision da01a629
Added by Benoit Parmentier over 11 years ago
| climate/research/oregon/interpolation/master_script_temp.R | ||
|---|---|---|
| 10 | 10 |
#STAGE 5: Output analyses: assessment of results for specific dates... |
| 11 | 11 |
# |
| 12 | 12 |
#AUTHOR: Benoit Parmentier |
| 13 |
#DATE: 06/05/2013
|
|
| 13 |
#DATE: 06/07/2013
|
|
| 14 | 14 |
|
| 15 | 15 |
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#363, TASK$568-- |
| 16 | 16 |
|
| ... | ... | |
| 56 | 56 |
clim_script <- file.path(script_path,"climatology_05142013.py") # LST climatology python script |
| 57 | 57 |
grass_setting_script <- file.path(script_path,"grass-setup.R") #Set up system shell environment for python+GRASS |
| 58 | 58 |
source(file.path(script_path,"download_and_produce_MODIS_LST_climatology_05302013.R")) |
| 59 |
source(file.path(script_path,"covariates_production_temperatures_05302013.R"))
|
|
| 59 |
source(file.path(script_path,"covariates_production_temperatures_06072013.R"))
|
|
| 60 | 60 |
source(file.path(script_path,"Database_stations_covariates_processing_function_05212013.R")) |
| 61 | 61 |
source(file.path(script_path,"GAM_fusion_analysis_raster_prediction_multisampling_06052013.R")) |
| 62 | 62 |
source(file.path(script_path,"results_interpolation_date_output_analyses_05062013.R")) |
| ... | ... | |
| 73 | 73 |
stages_to_run<-c(0,2,3,4,5) #May decide on antoher strategy later on... |
| 74 | 74 |
|
| 75 | 75 |
var<-"TMAX" # variable being interpolated |
| 76 |
out_prefix<-"_365d_kriging_day_lst_06052013" #User defined output prefix |
|
| 77 |
out_suffix<-"_OR_06052013" |
|
| 78 |
out_suffix_modis <-"_05302013" #use tiles produce previously |
|
| 76 |
#out_prefix<-"_365d_kriging_day_lst_06072013" #User defined output prefix |
|
| 77 |
out_prefix<-"_365d_GAM_CAI_all_lst_06072013" #User defined output prefix |
|
| 78 |
out_suffix<-"_VE_06072013" |
|
| 79 |
#out_suffix_modis <-"_05302013" #use tiles produce previously |
|
| 80 |
out_suffix_modis <-"_05242013" #use tiles produce previously |
|
| 79 | 81 |
|
| 80 | 82 |
#interpolation_method<-c("gam_fusion","gam_CAI","gam_daily") #other otpions to be added later
|
| 81 |
#interpolation_method<-c("gam_CAI") #other otpions to be added later
|
|
| 83 |
interpolation_method<-c("gam_CAI") #other otpions to be added later
|
|
| 82 | 84 |
#interpolation_method<-c("gam_fusion") #other otpions to be added later
|
| 83 | 85 |
#interpolation_method<-c("gam_daily") #other otpions to be added later
|
| 84 |
interpolation_method<-c("kriging_daily") #other otpions to be added later
|
|
| 86 |
#interpolation_method<-c("kriging_daily") #other otpions to be added later
|
|
| 85 | 87 |
|
| 86 |
#out_path <- paste("/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/output_data",
|
|
| 87 |
# out_prefix,"/",sep="") |
|
| 88 |
out_path<-"/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data" |
|
| 88 |
out_path <-"/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/output_data" |
|
| 89 |
#out_path<-"/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data" |
|
| 89 | 90 |
out_path <-paste(out_path,out_prefix,sep="") |
| 90 | 91 |
|
| 91 | 92 |
if (!file.exists(out_path)){
|
| ... | ... | |
| 103 | 104 |
#infile_reg_outline<- "/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/outline_venezuela_region__VE_01292013.shp" |
| 104 | 105 |
#infile_covariates<-"/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/covariates__venezuela_region_TMIN__VE_03192013.tif" #covariates stack for TMIN |
| 105 | 106 |
#infile_covariates<- "/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/covariates_Oregon_region_TMAX__OR_04052013.tif" #Oregon covar TMAX from earlier codes...for continuity |
| 106 |
#infile_reg_outline="" #input region outline defined by polygon: none for Venezuela
|
|
| 107 |
infile_reg_outline="" #input region outline defined by polygon: none for Venezuela |
|
| 107 | 108 |
#This is the shape file of outline of the study area #It is an input/output of the covariate script |
| 108 |
infile_reg_outline <- "/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/OR83M_state_outline.shp" #input region outline defined by polygon: Oregon |
|
| 109 |
#infile_reg_outline <- "/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/OR83M_state_outline.shp" #input region outline defined by polygon: Oregon
|
|
| 109 | 110 |
#infile_reg_outline <-"OR83M_state_outline.shp" #remove this parameter!!! |
| 110 |
#ref_rast_name<-"" #local raster name defining resolution, exent, local projection--. set on the fly??
|
|
| 111 |
ref_rast_name<-"" #local raster name defining resolution, exent, local projection--. set on the fly?? |
|
| 111 | 112 |
#this may be redundant with infile_reg_outline |
| 112 |
ref_rast_name<-"/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/mean_day244_rescaled.rst" #local raster name defining resolution, exent: oregon |
|
| 113 |
#ref_rast_name<-"/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/mean_day244_rescaled.rst" #local raster name defining resolution, exent: oregon
|
|
| 113 | 114 |
|
| 114 | 115 |
#covar_names see stage 2 |
| 115 | 116 |
|
| 116 |
#list_tiles_modis <- c("h11v08,h11v07,h12v07,h12v08,h10v07,h10v08") #tile for Venezuela and surrounding area
|
|
| 117 |
list_tiles_modis <- c("h08v04,h09v04") #tiles for Oregon
|
|
| 117 |
list_tiles_modis <- c("h11v08,h11v07,h12v07,h12v08,h10v07,h10v08") #tile for Venezuela and surrounding area
|
|
| 118 |
#list_tiles_modis <- c("h08v04,h09v04") #tiles for Oregon
|
|
| 118 | 119 |
|
| 119 |
#CRS_interp<-"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs";
|
|
| 120 |
CRS_interp <-"+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 +ellps=GRS80 +units=m +no_defs"; |
|
| 120 |
CRS_interp<-"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"; |
|
| 121 |
#CRS_interp <-"+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 +ellps=GRS80 +units=m +no_defs";
|
|
| 121 | 122 |
#"+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80" |
| 122 | 123 |
CRS_locs_WGS84<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84
|
| 123 |
#out_region_name<-"_venezuela_region" #generated on the fly
|
|
| 124 |
out_region_name<-"_oregon_region" #generated on the fly |
|
| 124 |
out_region_name<-"_venezuela_region" #generated on the fly |
|
| 125 |
#out_region_name<-"_oregon_region" #generated on the fly
|
|
| 125 | 126 |
|
| 126 | 127 |
#The names of covariates can be changed...these names should be output/input from covar script!!! |
| 127 | 128 |
rnames<-c("x","y","lon","lat","N","E","N_w","E_w","elev_s","slope","aspect","CANHEIGHT","DISTOC")
|
| ... | ... | |
| 186 | 187 |
range_years<-c("2010","2011") #right bound not included in the range!!
|
| 187 | 188 |
range_years_clim<-c("2000","2011") #right bound not included in the range!!
|
| 188 | 189 |
infile_ghncd_data <-"/home/layers/data/climate/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??
|
|
| 190 |
qc_flags_stations<-c("0","S") #flags allowed for screening after the query from the GHCND??
|
|
| 191 |
#qc_flags_stations<-c("0") #flags allowed for screening after the query from the GHCND??
|
|
| 191 | 192 |
|
| 192 | 193 |
#infile_covariates and infile_reg_outline defined in stage 2 or at the start of script... |
| 193 | 194 |
|
| ... | ... | |
| 229 | 230 |
|
| 230 | 231 |
#Models to run...this can be change for each run |
| 231 | 232 |
|
| 232 |
#list_models<-c("y_var ~ s(elev_s)",
|
|
| 233 |
# "y_var ~ s(LST)",
|
|
| 234 |
# "y_var ~ s(elev_s,LST)",
|
|
| 235 |
# "y_var ~ s(lat) + s(lon)+ s(elev_s)",
|
|
| 236 |
# "y_var ~ s(lat,lon,elev_s)",
|
|
| 237 |
# "y_var ~ s(lat,lon) + s(elev_s) + s(N_w,E_w) + s(LST)",
|
|
| 238 |
# "y_var ~ s(lat,lon) + s(elev_s) + s(N_w,E_w) + s(LST) + s(LC2)",
|
|
| 239 |
# "y_var ~ s(lat,lon) + s(elev_s) + s(N_w,E_w) + s(LST) + s(LC6)",
|
|
| 240 |
# "y_var ~ s(lat,lon) + s(elev_s) + s(N_w,E_w) + s(LST) + s(DISTOC)")
|
|
| 233 |
list_models<-c("y_var ~ s(elev_s)",
|
|
| 234 |
"y_var ~ s(LST)", |
|
| 235 |
"y_var ~ s(elev_s,LST)", |
|
| 236 |
"y_var ~ s(lat) + s(lon)+ s(elev_s)", |
|
| 237 |
"y_var ~ s(lat,lon,elev_s)", |
|
| 238 |
"y_var ~ s(lat,lon) + s(elev_s) + s(N_w,E_w) + s(LST)", |
|
| 239 |
"y_var ~ s(lat,lon) + s(elev_s) + s(N_w,E_w) + s(LST) + s(LC2)", |
|
| 240 |
"y_var ~ s(lat,lon) + s(elev_s) + s(N_w,E_w) + s(LST) + s(LC6)", |
|
| 241 |
"y_var ~ s(lat,lon) + s(elev_s) + s(N_w,E_w) + s(LST) + s(DISTOC)") |
|
| 241 | 242 |
#krmod2<-autoKrige(tmax~x_OR83M+y_OR83M,input_data=data_s,new_data=s_sgdf,data_variogram=data_s) |
| 242 |
list_models<-c("y_var ~ 1",
|
|
| 243 |
"y_var ~ x + y", |
|
| 244 |
"y_var ~ x + y + elev_s", |
|
| 245 |
"y_var ~ x + y + DISTOC", |
|
| 246 |
"y_var ~ x + y + elev_s + DISTOC", |
|
| 247 |
"y_var ~ x + y + N_w + E_w", |
|
| 248 |
"y_var ~ LST", |
|
| 249 |
"y_var ~ x + y + LST", |
|
| 250 |
"y_var ~ x + y + elev_s + LST") |
|
| 243 |
#list_models<-c("y_var ~ 1",
|
|
| 244 |
# "y_var ~ x + y",
|
|
| 245 |
# "y_var ~ x + y + elev_s",
|
|
| 246 |
# "y_var ~ x + y + DISTOC",
|
|
| 247 |
# "y_var ~ x + y + elev_s + DISTOC",
|
|
| 248 |
# "y_var ~ x + y + N_w + E_w",
|
|
| 249 |
# "y_var ~ LST",
|
|
| 250 |
# "y_var ~ x + y + LST",
|
|
| 251 |
# "y_var ~ x + y + elev_s + LST")
|
|
| 251 | 252 |
|
| 252 | 253 |
#Default name of LST avg to be matched |
| 253 | 254 |
lst_avg<-c("mm_01","mm_02","mm_03","mm_04","mm_05","mm_06","mm_07","mm_08","mm_09","mm_10","mm_11","mm_12")
|
Also available in: Unified diff
running GAM CAI with monthly stations averages for 2000-2010 time period