Revision 9af28b45
Added by Benoit Parmentier over 11 years ago
| climate/research/oregon/interpolation/covariates_production_temperatures.R | ||
|---|---|---|
| 50 | 50 |
infile5<-"Simard_Pinto_3DGlobalVeg_JGR.tif" #Canopy height |
| 51 | 51 |
list_tiles_modis = c('h11v08','h11v07','h12v07','h12v08','h10v07','h10v08') #tile for Venezuel and surrounding area
|
| 52 | 52 |
infile_reg_outline="" #input region outline defined by polygon |
| 53 |
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";
|
|
| 53 |
CRS_interp<-"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs";
|
|
| 54 | 54 |
CRS_locs_WGS84<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84
|
| 55 | 55 |
out_region_name<-"_venezuela_region" |
| 56 |
out_suffix<-"_VE_01292013"
|
|
| 56 |
out_suffix<-"_VE_02082013"
|
|
| 57 | 57 |
ref_rast_name<-"" #local raster name defining resolution, exent, local projection--. set on the fly?? |
| 58 | 58 |
#for the processing tile/region? This is a group fo six tiles for now. |
| 59 | 59 |
|
| ... | ... | |
| 204 | 204 |
#"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" ?? |
| 205 | 205 |
#reassign proj from modis tile to raster? there is a 10m diff in semi-axes...(a and b) |
| 206 | 206 |
|
| 207 |
#Screen LST for extreme values? |
|
| 208 |
#min_val<-(-15+273.16) #if values less than -15C then screen out (note the Kelvin units that will need to be changed later in all datasets) |
|
| 209 |
#LST[LST < (min_val)]<-NA |
|
| 210 |
|
|
| 207 | 211 |
|
| 208 | 212 |
######################################### |
| 209 | 213 |
##2) Crop and reproject Canopy height data |
| ... | ... | |
| 275 | 279 |
lc_reg_list[[i]]<-file.path(out_path,raster_name) |
| 276 | 280 |
} |
| 277 | 281 |
setwd(out_path) |
| 278 |
lc_reg_list<-mixedsort(list.files(pattern="reg_con.*.tif")) |
|
| 282 |
lc_reg_list<-mixedsort(list.files(pattern="^reg_con.*.tif"))
|
|
| 279 | 283 |
lc_reg_s<-stack(lc_reg_list) |
| 280 | 284 |
#lc_reg_s<-as.character(lc_reg_list) |
| 281 | 285 |
#Now combine forest classes...in LC1 forest, LC2, LC3, LC4 and LC6-urban...?? |
| 282 | 286 |
|
| 283 | 287 |
#create a local mask for the tile/processing region |
| 284 | 288 |
|
| 285 |
LC12<-raster(paste("reg_con_1km_class_12_",out_suffix,".tif",sep="")) #this is open water
|
|
| 289 |
#LC12<-raster(paste("reg_con_1km_class_12_",out_suffix,".tif",sep="")) #this is open water
|
|
| 290 |
LC12<-raster(lc_reg_s,layer=nlayers(lc_reg_s)) #this is open water |
|
| 291 |
|
|
| 286 | 292 |
LC_mask<-LC12 |
| 287 | 293 |
LC_mask[LC_mask==100]<-NA |
| 288 | 294 |
LC_mask <- LC_mask > 100 |
| ... | ... | |
| 339 | 345 |
"nobs_09","nobs_10","nobs_11","nobs_12") |
| 340 | 346 |
names(lst_s)<-lst_names |
| 341 | 347 |
s_raster<-addLayer(s_raster, lst_s) |
| 342 |
s_raster<-mask(s_raster,LC_mask,filename="test.tif") |
|
| 343 | 348 |
|
| 344 | 349 |
covar_names<-c(rnames,lc_names,lst_names) |
| 345 | 350 |
names(s_raster)<-covar_names |
| 346 | 351 |
#Write out stack of number of change |
| 347 | 352 |
data_name<-paste("covariates_",out_region_name,"_",sep="")
|
| 348 | 353 |
raster_name<-paste(data_name,out_suffix,".tif", sep="") |
| 349 |
writeRaster(s_raster, filename=raster_name,NAflag=-999,bylayer=FALSE,bandorder="BSQ",overwrite=TRUE) #Writing the data in a raster file format... |
|
| 354 |
#writeRaster(s_raster, filename=raster_name,NAflag=-999,bylayer=FALSE,bandorder="BSQ",overwrite=TRUE) #Writing the data in a raster file format... |
|
| 355 |
s_raster<-mask(s_raster,LC_mask,filename=raster_name, |
|
| 356 |
overwrite=TRUE,NAflag=-999,bylayer=FALSE,bandorder="BSQ") |
|
| 350 | 357 |
#using bil format more efficient?? |
| 351 | 358 |
|
| 352 | 359 |
####################################################### |
Also available in: Unified diff
Data covariates preparation, slight modifications and updates