Revision 70d0063c
Added by Benoit Parmentier almost 12 years ago
| climate/research/oregon/interpolation/covariates_production_temperatures.R | ||
|---|---|---|
| 10 | 10 |
# -MODIS LST: mean and obs |
| 11 | 11 |
#3) The output is a multiband file in tif format with projected covariates for the processing region/tile. |
| 12 | 12 |
#AUTHOR: Benoit Parmentier |
| 13 |
#DATE: 03/19/2013
|
|
| 13 |
#DATE: 03/21/2013
|
|
| 14 | 14 |
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#363-- |
| 15 | 15 |
|
| 16 | 16 |
##Comments and TODO: |
| ... | ... | |
| 22 | 22 |
|
| 23 | 23 |
################################################################################################## |
| 24 | 24 |
|
| 25 |
###Loading R library and packages |
|
| 26 |
library(RPostgreSQL) |
|
| 27 |
library(sp) # Spatial pacakge with class definition by Bivand et al. |
|
| 28 |
library(spdep) # Spatial pacakge with methods and spatial stat. by Bivand et al. |
|
| 29 |
library(rgdal) # GDAL wrapper for R, spatial utilities |
|
| 30 |
library(raster) |
|
| 31 |
library(gtools) |
|
| 32 |
library(rasterVis) |
|
| 33 |
library(graphics) |
|
| 34 |
library(grid) |
|
| 35 |
library(lattice) |
|
| 36 |
|
|
| 37 |
### Parameters and arguments |
|
| 38 |
|
|
| 39 |
##Paths to inputs and output |
|
| 40 |
var<-"TMAX" |
|
| 41 |
in_path <- "/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/input_data/" |
|
| 42 |
out_path<- "/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/output_data/" |
|
| 43 |
lc_path<-"/home/layers/data/land-cover/lc-consensus-global" |
|
| 44 |
#elev_path<-"/home/layers/data/terrain/dem-cgiar-srtm-1km-tif" |
|
| 45 |
#infile3<-"countries_sinusoidal_world.shp" |
|
| 46 |
#infile1<-"worldborder_sinusoidal.shp" |
|
| 47 |
infile_modis_grid<-"modis_sinusoidal_grid_world.shp" |
|
| 48 |
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 |
|
| 49 |
infile_canheight<-"Simard_Pinto_3DGlobalVeg_JGR.tif" #Canopy height |
|
| 50 |
list_tiles_modis = c('h11v08','h11v07','h12v07','h12v08','h10v07','h10v08') #tile for Venezuel and surrounding area
|
|
| 51 |
infile_reg_outline="" #input region outline defined by polygon |
|
| 52 |
CRS_interp<-"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"; |
|
| 53 |
CRS_locs_WGS84<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84
|
|
| 54 |
out_region_name<-"_venezuela_region" #generated on the fly |
|
| 55 |
out_suffix<-"_VE_03182013" |
|
| 56 |
ref_rast_name<-"" #local raster name defining resolution, exent, local projection--. set on the fly?? |
|
| 57 |
#for the processing tile/region? This is a group fo six tiles for now. |
|
| 58 |
|
|
| 59 |
#The names of covariates can be changed...these names should be output/input from covar script!!! |
|
| 60 |
rnames<-c("x","y","lon","lat","N","E","N_w","E_w","elev","slope","aspect","CANHEIGHT","DISTOC")
|
|
| 61 |
lc_names<-c("LC1","LC2","LC3","LC4","LC5","LC6","LC7","LC8","LC9","LC10","LC11","LC12")
|
|
| 62 |
lst_names<-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",
|
|
| 63 |
"nobs_01","nobs_02","nobs_03","nobs_04","nobs_05","nobs_06","nobs_07","nobs_08", |
|
| 64 |
"nobs_09","nobs_10","nobs_11","nobs_12") |
|
| 65 |
covar_names<-c(rnames,lc_names,lst_names) |
|
| 66 |
|
|
| 67 |
list_param_covar_production<-list(var,in_path,out_path,lc_path,infile_modis_grid,infile_elev,infile_canheight, |
|
| 68 |
list_tiles_modis,infile_reg_outline,CRS_interp,CRS_locs_WGS84,out_region_name, |
|
| 69 |
out_suffix,ref_rast_name,covar_names) |
|
| 70 |
|
|
| 71 |
names(list_param_covar_production)<-c("var","in_path","out_path","lc_path","infile_modis_grid","infile_elev","infile_canheight",
|
|
| 72 |
"list_tiles_modis","infile_reg_outline","CRS_interp","CRS_locs_WGS84","out_region_name", |
|
| 73 |
"out_suffix","ref_rast_name","covar_names") |
|
| 74 |
|
|
| 75 |
#LST_night_rast<-raster("mean_LST_Night_1km_h11v08_dec_11_03192013.tif")
|
|
| 76 |
#### Functions used in the script ### |
|
| 77 | 25 |
|
| 78 | 26 |
create_modis_tiles_region<-function(modis_grid,tiles){
|
| 79 | 27 |
#This functions returns a subset of tiles from the modis grdi. |
| ... | ... | |
| 145 | 93 |
# |
| 146 | 94 |
# |
| 147 | 95 |
# |
| 96 |
|
|
| 97 |
###Loading R library and packages |
|
| 98 |
library(RPostgreSQL) |
|
| 99 |
library(sp) # Spatial pacakge with class definition by Bivand et al. |
|
| 100 |
library(spdep) # Spatial pacakge with methods and spatial stat. by Bivand et al. |
|
| 101 |
library(rgdal) # GDAL wrapper for R, spatial utilities |
|
| 102 |
library(raster) |
|
| 103 |
library(gtools) |
|
| 104 |
library(rasterVis) |
|
| 105 |
library(graphics) |
|
| 106 |
library(grid) |
|
| 107 |
library(lattice) |
|
| 108 |
|
|
| 109 |
### Parameters and arguments |
|
| 110 |
|
|
| 148 | 111 |
########################################################### |
| 149 | 112 |
############ Main body: BEGIN--START OF THE SCRIPT ################### |
| 150 | 113 |
|
| 151 |
##### STEP 1: Reading region or tile information to set the study or processing region
|
|
| 114 |
##### STEP 1: PARSING ARGUMENTS
|
|
| 152 | 115 |
|
| 153 | 116 |
var<-list_param$var |
| 154 | 117 |
in_path <-list_param$in_path |
| 155 | 118 |
out_path<- list_param$out_path |
| 156 | 119 |
lc_path<-list_param$lc_path |
| 157 |
#elev_path<-"/home/layers/data/terrain/dem-cgiar-srtm-1km-tif" |
|
| 158 |
#infile3<-"countries_sinusoidal_world.shp" |
|
| 159 |
#infile1<-"worldborder_sinusoidal.shp" |
|
| 160 | 120 |
infile_modis_grid<-list_param$infile_modis_grid |
| 161 | 121 |
infile_elev<-list_param$infile_elev #this is the global file: replace later with the input produced by the DEM team |
| 162 | 122 |
infile_canheight<-list_param$infile_canheight #Canopy height |
| 163 |
list_tiles_modis<-list_param$list_tiles_modis #tile for Venezuel and surrounding area |
|
| 123 |
list_tiles_modis<-list_param$list_tiles_modis #tile for Venezuela and surrounding area
|
|
| 164 | 124 |
infile_reg_outline<-list_param$infile_reg_outline #input region outline defined by polygon |
| 165 | 125 |
CRS_interp<-list_param$CRS_interp #local projection system |
| 166 | 126 |
CRS_locs_WGS84<-list_param$CRS_locs_WGS84 # |
| 167 | 127 |
out_region_name<-list_param$out_region_name #generated on the fly |
| 168 | 128 |
out_suffix<-list_param$out_suffix |
| 169 | 129 |
ref_rast_name<-list_param$ref_rast_name #local raster name defining resolution, exent, local projection--. set on the fly?? |
| 170 |
#for the processing tile/region? This is a group fo six tiles for now. |
|
| 171 |
|
|
| 172 | 130 |
covar_names<-list_param$covar_names |
| 173 | 131 |
|
| 132 |
##### SET UP STUDY AREA #### |
|
| 133 |
|
|
| 174 | 134 |
setwd(in_path) |
| 175 | 135 |
|
| 176 | 136 |
filename<-sub(".shp","",infile_modis_grid) #Removing the extension from file.
|
| ... | ... | |
| 293 | 253 |
Nw<-sin(r2)*cos(r1) #Adding a variable to the dataframe |
| 294 | 254 |
Ew<-sin(r2)*sin(r1) #Adding variable to the dataframe. |
| 295 | 255 |
|
| 296 |
#topo_rast<-stack(STRM_reg,N,E,Nw,Ew) |
|
| 297 |
|
|
| 298 | 256 |
###################################### |
| 299 | 257 |
#4) LCC land cover |
| 300 | 258 |
|
| ... | ... | |
| 324 | 282 |
|
| 325 | 283 |
#LC12<-raster(paste("reg_con_1km_class_12_",out_suffix,".tif",sep="")) #this is open water
|
| 326 | 284 |
LC12<-raster(lc_reg_s,layer=nlayers(lc_reg_s)) #this is open water |
| 327 |
|
|
| 328 | 285 |
LC_mask<-LC12 |
| 329 | 286 |
LC_mask[LC_mask==100]<-NA |
| 330 | 287 |
LC_mask <- LC_mask > 100 |
| 331 |
#lc_reg_s<-mask(x=lc_reg_s,mask=LC_mask,filename=paste("reg_con_1km_all_classes_",out_suffix,".tif",sep=""),
|
|
| 332 |
# bandorder="BSQ",overwrite=TRUE) |
|
| 333 | 288 |
|
| 334 | 289 |
############################### |
| 335 | 290 |
#5) DISTOC, distance to coast: Would be useful to have a distance to coast layer ready... |
| ... | ... | |
| 365 | 320 |
|
| 366 | 321 |
################################ |
| 367 | 322 |
##Step 3: combine covariates in one stack for the next work flow stage |
| 368 |
#Create a stack in tif format... |
|
| 369 |
|
|
| 370 |
#? output name?? |
|
| 323 |
|
|
| 371 | 324 |
r<-stack(x,y,lon,lat,N,E,Nw,Ew,SRTM_reg,terrain_rast,CANHEIGHT,distoc_reg) |
| 372 | 325 |
#rnames<-c("x","y","lon","lat","N","E","N_w","E_w","elev","slope","aspect","CANHEIGHT","DISTOC")
|
| 373 |
#names(r)<-rnames |
|
| 374 | 326 |
s_raster<-r |
| 375 | 327 |
#Add landcover layers |
| 376 | 328 |
#lc_names<-c("LC1","LC2","LC3","LC4","LC5","LC6","LC7","LC8","LC9","LC10","LC11","LC12")
|
| 377 |
#names(lc_reg_s)<-lc_names #assign land cover names |
|
| 378 | 329 |
s_raster<-addLayer(s_raster, lc_reg_s) |
| 379 | 330 |
|
| 380 | 331 |
lst_s<-stack(c(as.character(mean_m_list),as.character(nobs_m_list))) |
| 381 |
#lst_names<-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",
|
|
| 382 |
# "nobs_01","nobs_02","nobs_03","nobs_04","nobs_05","nobs_06","nobs_07","nobs_08", |
|
| 383 |
# "nobs_09","nobs_10","nobs_11","nobs_12") |
|
| 384 |
#names(lst_s)<-lst_names |
|
| 332 |
|
|
| 385 | 333 |
s_raster<-addLayer(s_raster, lst_s) |
| 386 | 334 |
|
| 387 | 335 |
#covar_names<-c(rnames,lc_names,lst_names) |
| ... | ... | |
| 396 | 344 |
return(raster_name) |
| 397 | 345 |
} |
| 398 | 346 |
|
| 399 |
|
|
| 400 | 347 |
####################################################### |
| 401 |
################### END OF SCRIPT/FUNCTION ##################### |
|
| 348 |
################### END OF SCRIPT/FUNCTION ##################### |
|
Also available in: Unified diff
Covariates production raster modifications -transforming scrip into function