Revision 5e55c5eb
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: 05/11/2015
|
|
| 8 |
#MODIFIED ON: 05/18/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 |
| ... | ... | |
| 69 | 69 |
return(out_dir) |
| 70 | 70 |
} |
| 71 | 71 |
|
| 72 |
#Remove models that were not fitted from the list |
|
| 73 |
#All modesl that are "try-error" are removed |
|
| 74 |
remove_errors_list<-function(list_items){
|
|
| 75 |
|
|
| 76 |
#This function removes "error" items in a list |
|
| 77 |
list_tmp<-list_items |
|
| 78 |
if(is.null(names(list_tmp))){
|
|
| 79 |
names(list_tmp) <- paste("l",1:length(list_tmp),sep="_")
|
|
| 80 |
names(list_items) <- paste("l",1:length(list_tmp),sep="_")
|
|
| 81 |
} |
|
| 82 |
|
|
| 83 |
for(i in 1:length(list_items)){
|
|
| 84 |
if(inherits(list_items[[i]],"try-error")){
|
|
| 85 |
list_tmp[[i]]<-0 |
|
| 86 |
}else{
|
|
| 87 |
list_tmp[[i]]<-1 |
|
| 88 |
} |
|
| 89 |
} |
|
| 90 |
cnames<-names(list_tmp[list_tmp>0]) |
|
| 91 |
x <- list_items[match(cnames,names(list_items))] |
|
| 92 |
return(x) |
|
| 93 |
} |
|
| 94 |
|
|
| 95 |
|
|
| 96 |
plot_daily_mosaics <- function(i,list_param_plot_daily_mosaics){
|
|
| 97 |
#Purpose: |
|
| 98 |
#This functions mask mosaics files for a default range (-100,100 deg). |
|
| 99 |
#It produces a masked tif in a given dataType format (FLT4S) |
|
| 100 |
#It creates a figure of mosaiced region being interpolated. |
|
| 101 |
#Author: Benoit Parmentier |
|
| 102 |
#Parameters: |
|
| 103 |
#lf_m: list of files |
|
| 104 |
#reg_name:region name with tile size included |
|
| 105 |
#To do: |
|
| 106 |
#Add filenames |
|
| 107 |
#Add range |
|
| 108 |
#Add output dir |
|
| 109 |
#Add dataType_val option |
|
| 110 |
|
|
| 111 |
##### BEGIN ######## |
|
| 112 |
|
|
| 113 |
#Parse the list of parameters |
|
| 114 |
lf_m <- list_param_plot_daily_mosaics$lf_m |
|
| 115 |
reg_name <- list_param_plot_daily_mosaics$reg_name |
|
| 116 |
out_dir_str <- list_param_plot_daily_mosaics$out_dir_str |
|
| 117 |
out_suffix <- list_param_plot_daily_mosaics$out_suffix |
|
| 118 |
l_dates <- list_param_plot_daily_mosaics$l_dates |
|
| 119 |
|
|
| 120 |
|
|
| 121 |
#list_param_plot_daily_mosaics <- list(lf_m=lf_m,reg_name=reg_name,out_dir_str=out_dir_str) |
|
| 122 |
|
|
| 123 |
|
|
| 124 |
#rast_list <- vector("list",length=length(lf_m))
|
|
| 125 |
r_test<- raster(lf_m[i]) |
|
| 126 |
|
|
| 127 |
m <- c(-Inf, -100, NA, |
|
| 128 |
-100, 100, 1, |
|
| 129 |
100, Inf, NA) #can change the thresholds later |
|
| 130 |
rclmat <- matrix(m, ncol=3, byrow=TRUE) |
|
| 131 |
rc <- reclassify(r_test, rclmat,filename=paste("rc_tmp_",i,".tif",sep=""),dataType="FLT4S",overwrite=T)
|
|
| 132 |
file_name <- unlist(strsplit(basename(lf_m[i]),"_")) |
|
| 133 |
|
|
| 134 |
#date_proc <- file_name[7] #specific tot he current naming of files |
|
| 135 |
date_proc <- l_dates[i] |
|
| 136 |
#paste(raster_name[1:7],collapse="_") |
|
| 137 |
#add filename option later |
|
| 138 |
extension_str <- extension(filename(r_test)) |
|
| 139 |
raster_name_tmp <- gsub(extension_str,"",basename(filename(r_test))) |
|
| 140 |
raster_name <- file.path(out_dir_str,paste(raster_name_tmp,"_masked.tif",sep="")) |
|
| 141 |
r_pred <- mask(r_test,rc,filename=raster_name,overwrite=TRUE) |
|
| 142 |
|
|
| 143 |
res_pix <- 1200 |
|
| 144 |
#res_pix <- 480 |
|
| 145 |
|
|
| 146 |
col_mfrow <- 1 |
|
| 147 |
row_mfrow <- 1 |
|
| 148 |
|
|
| 149 |
png(filename=paste("Figure9_clim_mosaics_day_test","_",date_proc,"_",reg_name,"_",out_suffix,".png",sep=""),
|
|
| 150 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 151 |
|
|
| 152 |
plot(r_pred,main=paste("Predicted on ",date_proc , " ", reg_name,sep=""),cex.main=1.5)
|
|
| 153 |
dev.off() |
|
| 154 |
|
|
| 155 |
return(raster_name) |
|
| 156 |
|
|
| 157 |
} |
|
| 158 |
|
|
| 159 |
plot_screen_raster_val<-function(i,list_param){
|
|
| 160 |
##USAGE### |
|
| 161 |
#Screen raster list and produced plot as png. |
|
| 162 |
fname <-list_param$lf_raster_fname[i] |
|
| 163 |
screenRast <- list_param$screenRast |
|
| 164 |
l_dates<- list_param$l_dates |
|
| 165 |
out_dir_str <- list_param$out_dir_str |
|
| 166 |
prefix_str <-list_param$prefix_str |
|
| 167 |
out_suffix_str <- list_param$out_suffix_str |
|
| 168 |
|
|
| 169 |
### START SCRIPT #### |
|
| 170 |
date_proc <- l_dates[i] |
|
| 171 |
|
|
| 172 |
if(screenRast==TRUE){
|
|
| 173 |
r_test <- raster(fname) |
|
| 174 |
|
|
| 175 |
m <- c(-Inf, -100, NA, |
|
| 176 |
-100, 100, 1, |
|
| 177 |
100, Inf, NA) #can change the thresholds later |
|
| 178 |
rclmat <- matrix(m, ncol=3, byrow=TRUE) |
|
| 179 |
rc <- reclassify(r_test, rclmat,filename=paste("rc_tmp_",i,".tif",sep=""),dataType="FLT4S",overwrite=T)
|
|
| 180 |
#file_name <- unlist(strsplit(basename(lf_m[i]),"_")) |
|
| 181 |
extension_str <- extension(filename(r_test)) |
|
| 182 |
raster_name_tmp <- gsub(extension_str,"",basename(filename(r_test))) |
|
| 183 |
raster_name <- file.path(out_dir_str,paste(raster_name_tmp,"_masked.tif",sep="")) |
|
| 184 |
|
|
| 185 |
r_pred <- mask(r_test,rc,filename=raster_name,overwrite=TRUE) |
|
| 186 |
}else{
|
|
| 187 |
r_pred <- raster(fname) |
|
| 188 |
} |
|
| 189 |
|
|
| 190 |
#date_proc <- file_name[7] #specific tot he current naming of files |
|
| 191 |
#date_proc<- "2010010101" |
|
| 192 |
|
|
| 193 |
#paste(raster_name[1:7],collapse="_") |
|
| 194 |
#add filename option later |
|
| 195 |
|
|
| 196 |
res_pix <- 960 |
|
| 197 |
#res_pix <- 480 |
|
| 198 |
|
|
| 199 |
col_mfrow <- 1 |
|
| 200 |
row_mfrow <- 1 |
|
| 201 |
|
|
| 202 |
# png(filename=paste("Figure10_clim_world_mosaics_day_","_",date_proc,"_",tile_size,"_",out_suffix,".png",sep=""),
|
|
| 203 |
# width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 204 |
png(filename=paste(prefix_str,"_",date_proc,"_",tile_size,"_",out_suffix_str,".png",sep=""), |
|
| 205 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 206 |
|
|
| 207 |
plot(r_pred,main=paste("Predicted on ",date_proc , " ", tile_size,sep=""),cex.main=1.5)
|
|
| 208 |
dev.off() |
|
| 209 |
|
|
| 210 |
} |
|
| 211 |
|
|
| 212 | 72 |
create_weights_fun <- function(i, list_param){
|
| 213 | 73 |
#This function generates weights from a point location on a raster layer. |
| 214 | 74 |
#Note that the weights are normatlized on 0-1 scale using max and min values. |
| ... | ... | |
| 216 | 76 |
#lf: list of raster files |
| 217 | 77 |
#df_points: |
| 218 | 78 |
#Outputs: |
| 219 |
# |
|
| 79 |
#TODO: add options to generate weights from edges/reference image rather than reference centroids points...
|
|
| 220 | 80 |
############ |
| 221 | 81 |
|
| 222 | 82 |
lf <- list_param$lf |
| ... | ... | |
| 233 | 93 |
cell_ID <- cellFromXY(r_init,xy=df_centroids[i,]) |
| 234 | 94 |
r_init[cell_ID] <- df_centroids$ID[i] |
| 235 | 95 |
|
| 96 |
#change here to distance from edges... |
|
| 236 | 97 |
r_dist <- distance(r_init) |
| 237 | 98 |
min_val <- cellStats(r_dist,min) |
| 238 | 99 |
max_val <- cellStats(r_dist,max) |
| ... | ... | |
| 314 | 175 |
return(rast_list) |
| 315 | 176 |
} |
| 316 | 177 |
|
| 178 |
raster_match <- function(i,list_param){
|
|
| 179 |
### Read in parameters/arguments |
|
| 180 |
lf_files <- list_param$lf_files |
|
| 181 |
rast_ref <- list_param$rast_ref #name of reference file |
|
| 182 |
file_format <- list_param$file_format #".tif",".rst" or others |
|
| 183 |
out_suffix <- list_param$out_suffix |
|
| 184 |
out_dir_str <- list_param$out_dir_str |
|
| 185 |
|
|
| 186 |
### START SCRIPT ## |
|
| 187 |
|
|
| 188 |
r_m <- raster(rast_ref) #ref image with resolution and extent to match |
|
| 189 |
|
|
| 190 |
set1f <- function(x){rep(NA, x)}
|
|
| 191 |
|
|
| 192 |
inFilename <- lf_files[i] |
|
| 193 |
|
|
| 194 |
extension_str <- extension(inFilename) |
|
| 195 |
raster_name_tmp <- gsub(extension_str,"",basename(inFilename)) |
|
| 196 |
#outFilename <- file.path(out_dir,paste(raster_name_tmp,"_","m_",out_suffix,file_format,sep="")) #for use in function later... |
|
| 197 |
|
|
| 198 |
raster_name <- file.path(out_dir_str,paste(raster_name_tmp,"_","m_",out_suffix,file_format,sep=""))#output file |
|
| 199 |
r_ref <- init(r_m, fun=set1f, filename=raster_name, overwrite=TRUE) |
|
| 200 |
#NAvalue(r_ref) <- -9999 |
|
| 201 |
|
|
| 202 |
cmd_str <- paste("/usr/bin/gdalwarp",inFilename,raster_name,sep=" ")
|
|
| 203 |
#gdalwarp -t_srs '+proj=utm +zone=11 +datum=WGS84' raw_spot.tif utm11.tif |
|
| 204 |
|
|
| 205 |
system(cmd_str) |
|
| 206 |
|
|
| 207 |
##return name of file created |
|
| 208 |
return(raster_name) |
|
| 209 |
} |
|
| 210 |
|
|
| 211 |
|
|
| 317 | 212 |
############################################ |
| 318 | 213 |
#### Parameters and constants |
| 319 | 214 |
|
| 320 |
#on ATLAS |
|
| 215 |
#Data is on ATLAS
|
|
| 321 | 216 |
|
| 322 |
#in_dir <- "~/Data/IPLANT_project/mosaicing_data_test/reg2" |
|
| 323 |
in_dir <- "~/Data/IPLANT_project/mosaicing_data_test/reg1" |
|
| 217 |
in_dir <- "/data/project/layers/commons/NEX_data/mosaicing_data_test" #reg1 is North America and reg5 is Africa |
|
| 218 |
#in_dir <- "/data/project/layers/commons/NEX_data/mosaicing_data_test/reg1" |
|
| 219 |
#in_dir <- "/data/project/layers/commons/NEX_data/mosaicing_data_test/reg2" #Europe |
|
| 220 |
#in_dir <- "/data/project/layers/commons/NEX_data/mosaicing_data_test/reg5" #Africa |
|
| 324 | 221 |
|
| 325 | 222 |
y_var_name <- "dailyTmax" #PARAM1 |
| 326 | 223 |
interpolation_method <- c("gam_CAI") #PARAM2
|
| 327 |
out_suffix <- "mosaic_run10_1500x4500_global_analyses_03252015" #PARAM3
|
|
| 224 |
out_suffix <- "mosaic_run10_1500x4500_global_analyses_05172015" #PARAM3
|
|
| 328 | 225 |
out_dir <- in_dir #PARAM4 |
| 329 | 226 |
create_out_dir_param <- TRUE #PARAM 5 |
| 330 | 227 |
|
| 331 | 228 |
mosaic_plot <- FALSE #PARAM6 |
| 332 | 229 |
|
| 333 | 230 |
#if daily mosaics NULL then mosaicas all days of the year |
| 334 |
day_to_mosaic <- c("20100101","20100102","20100103","20100104","20100105",
|
|
| 335 |
"20100301","20100302","20100303","20100304","20100305", |
|
| 336 |
"20100501","20100502","20100503","20100504","20100505", |
|
| 337 |
"20100701","20100702","20100703","20100704","20100705", |
|
| 338 |
"20100901","20100902","20100903","20100904","20100905", |
|
| 339 |
"20101101","20101102","20101103","20101104","20101105") #PARAM7 |
|
| 231 |
day_to_mosaic <- c("20100831",
|
|
| 232 |
"20100901") #PARAM7 |
|
| 340 | 233 |
|
| 341 | 234 |
#CRS_locs_WGS84 <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84
|
| 342 | 235 |
CRS_WGS84 <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84 #CONSTANT1
|
| 343 | 236 |
CRS_locs_WGS84<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84
|
| 344 | 237 |
|
| 345 | 238 |
proj_str<- CRS_WGS84 #PARAM 8 #check this parameter |
| 346 |
file_format <- ".rst" #PARAM 9
|
|
| 239 |
file_format <- ".tif" #PARAM 9
|
|
| 347 | 240 |
NA_value <- -9999 #PARAM10 |
| 348 | 241 |
NA_flag_val <- NA_value |
| 349 | 242 |
|
| 350 | 243 |
tile_size <- "1500x4500" #PARAM 11 |
| 351 | 244 |
mulitple_region <- TRUE #PARAM 12 |
| 352 | 245 |
|
| 353 |
region_name <- "world" #PARAM 13
|
|
| 246 |
region_name <- "reg1" #PARAM 13 #reg1 is North America
|
|
| 354 | 247 |
plot_region <- FALSE |
| 355 | 248 |
|
| 356 | 249 |
########################## START SCRIPT ############################## |
| 357 | 250 |
|
| 358 | 251 |
|
| 359 |
####### PART 1: Read in data ######## |
|
| 252 |
####### PART 1: Read in data and process data ########
|
|
| 360 | 253 |
|
| 254 |
in_dir <- file.path(in_dir,region_name) |
|
| 255 |
out_dir <- in_dir |
|
| 361 | 256 |
if(create_out_dir_param==TRUE){
|
| 362 | 257 |
out_dir <- create_dir_fun(out_dir,out_suffix) |
| 363 | 258 |
setwd(out_dir) |
| ... | ... | |
| 367 | 262 |
|
| 368 | 263 |
setwd(out_dir) |
| 369 | 264 |
|
| 370 |
###Table 1: Average accuracy metrics |
|
| 371 |
################## PLOTTING WORLD MOSAICS ################ |
|
| 372 | 265 |
|
| 373 |
lf_mosaic_pred_1500x4500 <-list.files(path=file.path(in_dir),
|
|
| 374 |
pattern=paste(".*.tif$",sep=""),full.names=T)
|
|
| 266 |
lf_mosaic <-list.files(path=file.path(in_dir), |
|
| 267 |
pattern=paste(".*.",day_to_mosaic[2],".*.tif$",sep=""),full.names=T) #choosing date 2...20100901
|
|
| 375 | 268 |
|
| 376 |
r1 <- raster(lf_mosaic_pred_1500x4500[1])
|
|
| 377 |
r2 <- raster(lf_mosaic_pred_1500x4500[2])
|
|
| 269 |
r1 <- raster(lf_mosaic[1]) |
|
| 270 |
r2 <- raster(lf_mosaic[2]) |
|
| 378 | 271 |
|
| 379 | 272 |
plot(r1) |
| 380 | 273 |
plot(r2) |
| 381 | 274 |
|
| 382 |
lf <- sub(".tif","",lf_mosaic_pred_1500x4500)
|
|
| 275 |
lf <- sub(".tif","",lf_mosaic)
|
|
| 383 | 276 |
tx<-strsplit(as.character(lf),"_") |
| 384 | 277 |
|
| 385 | 278 |
lat<- as.character(lapply(1:length(tx),function(i,x){x[[i]][13]},x=tx))
|
| 386 | 279 |
long<- as.character(lapply(1:length(tx),function(i,x){x[[i]][14]},x=tx))
|
| 387 | 280 |
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
|
| 388 | 281 |
|
| 282 |
#Produce data.frame with centroids of each tiles... |
|
| 283 |
|
|
| 389 | 284 |
df_centroids <- data.frame(long=as.numeric(long),lat=as.numeric(lat)) |
| 390 | 285 |
df_centroids$ID <- as.numeric(1:nrow(df_centroids)) |
| 391 |
# |
|
| 392 |
#extract(r1,) |
|
| 393 | 286 |
coordinates(df_centroids) <- cbind(df_centroids$long,df_centroids$lat) |
| 394 | 287 |
proj4string(df_centroids) <- projection(r1) |
| 395 |
#centroid distance |
|
| 396 |
#c1 <- gCentroid(x,byid=TRUE) |
|
| 397 |
#pt <- gCentroid(shp1) |
|
| 398 | 288 |
|
| 399 |
#Make this a function later... |
|
| 400 |
#then distance... |
|
| 289 |
############### |
|
| 290 |
### PART 2: prepare weights using tile rasters ############ |
|
| 291 |
|
|
| 401 | 292 |
out_dir_str <- out_dir |
| 402 |
lf_r_weights <- vector("list",length=length(lf_mosaic_pred_1500x4500))
|
|
| 293 |
lf_r_weights <- vector("list",length=length(lf_mosaic))
|
|
| 403 | 294 |
|
| 404 |
list_param_create_weights <- list(lf_mosaic_pred_1500x4500,df_centroids,out_dir_str)
|
|
| 295 |
list_param_create_weights <- list(lf_mosaic,df_centroids,out_dir_str) |
|
| 405 | 296 |
names(list_param_create_weights) <- c("lf","df_points","out_dir_str")
|
| 406 |
num_cores <- 6
|
|
| 297 |
num_cores <- 11
|
|
| 407 | 298 |
|
| 408 | 299 |
#debug(create_weights_fun) |
| 409 | 300 |
#weights_obj <- create_weights_fun(1,list_param=list_param_create_weights) |
| 410 | 301 |
|
| 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) |
|
| 412 |
|
|
| 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" |
|
| 415 |
|
|
| 416 |
#"r_weights","r_weights_prod" |
|
| 417 |
|
|
| 418 |
#list_r_weights <- lapply(1:length(weights_obj_list), FUN=function(i,x){raster(x[[i]]$r_weights)},x=weights_obj_list)
|
|
| 419 |
#list_r_weights_prod <- lapply(1:length(weights_obj_list), FUN=function(i,x){raster(x[[i]]$r_weights_prod)},x=weights_obj_list)
|
|
| 302 |
#This is the function creating the weights by tile. Distance from the centroids needs to be change from distance to |
|
| 303 |
#the edges...can use rows and columsn to set edges to 1 and 0 for the others. |
|
| 304 |
weights_obj_list <- mclapply(1:length(lf_mosaic),FUN=create_weights_fun,list_param=list_param_create_weights,mc.preschedule=FALSE,mc.cores = num_cores) |
|
| 420 | 305 |
|
| 421 | 306 |
list_r_weights <- lapply(1:length(weights_obj_list), FUN=function(i,x){x[[i]]$r_weights},x=weights_obj_list)
|
| 422 | 307 |
list_r_weights_prod <- lapply(1:length(weights_obj_list), FUN=function(i,x){x[[i]]$r_weights_prod},x=weights_obj_list)
|
| 423 | 308 |
|
| 424 |
#r_weights_sum <- do.call(overlay,list_args_weights) #sum |
|
| 425 |
#r_weights_sum <- do.call(overlay,list_args_weights) #sum |
|
| 426 |
|
|
| 427 |
#r_test_w <-do.call(overlay,list_args_w) #prod |
|
| 428 |
|
|
| 429 |
|
|
| 430 |
###Mosaic with do.call... |
|
| 431 |
#rasters1.mosaicargs <- rasters1 |
|
| 432 |
#rasters1.mosaicargs$fun <- mean |
|
| 433 |
#mos2 <- do.call(mosaic, rasters1.mosaicargs) |
|
| 434 |
|
|
| 435 | 309 |
#Then scale on 1 to zero? or 0 to 1000 |
| 436 | 310 |
#e.g. for a specific pixel |
| 437 | 311 |
#weight_sum=0.2 +0.4 +0.4+0.2=1.2 |
| 438 | 312 |
#val_w_sum= (0.2*val1+0.4*val2+0.4*val3+0.2*val4) |
| 439 |
#no_valid = |
|
| 313 |
#no_valid = number of valid observation, needs to be added in the function
|
|
| 440 | 314 |
#m_val= sum_val/weight_sum #mean value |
| 441 | 315 |
# |
| 442 | 316 |
|
| 443 |
#in raster term |
|
| 444 |
#r_weights_sum <- ... |
|
| 445 |
#r_val_w_sum <- |
|
| 446 |
# |
|
| 447 |
#mosaic_list_var <- list(list_r_weights) |
|
| 448 |
#mosaic_list_var <- list(lf_mosaic_pred_1500x4500) |
|
| 449 |
#out_rastnames_var <- l_out_rastnames_var[[i]] |
|
| 450 |
#out_rastnames_var <- c("reg1_mosaic_mean.tif")
|
|
| 451 |
|
|
| 317 |
## Rasters tiles vary slightly in resolution, they need to be matched for the mosaic. Resolve issue in the |
|
| 318 |
#mosaic funciton using gdal_merge to compute a reference image to mach. |
|
| 319 |
#outrastnames <- "reg1_mosaic_weights.tif" |
|
| 452 | 320 |
#list_param_mosaic <- list(list_r_weights,out_dir,outrastnames,file_format,NA_flag_val,out_suffix) |
| 321 |
#r1_projected <- projectRaster(raster(list_r_weights[[1]]),r) |
|
| 453 | 322 |
|
| 454 |
#file_format <- ".tif" |
|
| 455 |
#NA_flag_val <- -9999 |
|
| 456 |
|
|
| 457 |
#j<-1 #date index for loop |
|
| 458 |
#list_param_mosaic<-list(j,mosaic_list_var,out_rastnames_var,out_dir,file_format,NA_flag_val) |
|
| 459 |
#names(list_param_mosaic)<-c("j","mosaic_list","out_rastnames","out_path","file_format","NA_flag_val")
|
|
| 460 |
#undebug(mosaic_m_raster_list) |
|
| 461 |
#r_mosaiced <- mosaic_m_raster_list(1,list_param_mosaic) |
|
| 462 |
|
|
| 463 |
|
|
| 464 |
#list_var_mosaiced <- mclapply(1:2,FUN=mosaic_m_raster_list,list_param=list_param_mosaic,mc.preschedule=FALSE,mc.cores = 2) |
|
| 465 |
#list_var_mosaiced <- mclapply(1,FUN=mosaic_m_raster_list,list_param=list_param_mosaic,mc.preschedule=FALSE,mc.cores = 1) |
|
| 466 |
#list_var_mosaiced <- mclapply(1:1,FUN=mosaic_m_raster_list,list_param=list_param_mosaic,mc.preschedule=FALSE,mc.cores = 1) |
|
| 467 |
#list_var_mosaiced <- mclapply(1:365,FUN=mosaic_m_raster_list,list_param=list_param_mosaic,mc.preschedule=FALSE,mc.cores = 2) |
|
| 468 |
|
|
| 469 |
outrastnames <- "reg1_mosaic_weights.tif" |
|
| 470 |
|
|
| 471 |
list_param_mosaic <- list(list_r_weights,out_dir,outrastnames,file_format,NA_flag_val,out_suffix) |
|
| 472 |
|
|
| 473 |
r1_projected <- projectRaster(raster(list_r_weights[[1]]),r) |
|
| 474 |
|
|
| 475 |
cmd_str <- paste("python","/usr/bin/gdal_merge.py","-o avg.tif",paste(lf_mosaic_pred_1500x4500,collapse=" "))
|
|
| 323 |
cmd_str <- paste("python","/usr/bin/gdal_merge.py","-o avg.tif",paste(lf_mosaic,collapse=" "))
|
|
| 476 | 324 |
system(cmd_str) |
| 477 | 325 |
|
| 478 |
##Create raster image for weights matching resolution and extent to the mosaic |
|
| 479 |
lf_files <- list_r_weights |
|
| 480 |
r_m <- raster("avg.tif")
|
|
| 481 |
|
|
| 482 |
for (i in 1:length(lf_files)){
|
|
| 483 |
set1f <- function(x){rep(NA, x)}
|
|
| 484 |
r_ref <- init(r_m, fun=set1f, filename=paste('r_weights_m_',i,'.tif',sep=""), overwrite=TRUE)
|
|
| 485 |
#NAvalue(r_ref) <- -9999 |
|
| 486 |
|
|
| 487 |
cmd_str <- paste("/usr/bin/gdalwarp",list_r_weights[[i]],paste('r_weights_m_',i,'.tif',sep=""),sep=" ")
|
|
| 488 |
#gdalwarp -t_srs '+proj=utm +zone=11 +datum=WGS84' raw_spot.tif utm11.tif |
|
| 489 |
|
|
| 490 |
system(cmd_str) |
|
| 491 |
#r_ref <-raster("r_ref.tif")
|
|
| 492 |
#plot(r_ref) |
|
| 493 |
} |
|
| 494 |
|
|
| 495 |
##Create raster image for prod weights matching resolution and extent to the mosaic |
|
| 326 |
## Create raster image for original predicted images with matching resolution and extent to the mosaic (reference image) |
|
| 496 | 327 |
|
| 497 |
lf_files <- list_r_weights_prod
|
|
| 498 |
r_m <- raster("avg.tif")
|
|
| 328 |
rast_ref <- file.path(out_dir,"avg.tif")
|
|
| 329 |
lf_files <- unlist(list_r_weights)
|
|
| 499 | 330 |
|
| 500 |
for (i in 1:length(lf_files)){
|
|
| 501 |
set1f <- function(x){rep(NA, x)}
|
|
| 502 |
r_ref <- init(r_m, fun=set1f, filename=paste('r_weights_prod_m_',i,'.tif',sep=""), overwrite=TRUE)
|
|
| 503 |
#NAvalue(r_ref) <- -9999 |
|
| 331 |
list_param_raster_match <- list(lf_files,rast_ref,file_format,out_suffix,out_dir) |
|
| 332 |
names(list_param_raster_match) <- c("lf_files","rast_ref","file_format","out_suffix","out_dir_str")
|
|
| 504 | 333 |
|
| 505 |
cmd_str <- paste("/usr/bin/gdalwarp",lf_files[[i]],paste('r_weights_prod_m_',i,'.tif',sep=""),sep=" ")
|
|
| 506 |
#gdalwarp -t_srs '+proj=utm +zone=11 +datum=WGS84' raw_spot.tif utm11.tif
|
|
| 334 |
#debug(raster_match)
|
|
| 335 |
#r_test <- raster(raster_match(1,list_param_raster_match))
|
|
| 507 | 336 |
|
| 508 |
system(cmd_str) |
|
| 509 |
#r_ref <-raster("r_ref.tif")
|
|
| 510 |
#plot(r_ref) |
|
| 511 |
} |
|
| 337 |
list_weights_m <- mclapply(1:length(lf_files),FUN=raster_match,list_param=list_param_raster_match,mc.preschedule=FALSE,mc.cores = num_cores) |
|
| 512 | 338 |
|
| 513 |
## Create raster image for original predicted images with matching resolution and extent to the mosaic (reference image) |
|
| 514 |
lf_files <- lf_mosaic_pred_1500x4500 |
|
| 515 |
lf_out <- vector("list",length(lf_files))
|
|
| 516 |
r_m <- raster("avg.tif") #ref image
|
|
| 339 |
lf_files <- unlist(list_r_weights_prod) |
|
| 340 |
list_param_raster_match <- list(lf_files,rast_ref,file_format,out_suffix,out_dir) |
|
| 341 |
names(list_param_raster_match) <- c("lf_files","rast_ref","file_format","out_suffix","out_dir_str")
|
|
| 517 | 342 |
|
| 518 |
#Make this a full function... |
|
| 519 |
for (i in 1:length(lf_files)){
|
|
| 520 |
set1f <- function(x){rep(NA, x)}
|
|
| 521 |
|
|
| 522 |
inFilename <- lf_files[i] |
|
| 523 |
extension_str <- extension(inFilename) |
|
| 524 |
raster_name_tmp <- gsub(extension_str,"",basename(inFilename)) |
|
| 525 |
#outFilename <- file.path(out_dir,paste(raster_name_tmp,"_","m_",out_suffix,file_format,sep="")) #for use in function later... |
|
| 526 |
|
|
| 527 |
r_ref <- init(r_m, fun=set1f, filename=paste(raster_name_tmp,"_","m_",out_suffix,file_format,sep=""), overwrite=TRUE) |
|
| 528 |
#NAvalue(r_ref) <- -9999 |
|
| 343 |
num_cores <-11 |
|
| 344 |
list_weights_prod_m <- mclapply(1:length(lf_files),FUN=raster_match,list_param=list_param_raster_match,mc.preschedule=FALSE,mc.cores = num_cores) |
|
| 529 | 345 |
|
| 530 |
cmd_str <- paste("/usr/bin/gdalwarp",inFilename,paste(raster_name_tmp,"_","m_",out_suffix,file_format,sep=""),sep=" ")
|
|
| 531 |
#gdalwarp -t_srs '+proj=utm +zone=11 +datum=WGS84' raw_spot.tif utm11.tif |
|
| 346 |
#macth file to mosaic extent using the original predictions |
|
| 347 |
lf_files <- lf_mosaic |
|
| 348 |
list_param_raster_match <- list(lf_files,rast_ref,file_format,out_suffix,out_dir) |
|
| 349 |
names(list_param_raster_match) <- c("lf_files","rast_ref","file_format","out_suffix","out_dir_str")
|
|
| 532 | 350 |
|
| 533 |
system(cmd_str) |
|
| 534 |
#r_ref <-raster("r_ref.tif")
|
|
| 535 |
#plot(r_ref) |
|
| 536 |
lf_out <- "" |
|
| 537 |
} |
|
| 351 |
list_pred_m <- mclapply(1:length(lf_files),FUN=raster_match,list_param=list_param_raster_match,mc.preschedule=FALSE,mc.cores = num_cores) |
|
| 538 | 352 |
|
| 539 |
list_args_weights_prod <- (mixedsort(list.files(pattern="r_weights_prod_m_.*.tif")))
|
|
| 540 |
list_args_weights_prod <- lapply(1:length(list_args_weights_prod), FUN=function(i,x){raster(x[[i]])},x=list_args_weights_prod)
|
|
| 353 |
#####################
|
|
| 354 |
###### PART 3: compute the weighted mean with the mosaic function #####
|
|
| 541 | 355 |
|
| 542 |
list_args_weights <- (mixedsort(list.files(pattern="r_weights_m_.*.tif"))) |
|
| 356 |
#get the list of weights into raster objects |
|
| 357 |
list_args_weights <- list_weights_m |
|
| 358 |
#list_args_weights <- (mixedsort(list.files(pattern="r_weights_m_.*.tif"))) |
|
| 543 | 359 |
list_args_weights <- lapply(1:length(list_args_weights), FUN=function(i,x){raster(x[[i]])},x=list_args_weights)
|
| 544 | 360 |
|
| 545 |
#lf_mosaic_pred_1500x4500 <-list.files(path=file.path(in_dir), |
|
| 546 |
# pattern=paste(".*.tif$",sep=""),full.names=T)
|
|
| 361 |
#get the list of weights product into raster objects |
|
| 362 |
list_args_weights_prod <- list_weights_prod_m |
|
| 363 |
#list_args_weights_prod <- (mixedsort(list.files(pattern="r_weights_prod_m_.*.tif"))) |
|
| 364 |
list_args_weights_prod <- lapply(1:length(list_args_weights_prod), FUN=function(i,x){raster(x[[i]])},x=list_args_weights_prod)
|
|
| 547 | 365 |
|
| 548 |
list_args_pred_m <- (mixedsort(list.files(pattern="^gam_CAI.*.m_mosaic_run10_1500x4500_global_analyses_03252015.tif"))) |
|
| 366 |
#get the original predicted image to raster (matched previously to the mosaic extent) |
|
| 367 |
list_args_pred_m <- list_pred_m |
|
| 368 |
#list_args_pred_m <- (mixedsort(list.files(pattern="^gam_CAI.*.m_mosaic_run10_1500x4500_global_analyses_03252015.tif"))) |
|
| 549 | 369 |
list_args_pred_m <- lapply(1:length(list_args_pred_m), FUN=function(i,x){raster(x[[i]])},x=list_args_pred_m)
|
| 550 | 370 |
|
| 551 |
list_args_weights$fun <- "sum" |
|
| 552 |
#list_args_weights$fun <- "mean" |
|
| 371 |
list_args_weights$fun <- "sum" #we want the sum to compute the weighted mean |
|
| 553 | 372 |
list_args_weights$na.rm <- TRUE |
| 554 |
#list_args_weights$tolerance <- 1 |
|
| 555 | 373 |
|
| 556 | 374 |
list_args_weights_prod$fun <- "sum" |
| 557 | 375 |
list_args_weights_prod$na.rm <- TRUE |
| ... | ... | |
| 559 | 377 |
list_args_pred_m$fun <- "mean" |
| 560 | 378 |
list_args_pred_m$na.rm <- TRUE |
| 561 | 379 |
|
| 562 |
#l#ist_args_weights_prod$fun <- "sum" |
|
| 563 |
#list_args_weights_prod$na.rm <- TRUE |
|
| 564 |
#list_args_weights_prod$na.rm <- TRUE |
|
| 565 |
|
|
| 566 |
r_weights_sum <- do.call(mosaic,list_args_weights) #sum |
|
| 567 |
r_prod_sum <- do.call(mosaic,list_args_weights_prod) #sum |
|
| 380 |
#Mosaic files |
|
| 381 |
r_weights_sum <- do.call(mosaic,list_args_weights) #weights sum image mosaiced |
|
| 382 |
r_prod_sum <- do.call(mosaic,list_args_weights_prod) #weights sum product image mosacied |
|
| 568 | 383 |
|
| 569 | 384 |
r_m_weighted_mean <- r_prod_sum/r_weights_sum #this is the mosaic using weighted mean... |
| 570 | 385 |
|
| 571 |
#extension_str <- extension(lf[i]) |
|
| 572 |
#raster_name_tmp <- gsub(extension_str,"",basename(lf[i])) |
|
| 573 | 386 |
raster_name <- file.path(out_dir,paste("r_m_weighted_mean",out_suffix,".tif",sep=""))
|
| 574 | 387 |
writeRaster(r_m_weighted_mean, NAflag=NA_flag_val,filename=raster_name,overwrite=TRUE) |
| 575 | 388 |
|
| 576 |
r_m_mean <- do.call(mosaic,list_args_pred_m) #this is unweighted mean |
|
| 389 |
r_m_mean <- do.call(mosaic,list_args_pred_m) #this is unweighted mean from the predicted raster |
|
| 390 |
|
|
| 391 |
raster_name <- file.path(out_dir,paste("r_m_mean",out_suffix,".tif",sep=""))
|
|
| 392 |
writeRaster(r_m_mean, NAflag=NA_flag_val,filename=raster_name,overwrite=TRUE) #unweighted mean |
|
| 577 | 393 |
|
| 578 | 394 |
res_pix <- 1200 |
| 579 | 395 |
col_mfrow <- 1 |
| 580 | 396 |
row_mfrow <- 0.8 |
| 581 | 397 |
|
| 582 |
png(filename=paste("Figure2_mean_for_region_",region_names[1],"_",out_suffix,".png",sep=""),
|
|
| 398 |
png(filename=paste("Figure2_mean_for_region_",region_name,"_",out_suffix,".png",sep=""),
|
|
| 583 | 399 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
| 584 | 400 |
|
| 585 | 401 |
plot(r_m_mean) |
| ... | ... | |
| 590 | 406 |
col_mfrow <- 1 |
| 591 | 407 |
row_mfrow <- 0.8 |
| 592 | 408 |
|
| 593 |
png(filename=paste("Figure2_weigthed_mean_for_region_",region_names[1],"_",out_suffix,".png",sep=""),
|
|
| 409 |
png(filename=paste("Figure2_weigthed_mean_for_region_",region_name,"_",out_suffix,".png",sep=""),
|
|
| 594 | 410 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
| 595 | 411 |
|
| 596 | 412 |
plot(r_m_weighted_mean) |
| ... | ... | |
| 600 | 416 |
|
| 601 | 417 |
##################### END OF SCRIPT ###################### |
| 602 | 418 |
|
| 603 |
|
|
| 604 | 419 |
################################################# |
| 605 |
#Ok testing on fake data: |
|
| 420 |
#Ok testing on fake data to experiment and check methods:
|
|
| 606 | 421 |
|
| 607 | 422 |
##Quick function to generate test dataset |
| 608 | 423 |
# make_raster_from_lf <- function(i,list_lf,r_ref){
|
Also available in: Unified diff
mosaicing test script, major clean up to share code