Revision c79a8453
Added by Benoit Parmentier over 9 years ago
| climate/research/oregon/interpolation/global_run_scalingup_mosaicing.R | ||
|---|---|---|
| 319 | 319 |
|
| 320 | 320 |
#on ATLAS |
| 321 | 321 |
|
| 322 |
in_dir <- "~/Data/IPLANT_project/mosaicing_data_test/reg2" |
|
| 322 |
#in_dir <- "~/Data/IPLANT_project/mosaicing_data_test/reg2" |
|
| 323 |
in_dir <- "~/Data/IPLANT_project/mosaicing_data_test/reg1" |
|
| 323 | 324 |
|
| 324 | 325 |
y_var_name <- "dailyTmax" #PARAM1 |
| 325 | 326 |
interpolation_method <- c("gam_CAI") #PARAM2
|
| ... | ... | |
| 374 | 375 |
|
| 375 | 376 |
r1 <- raster(lf_mosaic_pred_1500x4500[1]) |
| 376 | 377 |
r2 <- raster(lf_mosaic_pred_1500x4500[2]) |
| 377 |
|
|
| 378 |
r5 |
|
| 379 |
plot(r1) |
|
| 380 |
plot(r2) |
|
| 378 | 381 |
|
| 379 | 382 |
lf <- sub(".tif","",lf_mosaic_pred_1500x4500)
|
| 380 | 383 |
tx<-strsplit(as.character(lf),"_") |
| ... | ... | |
| 386 | 389 |
df_centroids <- data.frame(long=as.numeric(long),lat=as.numeric(lat)) |
| 387 | 390 |
df_centroids$ID <- as.numeric(1:nrow(df_centroids)) |
| 388 | 391 |
# |
| 389 |
extract(r1,) |
|
| 392 |
#extract(r1,)
|
|
| 390 | 393 |
coordinates(df_centroids) <- cbind(df_centroids$long,df_centroids$lat) |
| 391 | 394 |
proj4string(df_centroids) <- projection(r1) |
| 392 | 395 |
#centroid distance |
| ... | ... | |
| 403 | 406 |
num_cores <- 6 |
| 404 | 407 |
|
| 405 | 408 |
#debug(create_weights_fun) |
| 406 |
weights_obj <- create_weights_fun(1,list_param=list_param_create_weights) |
|
| 409 |
#weights_obj <- create_weights_fun(1,list_param=list_param_create_weights)
|
|
| 407 | 410 |
|
| 408 | 411 |
weights_obj_list <- mclapply(1:length(lf_mosaic_pred_1500x4500),FUN=create_weights_fun,list_param=list_param_create_weights,mc.preschedule=FALSE,mc.cores = num_cores) |
| 409 | 412 |
|
| 410 |
list_args_weights_prod <- lapply(1:length(weights_obj_list), FUN=function(x){raster(x[[i]]$r_weights_prod})}
|
|
| 411 |
list_args_weights_prod$fun <- "sum" |
|
| 413 |
#list_args_weights_prod <- lapply(1:length(weights_obj_list), FUN=function(x){raster(x[[i]]$r_weights_prod})}
|
|
| 414 |
#list_args_weights_prod$fun <- "sum"
|
|
| 412 | 415 |
|
| 413 | 416 |
|
| 414 | 417 |
#"r_weights","r_weights_prod" |
| ... | ... | |
| 419 | 422 |
list_r_weights <- lapply(1:length(weights_obj_list), FUN=function(i,x){x[[i]]$r_weights},x=weights_obj_list)
|
| 420 | 423 |
list_r_weights_prod <- lapply(1:length(weights_obj_list), FUN=function(i,x){x[[i]]$r_weights_prod},x=weights_obj_list)
|
| 421 | 424 |
|
| 422 |
list_args_weights <- list_r_weights |
|
| 423 |
list_args_weights_prod <- list_r_weights_prod |
|
| 424 |
|
|
| 425 |
list_args_weights$fun <- "sum" |
|
| 426 |
#list_args_weights$fun <- "mean" |
|
| 427 |
|
|
| 428 |
list_args_weights$na.rm <- TRUE |
|
| 429 |
list_args_weights$tolerance <- 1 |
|
| 430 |
|
|
| 431 |
list_args_weights_prod$fun <- "sum" |
|
| 432 |
list_args_weights_prod$na.rm <- TRUE |
|
| 433 |
list_args_weights_prod$na.rm <- TRUE |
|
| 434 |
|
|
| 435 |
r_weights_sum <- do.call(mosaic,list_args_weights) #sum |
|
| 436 |
r_prod_sum <- do.call(mosaic,list_args_weights_prod) #sum |
|
| 437 |
|
|
| 438 |
r_weighted_mean <- r_weights_sum/r_prod_sum |
|
| 439 | 425 |
|
| 440 | 426 |
|
| 441 | 427 |
#r_weights_sum <- do.call(overlay,list_args_weights) #sum |
| ... | ... | |
| 462 | 448 |
#r_val_w_sum <- |
| 463 | 449 |
# |
| 464 | 450 |
mosaic_list_var <- list(list_r_weights) |
| 451 |
mosaic_list_var <- list(lf_mosaic_pred_1500x4500) |
|
| 465 | 452 |
#out_rastnames_var <- l_out_rastnames_var[[i]] |
| 466 |
out_rastnames_var <- c("reg2_mosaic_weights.tif")
|
|
| 453 |
out_rastnames_var <- c("reg1_mosaic_mean.tif")
|
|
| 467 | 454 |
|
| 468 | 455 |
#list_param_mosaic <- list(list_r_weights,out_dir,outrastnames,file_format,NA_flag_val,out_suffix) |
| 469 | 456 |
|
| ... | ... | |
| 473 | 460 |
j<-1 #date index for loop |
| 474 | 461 |
list_param_mosaic<-list(j,mosaic_list_var,out_rastnames_var,out_dir,file_format,NA_flag_val) |
| 475 | 462 |
names(list_param_mosaic)<-c("j","mosaic_list","out_rastnames","out_path","file_format","NA_flag_val")
|
| 476 |
debug(mosaic_m_raster_list) |
|
| 477 |
mosaic_m_raster_list(1,list_param_mosaic) |
|
| 463 |
#undebug(mosaic_m_raster_list)
|
|
| 464 |
r_mosaiced <- mosaic_m_raster_list(1,list_param_mosaic)
|
|
| 478 | 465 |
|
| 479 | 466 |
|
| 480 | 467 |
#list_var_mosaiced <- mclapply(1:2,FUN=mosaic_m_raster_list,list_param=list_param_mosaic,mc.preschedule=FALSE,mc.cores = 2) |
| ... | ... | |
| 482 | 469 |
#list_var_mosaiced <- mclapply(1:1,FUN=mosaic_m_raster_list,list_param=list_param_mosaic,mc.preschedule=FALSE,mc.cores = 1) |
| 483 | 470 |
#list_var_mosaiced <- mclapply(1:365,FUN=mosaic_m_raster_list,list_param=list_param_mosaic,mc.preschedule=FALSE,mc.cores = 2) |
| 484 | 471 |
|
| 485 |
outrastnames <- "reg2_mosaic_weights.tif"
|
|
| 472 |
outrastnames <- "reg1_mosaic_weights.tif"
|
|
| 486 | 473 |
|
| 487 | 474 |
list_param_mosaic <- list(list_r_weights,out_dir,outrastnames,file_format,NA_flag_val,out_suffix) |
| 488 | 475 |
|
| 489 |
#mosaic_list<-list_param$mosaic_list |
|
| 490 |
#out_path<-list_param$out_path |
|
| 491 |
# out_names<-list_param$out_rastnames |
|
| 492 |
# file_format <- list_param$file_format |
|
| 493 |
# NA_flag_val <- list_param$NA_flag_val |
|
| 494 |
# out_suffix <- list_param$out_suffix |
|
| 476 |
r1_projected <- projectRaster(raster(list_r_weights[[1]]),r) |
|
| 477 |
|
|
| 478 |
cmd_str <- paste("python","/usr/bin/gdal_merge.py","-o avg.tif",paste(lf_mosaic_pred_1500x4500,collapse=" "))
|
|
| 479 |
system(cmd_str) |
|
| 480 |
|
|
| 481 |
lf_files <- list_r_weights |
|
| 482 |
r_m <- raster("avg.tif")
|
|
| 483 |
|
|
| 484 |
for (i in 1:length(lf_files)){
|
|
| 485 |
set1f <- function(x){rep(NA, x)}
|
|
| 486 |
r_ref <- init(r_m, fun=set1f, filename=paste('r_weigts_m_',i,'.tif',sep=""), overwrite=TRUE)
|
|
| 487 |
#NAvalue(r_ref) <- -9999 |
|
| 488 |
|
|
| 489 |
cmd_str <- paste("/usr/bin/gdalwarp",list_r_weights[[i]],paste('r_weights_m_',i,'.tif',sep=""),sep=" ")
|
|
| 490 |
#gdalwarp -t_srs '+proj=utm +zone=11 +datum=WGS84' raw_spot.tif utm11.tif |
|
| 491 |
|
|
| 492 |
system(cmd_str) |
|
| 493 |
#r_ref <-raster("r_ref.tif")
|
|
| 494 |
#plot(r_ref) |
|
| 495 |
} |
|
| 496 |
|
|
| 497 |
list_args_weights <- (mixedsort(list.files(pattern="r_weigts_m_.*.tif"))) |
|
| 498 |
list_args_weights <- lapply(1:length(list_args_weights), FUN=function(i,x){raster(x[[i]])},x=list_args_weights)
|
|
| 499 |
|
|
| 500 |
#list_args_weights_prod <- list_r_weights_prod |
|
| 501 |
|
|
| 502 |
list_args_weights$fun <- "sum" |
|
| 503 |
#list_args_weights$fun <- "mean" |
|
| 504 |
|
|
| 505 |
list_args_weights$na.rm <- TRUE |
|
| 506 |
#list_args_weights$tolerance <- 1 |
|
| 507 |
|
|
| 508 |
#l#ist_args_weights_prod$fun <- "sum" |
|
| 509 |
#list_args_weights_prod$na.rm <- TRUE |
|
| 510 |
#list_args_weights_prod$na.rm <- TRUE |
|
| 511 |
|
|
| 512 |
r_weights_sum <- do.call(mosaic,list_args_weights) #sum |
|
| 513 |
#r_prod_sum <- do.call(mosaic,list_args_weights_prod) #sum |
|
| 514 |
|
|
| 515 |
r_weighted_mean <- r_weights_sum/r_prod_sum |
|
| 495 | 516 |
|
| 496 | 517 |
################################################# |
| 497 | 518 |
#Ok testing on fake data: |
Also available in: Unified diff
scaling up mosaicing experiment, using gdal_merge and gdalwarp to resolve mosaicing issue