Revision 18d48828
Added by Benoit Parmentier over 11 years ago
| climate/research/oregon/interpolation/master_script_temp.R | ||
|---|---|---|
| 3 | 3 |
# |
| 4 | 4 |
##This script produces intperpolated surface of TMIN and TMAX for specified processing region(s) given sets |
| 5 | 5 |
#of inputs and parameters. |
| 6 |
#STAGE 1: LST climatology calculation |
|
| 7 |
#STAGE 2: Covariates preparation: aspect, land cover, distance to coast etc.
|
|
| 8 |
#STAGE 3: Data preparation: database query, extraction of covariates from stack
|
|
| 9 |
#STAGE 4: Raster prediction: run interpolation method -- gam fusion, gam CAI, ...
|
|
| 10 |
#STAGE 5: Output analyses-visualization of results for specific dates...
|
|
| 6 |
#STAGE 1: LST climatology downloading and/or calculation
|
|
| 7 |
#STAGE 2: Covariates preparation for study/processing area: calculation of covariates (spect,land cover,etc.) and reprojection
|
|
| 8 |
#STAGE 3: Data preparation: meteorological station database query and extraction of covariates values from raster brick
|
|
| 9 |
#STAGE 4: Raster prediction: run interpolation method (-- gam fusion, gam CAI, ...) and perform validation
|
|
| 10 |
#STAGE 5: Output analyses: assessment of results for specific dates...
|
|
| 11 | 11 |
# |
| 12 | 12 |
#AUTHOR: Benoit Parmentier |
| 13 |
#DATE: 05/24/2013
|
|
| 13 |
#DATE: 05/27/2013
|
|
| 14 | 14 |
|
| 15 | 15 |
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#363, TASK$568-- |
| 16 | 16 |
|
| ... | ... | |
| 47 | 47 |
|
| 48 | 48 |
##SCRIPT USED FOR THE PREDICTIONS: Source or list all scripts here to avoid confusion on versions being run!!!! |
| 49 | 49 |
|
| 50 |
#source(file.path(script_path,"master_script_temp_05212013.R")) #Master script can be run directly...
|
|
| 50 |
#source(file.path(script_path,"master_script_temp_05272013.R")) #Master script can be run directly...
|
|
| 51 | 51 |
|
| 52 | 52 |
#CALLED FROM MASTER SCRIPT: |
| 53 | 53 |
|
| ... | ... | |
| 55 | 55 |
clim_script <- file.path(script_path,"climatology_05142013.py") # LST climatology python script |
| 56 | 56 |
grass_setting_script <- file.path(script_path,"grass-setup.R") #Set up system shell environment for python+GRASS |
| 57 | 57 |
#source(file.path(script_path,"download_and_produce_MODIS_LST_climatology_05162013.R")) |
| 58 |
source(file.path(script_path,"covariates_production_temperatures_05132013.R"))
|
|
| 58 |
source(file.path(script_path,"covariates_production_temperatures_05272013.R"))
|
|
| 59 | 59 |
source(file.path(script_path,"Database_stations_covariates_processing_function_05212013.R")) |
| 60 | 60 |
source(file.path(script_path,"GAM_fusion_analysis_raster_prediction_multisampling_05212013.R")) |
| 61 | 61 |
source(file.path(script_path,"results_interpolation_date_output_analyses_05062013.R")) |
| ... | ... | |
| 64 | 64 |
#FUNCTIONS CALLED FROM GAM ANALYSIS RASTER PREDICTION ARE FOUND IN... |
| 65 | 65 |
|
| 66 | 66 |
source(file.path(script_path,"sampling_script_functions_03122013.R")) |
| 67 |
source(file.path(script_path,"GAM_fusion_function_multisampling_05062013.R")) #Include GAM_CAI
|
|
| 67 |
source(file.path(script_path,"GAM_fusion_function_multisampling_05212013.R")) #Include GAM_CAI
|
|
| 68 | 68 |
source(file.path(script_path,"GAM_fusion_function_multisampling_validation_metrics_05062013.R")) |
| 69 | 69 |
|
| 70 | 70 |
#stages_to_run<-c(1,2,3,4,5) #May decide on antoher strategy later on... |
| 71 |
stages_to_run<-c(0,0,3,4,5) #May decide on antoher strategy later on...
|
|
| 71 |
stages_to_run<-c(0,2,3,4,5) #May decide on antoher strategy later on...
|
|
| 72 | 72 |
|
| 73 | 73 |
var<-"TMAX" # variable being interpolated |
| 74 | 74 |
out_prefix<-"_365d_GAM_fus_all_lst_05242013" #User defined output prefix |
| 75 |
out_suffix<-"_VE_05242013" |
|
| 76 |
out_suffix_modis <-"_05242013" |
|
| 77 |
|
|
| 75 | 78 |
#interpolation_method<-c("gam_fusion","gam_CAI") #other otpions to be added later
|
| 76 | 79 |
#interpolation_method<-c("gam_CAI") #other otpions to be added later
|
| 77 | 80 |
interpolation_method<-c("gam_fusion") #other otpions to be added later
|
| ... | ... | |
| 88 | 91 |
} |
| 89 | 92 |
|
| 90 | 93 |
lc_path<-"/home/layers/data/land-cover/lc-consensus-global" |
| 91 |
infile_modis_grid<-"/home/layers/commons/modis/modis_sinusoidal/modis_sinusoidal_grid_world.shp" #Give path!!! NEED TO CHANGE THIS... |
|
| 92 |
infile_elev<-"/home/layers/data/terrain/dem-cgiar-srtm-1km-tif/srtm_1km.tif" #this is the global file: replace later with the input produced by the DEM team |
|
| 93 |
infile_canheight<-"/home/layers/data/land-cover/treeheight-simard2011/Simard_Pinto_3DGlobalVeg_JGR.tif" #Canopy height |
|
| 94 |
infile_modis_grid<-"/data/project/layers/commons/modis/modis_sinusoidal/modis_sinusoidal_grid_world.shp" #modis grid tiling system, global |
|
| 95 |
infile_elev<-"/home/layers/data/terrain/dem-cgiar-srtm-1km-tif/srtm_1km.tif" #elevation at 1km, global extent to be replaced by the new fused product |
|
| 96 |
infile_canheight<-"/home/layers/data/land-cover/treeheight-simard2011/Simard_Pinto_3DGlobalVeg_JGR.tif" #Canopy height, global extent |
|
| 97 |
infile_distoc <- "/data/project/layers/commons/distance_to_coast/GMT_intermediate_coast_distance_01d_rev.tif" #distance to coast, global extent at 0.01 deg |
|
| 98 |
#infile_covariates<-"/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/covariates__venezuela_region__VE_01292013.tif" #this is an output from covariate script and used in stage 3 and stage 4 |
|
| 99 |
#infile_reg_outline<- "/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/outline_venezuela_region__VE_01292013.shp" |
|
| 100 |
#infile_covariates<-"/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/covariates__venezuela_region_TMIN__VE_03192013.tif" #covariates stack for TMIN |
|
| 101 |
#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 |
|
| 102 |
#infile_reg_outline="" #input region outline defined by polygon: none for Venezuela |
|
| 103 |
#This is the shape file of outline of the study area #It is an input/output of the covariate script |
|
| 104 |
#infile_reg_outline <- "OR83M_state_outline.shp" #input region outline defined by polygon: Oregon |
|
| 105 |
#infile_reg_outline <-"OR83M_state_outline.shp" #remove this parameter!!! |
|
| 106 |
#covar_names see stage 2 |
|
| 107 |
|
|
| 108 |
#lc_path<-"/home/layers/data/land-cover/lc-consensus-global" |
|
| 109 |
#infile_modis_grid<-"/home/layers/commons/modis/modis_sinusoidal/modis_sinusoidal_grid_world.shp" #Give path!!! NEED TO CHANGE THIS... |
|
| 110 |
#infile_elev<-"/home/layers/data/terrain/dem-cgiar-srtm-1km-tif/srtm_1km.tif" #this is the global file: replace later with the input produced by the DEM team |
|
| 111 |
#infile_canheight<-"/home/layers/data/land-cover/treeheight-simard2011/Simard_Pinto_3DGlobalVeg_JGR.tif" #Canopy height |
|
| 112 |
|
|
| 94 | 113 |
list_tiles_modis <- c("h11v08,h11v07,h12v07,h12v08,h10v07,h10v08") #tile for Venezuela and surrounding area
|
| 95 | 114 |
#list_tiles_modis <- c("h08v04","h09v04") #tiles for Oregon
|
| 96 | 115 |
|
| ... | ... | |
| 101 | 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"; |
| 102 | 121 |
CRS_locs_WGS84<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84
|
| 103 | 122 |
out_region_name<-"_venezuela_region" #generated on the fly |
| 104 |
out_suffix<-"_VE_05242013" |
|
| 105 |
ref_rast_name<-"" #local raster name defining resolution, exent, local projection--. set on the fly?? |
|
| 123 |
#out_suffix<-"_VE_05252013" |
|
| 124 |
ref_rast_name<-"" #local raster name defining resolution, exent, local projection--. set on the fly?? |
|
| 125 |
#this may be redundant with infile_reg_outline |
|
| 106 | 126 |
#ref_rast_name<-"mean_day244_rescaled.rst" #local raster name defining resolution, exent: oregon |
| 107 | 127 |
|
| 108 | 128 |
#The names of covariates can be changed...these names should be output/input from covar script!!! |
| ... | ... | |
| 122 | 142 |
end_year = "2010" |
| 123 | 143 |
#end_year = "2002" #for testing (year included?) |
| 124 | 144 |
#hdfdir = '/home/layers/commons/modis/MOD11A1_tiles' #destination file where hdf files are stored locally after download. |
| 125 |
hdfdir = '/home/parmentier/Data/IPLANT_project/MOD11A1_tiles' |
|
| 145 |
#hdfdir = '/home/parmentier/Data/IPLANT_project/MOD11A1_tiles' |
|
| 146 |
hdfdir <- "/data/project/layers/commons/modis/MOD11A1_tiles" |
|
| 126 | 147 |
download=1 |
| 127 | 148 |
clim_calc=0 |
| 128 |
out_suffix_modis="_05242013"
|
|
| 149 |
#out_suffix_modis="_05252013"
|
|
| 129 | 150 |
#end_month= "12" |
| 130 | 151 |
#start_month= "1" |
| 131 | 152 |
|
| ... | ... | |
| 145 | 166 |
|
| 146 | 167 |
############ STAGE 2: Covariate production ################ |
| 147 | 168 |
|
| 169 |
#list of 17 parameters |
|
| 148 | 170 |
list_param_covar_production<-list(var,in_path,out_path,lc_path,infile_modis_grid,infile_elev,infile_canheight, |
| 149 |
list_tiles_modis,infile_reg_outline,CRS_interp,CRS_locs_WGS84,out_region_name, |
|
| 150 |
out_suffix,ref_rast_name,hdfdir,covar_names) |
|
| 171 |
infile_distoc,list_tiles_modis,infile_reg_outline,CRS_interp,CRS_locs_WGS84,out_region_name,
|
|
| 172 |
out_suffix,out_suffix_modis,ref_rast_name,hdfdir,covar_names)
|
|
| 151 | 173 |
|
| 152 | 174 |
names(list_param_covar_production)<-c("var","in_path","out_path","lc_path","infile_modis_grid","infile_elev","infile_canheight",
|
| 153 |
"list_tiles_modis","infile_reg_outline","CRS_interp","CRS_locs_WGS84","out_region_name", |
|
| 154 |
"out_suffix","ref_rast_name","hdfdir","covar_names") |
|
| 175 |
"infile_distoc","list_tiles_modis","infile_reg_outline","CRS_interp","CRS_locs_WGS84","out_region_name",
|
|
| 176 |
"out_suffix","out_suffix_modis","ref_rast_name","hdfdir","covar_names")
|
|
| 155 | 177 |
|
| 156 | 178 |
## Modify to store infile_covar_brick in output folder!!! |
| 157 | 179 |
if (stages_to_run[2]==2){
|
| 158 |
#Transform into function... |
|
| 159 |
#/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/covariates_production_temperatures_03212013.R |
|
| 160 |
infile_covar_brick<-covariates_production_temperature(list_param_covar_production) |
|
| 161 |
}else{
|
|
| 162 |
#Provide brick of covariates if stage 2 is not run |
|
| 163 |
infile_covariates<-"/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/covariates__venezuela_region__VE_01292013.tif" #this is an output from covariate script and used in stage 3 and stage 4 |
|
| 164 |
infile_reg_outline<- "/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/outline_venezuela_region__VE_01292013.shp" |
|
| 165 |
#infile_covariates<-"/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/covariates__venezuela_region_TMIN__VE_03192013.tif" #covariates stack for TMIN |
|
| 166 |
#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 |
|
| 180 |
covar_obj <- covariates_production_temperature(list_param_covar_production) |
|
| 181 |
infile_covariates <- covar_obj$infile_covariates |
|
| 182 |
infile_reg_outline <- covar_obj$infile_reg_outline |
|
| 167 | 183 |
} |
| 168 | 184 |
|
| 185 |
#Note that if stages_to_run[2]!=2, then use values defined at the beginning of the script for infile_covariates and infile_reg_outline |
|
| 186 |
|
|
| 169 | 187 |
############# STAGE 3: Data preparation ############### |
| 170 | 188 |
|
| 171 |
#infile_reg_outline="" #input region outline defined by polygon: none for Venezuela |
|
| 172 |
#This is the shape file of outline of the study area #It is an input/output of the covariate script |
|
| 173 |
#infile_reg_outline <- "OR83M_state_outline.shp" #input region outline defined by polygon: Oregon |
|
| 174 |
#infile_reg_outline <-"OR83M_state_outline.shp" #remove this parameter!!! |
|
| 175 |
#covar_names see stage 2 |
|
| 176 | 189 |
|
| 177 | 190 |
#specific to this stage |
| 178 | 191 |
db.name <- "ghcn" # name of the Postgres database |
| ... | ... | |
| 180 | 193 |
range_years_clim<-c("2000","2011") #right bound not included in the range!!
|
| 181 | 194 |
infile_ghncd_data <-"/home/layers/data/climate/ghcn/v2.92-upd-2012052822/ghcnd-stations.txt" #This is the textfile of station locations from GHCND |
| 182 | 195 |
qc_flags_stations<-c("0","S") #flags allowed for screening after the query from the GHCND??
|
| 196 |
#infile_covariates and infile_reg_outline defined in stage 2 or at the start of script... |
|
| 183 | 197 |
|
| 184 | 198 |
#list of 12 parameters for input in the function... |
| 185 | 199 |
|
| ... | ... | |
| 227 | 241 |
"y_var ~ s(lat,lon) + s(elev_s) + s(N_w,E_w) + s(LST) + s(LC6)") |
| 228 | 242 |
# "y_var ~ s(lat,lon) + s(elev_s) + s(N_w,E_w) + s(LST) + s(DISTOC)") |
| 229 | 243 |
|
| 230 |
|
|
| 231 | 244 |
#Default name of LST avg to be matched |
| 232 | 245 |
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")
|
| 233 | 246 |
|
Also available in: Unified diff
master script, covariates production automation first changes