Revision 1a94c1db
Added by Benoit Parmentier over 9 years ago
| climate/research/oregon/interpolation/global_run_scalingup_mosaicing.R | ||
|---|---|---|
| 5 | 5 |
#Analyses, figures, tables and data are also produced in the script. |
| 6 | 6 |
#AUTHOR: Benoit Parmentier |
| 7 | 7 |
#CREATED ON: 04/14/2015 |
| 8 |
#MODIFIED ON: 04/14/2015
|
|
| 8 |
#MODIFIED ON: 05/07/2015
|
|
| 9 | 9 |
#Version: 4 |
| 10 | 10 |
#PROJECT: Environmental Layers project |
| 11 | 11 |
#COMMENTS: analyses for run 10 global analyses,all regions 1500x4500km and other tiles |
| ... | ... | |
| 267 | 267 |
lf_mosaic_pred_1500x4500 <-list.files(path=file.path(in_dir), |
| 268 | 268 |
pattern=paste(".*.tif$",sep=""),full.names=T)
|
| 269 | 269 |
|
| 270 |
r1 <- raster(lf_mosaic_pred_1500x4500[11])
|
|
| 270 |
r1 <- raster(lf_mosaic_pred_1500x4500[1]) |
|
| 271 | 271 |
r2 <- raster(lf_mosaic_pred_1500x4500[2]) |
| 272 | 272 |
|
| 273 |
assignVal<- function(x){rep(1,x)}
|
|
| 274 |
r_outline <- init(r2,fun=assignVal,filename="test.rst",overwrite=T) |
|
| 275 |
r_pol <- rasterToPolygons(r_outline,dissolve=T) |
|
| 276 | 273 |
|
| 274 |
lf <- sub(".tif","",lf_mosaic_pred_1500x4500)
|
|
| 275 |
tx<-strsplit(as.character(lf),"_") |
|
| 276 |
|
|
| 277 |
lat<- as.character(lapply(1:length(tx),function(i,x){x[[i]][13]},x=tx))
|
|
| 278 |
long<- as.character(lapply(1:length(tx),function(i,x){x[[i]][14]},x=tx))
|
|
| 279 |
lat <- as.character(lapply(1:length(lat),function(i,x){substr(x[[i]],2,nchar(x[i]))},x=lat)) #first number not in the coordinates
|
|
| 280 |
|
|
| 281 |
df_centroids <- data.frame(long=as.numeric(long),lat=as.numeric(lat)) |
|
| 282 |
df_centroids$ID <- as.numeric(1:nrow(df_centroids)) |
|
| 283 |
# |
|
| 284 |
extract(r1,) |
|
| 285 |
coordinates(df_centroids) <- cbind(df_centroids$long,df_centroids$lat) |
|
| 286 |
proj4string(df_centroids) <- projection(r1) |
|
| 277 | 287 |
#centroid distance |
| 278 |
c1 <- gCentroid(x,byid=TRUE) |
|
| 279 |
pt <- gCentroid(shp1) |
|
| 288 |
#c1 <- gCentroid(x,byid=TRUE)
|
|
| 289 |
#pt <- gCentroid(shp1)
|
|
| 280 | 290 |
|
| 291 |
#Make this a function later... |
|
| 281 | 292 |
#then distance... |
| 282 |
gDistance |
|
| 293 |
out_dir_str <- out_dir |
|
| 294 |
lf_r_weights <- vector("list",length=length(lf_mosaic_pred_1500x4500))
|
|
| 295 |
|
|
| 296 |
for(i in 1:length(lf)){
|
|
| 297 |
# |
|
| 298 |
r1 <- raster(lf_mosaic_pred_1500x4500[i]) |
|
| 299 |
|
|
| 300 |
set1f <- function(x){rep(NA, x)}
|
|
| 301 |
r_init <- init(r1, fun=set1f, filename='test.grd', overwrite=TRUE) |
|
| 283 | 302 |
|
| 303 |
cell_ID <- cellFromXY(r_init,xy=df_centroids[i,]) |
|
| 304 |
r_init[cell_ID] <- df_centroids$ID[i] |
|
| 305 |
|
|
| 306 |
r_dist <- distance(r_init) |
|
| 307 |
min_val <- cellStats(r_dist,min) |
|
| 308 |
max_val <- cellStats(r_dist,max) |
|
| 309 |
r_dist_normalized <- abs(r_dist - max_val)/ (max_val - min_val) |
|
| 310 |
|
|
| 311 |
extension_str <- extension(lf_mosaic_pred_1500x4500[i]) |
|
| 312 |
raster_name_tmp <- gsub(extension_str,"",basename(lf_mosaic_pred_1500x4500[i])) |
|
| 313 |
raster_name <- file.path(out_dir_str,paste(raster_name_tmp,"_weights.tif",sep="")) |
|
| 314 |
writeRaster(r_dist_normalized, NAflag=NA_flag_val,filename=raster_name,overwrite=TRUE) |
|
| 315 |
|
|
| 316 |
lf_r_weights[[i]] <- raster_name |
|
| 317 |
} |
|
| 284 | 318 |
|
| 285 | 319 |
#Then scale on 1 to zero? or 0 to 1000 |
| 286 | 320 |
#e.g. for a specific pixel |
| 287 | 321 |
#weight_sum=0.2 +0.4 +0.4+0.2=1.2 |
| 288 |
#sum_val (0.2*val1+0.4*val2+0.4*val3+0.2*val4) |
|
| 289 |
#m_val= /weight_sum |
|
| 322 |
#val_w_sum= (0.2*val1+0.4*val2+0.4*val3+0.2*val4) |
|
| 323 |
#no_valid = |
|
| 324 |
#m_val= sum_val/weight_sum #mean value |
|
| 325 |
# |
|
| 326 |
|
|
| 327 |
#in raster term |
|
| 328 |
#r_weights_sum <- ... |
|
| 329 |
#r_val_w_sum <- |
|
| 330 |
# |
|
| 290 | 331 |
|
| 332 |
|
|
| 333 |
r1 <- raster(c) |
|
| 291 | 334 |
#can do mosaic with sum?? for both weighted sum and val |
| 292 | 335 |
# |
| 293 | 336 |
#can use gdal calc... |
| 294 | 337 |
|
| 295 |
r_m <- r1 + r2 |
|
| 296 |
name_method <- paste(interpolation_method,"_",y_var_name,"_",sep="") |
|
| 338 |
#r_m <- r1 + r2
|
|
| 339 |
#name_method <- paste(interpolation_method,"_",y_var_name,"_",sep="")
|
|
| 297 | 340 |
##Use python code written by Alberto Guzman |
| 298 | 341 |
|
| 299 | 342 |
#system("MODULEPATH=$MODULEPATH:/nex/modules/files")
|
| ... | ... | |
| 302 | 345 |
lf2 <- lf_world_pred_1500x4500 |
| 303 | 346 |
|
| 304 | 347 |
#module_path <- "" |
| 305 |
module_path <- "/nobackupp6/aguzman4/climateLayers/sharedCode/" |
|
| 348 |
#module_path <- "/nobackupp6/aguzman4/climateLayers/sharedCode/"
|
|
| 306 | 349 |
#sh /nobackupp6/aguzman4/climateLayers/sharedCode/gdalCalDiff.sh file1.tif file2.tif output.tif |
| 307 | 350 |
#/nobackupp6/aguzman4/climateLayers/sharedCode/mosaicUsingGdalMerge.py |
| 308 | 351 |
#l_dates <- paste(day_to_mosaic,collapse=",",sep=" ") |
| 309 |
l_dates <- paste(day_to_mosaic,collapse=",") |
|
| 352 |
#l_dates <- paste(day_to_mosaic,collapse=",")
|
|
| 310 | 353 |
## use region 2 first |
| 311 |
lf_out <- paste("diff_world_","1000_3000","by1500_4500_","mod1","_",l_dates,out_suffix,"_",file_format,sep="")
|
|
| 312 |
|
|
| 313 |
|
|
| 314 |
for (i in 1:length(lf_out)){
|
|
| 315 |
out_file <- lf_out[i] |
|
| 316 |
in_file1 <- lf1[i] |
|
| 317 |
in_file2 <- lf2[i] |
|
| 318 |
|
|
| 319 |
cmd_str <- paste("sh", file.path(module_path,"gdalCalDiff.sh"),
|
|
| 320 |
in_file1, |
|
| 321 |
in_file2, |
|
| 322 |
out_file,sep=" ") |
|
| 323 |
system(cmd_str) |
|
| 324 |
|
|
| 325 |
} |
|
| 354 |
#lf_out <- paste("diff_world_","1000_3000","by1500_4500_","mod1","_",l_dates,out_suffix,"_",file_format,sep="")
|
|
| 355 |
|
|
| 356 |
|
|
| 357 |
#for (i in 1:length(lf_out)){
|
|
| 358 |
# out_file <- lf_out[i]
|
|
| 359 |
# in_file1 <- lf1[i]
|
|
| 360 |
# in_file2 <- lf2[i]
|
|
| 361 |
# |
|
| 362 |
# cmd_str <- paste("sh", file.path(module_path,"gdalCalDiff.sh"),
|
|
| 363 |
# in_file1,
|
|
| 364 |
# in_file2,
|
|
| 365 |
# out_file,sep=" ")
|
|
| 366 |
# system(cmd_str)
|
|
| 367 |
# |
|
| 368 |
#}
|
|
| 326 | 369 |
|
| 327 | 370 |
##################### END OF SCRIPT ###################### |
Also available in: Unified diff
scaling up mosacing experiment, edits