Homec4science
Diffusion minerve-d3 (master) Browse

R10424/Code/fct_get_data_pixel_from_raster.R
db97b68e7e43master

Last Change
Raw File

fct_get_data_pixel_from_raster.R
View Options

get_data_pixel_from_raster <- function(hours_to_consider,product_j,sensor_var,category_var,unit_var,interpolation_var,sp_res,suffix,extension,request_base) {
# This function contains 3 parts
# 1) Downloading of JSON files
# 2) Read all the downloaded files and create a 3D matrix of all the data
# 3) Write files corresponding to time series of each pixel
period_start <- strsplit(as.character(config$compute_time_down),"T")[[1]][1]
period_end <- strsplit(as.character(config$compute_time_up),"T")[[1]][1]
# Log
log_file <- file.path(config$dir$DATA,paste0(period_start,"TO",period_end),'Logs',paste0('log_',product_j,'_',config$catchment_ID,'.txt'))
dir.create(file.path(config$dir$DATA,paste0(period_start,"TO",period_end),'Logs'), showWarnings = FALSE, recursive = TRUE)
if (!file.exists(log_file)){
sink(file = log_file, type = "output", append = TRUE)
print.noquote('Start acquiring RAW data.')
# Time stamps definition
seq_posixct <- seq(as.POSIXct(config$compute_time_down, "%Y-%m-%dT%H:%M:%SZ",tz= 'UTC')-extension*3600,
as.POSIXct(config$compute_time_up, "%Y-%m-%dT%H:%M:%SZ",tz= 'UTC')+extension*3600,"hours")
pos_to_consider <- as.POSIXlt(seq_posixct)$hour %in% hours_to_consider
seq_posixct <- seq_posixct[pos_to_consider]
time_stamp_list <- format(seq_posixct,"%Y-%m-%dT%H:%M:%SZ")
request_txt <- gsub("productId",product_j,request_base)
# Preparation of the request
if (grepl("list_of_coords",request_base)) {
txt_coords_list <- c()
for (j in seq(config$list_coords)){
# Switch coordinates to upper left corner dependig on spatial resolution of the product
coords_j <- upper_left_coordinate(c(config$list_coords[[j]][1],config$list_coords[[j]][2]),sp_res)
txt_coords_list <- paste0(txt_coords_list,coords_j[1],',',coords_j[2],';')
}
# Remove the last semi-column char
config$txt_coords_list <- gsub('.{1}$', '', txt_coords_list)
# Request for list of coordinates
request_txt <- gsub("list_of_coords",config$txt_coords_list,request_txt)
sp_res <- 0 # This is to avoid shift that are required in the data preparation for Snow_Analyzer than are not desirable fr SWE_Analyzer
} else {
# Request for pixels contained in a catchment
request_txt <- gsub("catchmentId",config$catchment_ID,request_txt)
}
dir.create(file.path(config$dir$DATA,'RAW',product_j,config$catchment_ID), showWarnings = FALSE, recursive = TRUE)
# Download the RAW files
worker_acquisition <- function(inc) {
# inc <- 1
# Prepare the request for the time stamp
time_stamp <- time_stamp_list[inc]
file_name <- file.path(config$dir$DATA,'RAW',product_j,config$catchment_ID,paste0(product_j,'_',config$catchment_ID,'_',gsub(':00:00Z','',time_stamp),'.json'))
print(file_name)
# Download the file only if it does not exist
if (!file.exists(file_name)) {
tryCatch({
# Download the data
request_txt <- gsub("ENDDATE",time_stamp,request_txt)
print(request_txt)
header_opts <- paste("--header \"Authorization: Bearer",token,"\"")
r <- GET(request_txt, add_headers(Authorization = paste("Bearer",token)), write_disk(file_name, overwrite = TRUE))
}, error = function(err) {
print(err)
})
}
# sink()
}
cl <- makeCluster(max(1,detectCores()-config$cores2exclude_acquisition))
clusterEvalQ(cl,library(jsonlite))
clusterEvalQ(cl,library(httr))
clusterExport(cl,c("config","request_txt","time_stamp_list","product_j","token"), envir = environment())
worker_status <- parLapply(cl,1:length(time_stamp_list),worker_acquisition)
stopCluster(cl)
# Wait for file to be written
Sys.sleep(2)
# Remove empty RAW downloaded files
worker_acquisition <- function(inc) {
sink(file = log_file, type = "output", append = TRUE)
# inc <- 1
# Prepare the request for the time stamp
time_stamp <- time_stamp_list[inc]
file_name <- file.path(config$dir$DATA,'RAW',product_j,config$catchment_ID,paste0(product_j,'_',config$catchment_ID,'_',gsub(':00:00Z','',time_stamp),'.json'))
# Check length of all files in the folder and remove all the empty files.
if (file.exists(file_name)) {
tryCatch({
print.noquote(paste0("Removing empty file ",file_name))
if (file.info(file_name)$size == 0){
file.remove(file_name)
}
}, error = function(err) {
print(err)
})
}
sink()
}
cl <- makeCluster(max(1,detectCores()-config$cores2exclude_acquisition))
clusterExport(cl,c("config","time_stamp_list","product_j","log_file"), envir = environment())
worker_status <- parLapply(cl,1:length(time_stamp_list),worker_acquisition)
stopCluster(cl)
print.noquote('Finished acquiring RAW data.')
#########################################
# LOOP (i) OVER ALL TIME STAMPS TO READ THE RAW FILES
#########################################
print.noquote('Start loading RAW data.')
# List of downloaded files
FileList <- list.files(file.path(config$dir$DATA,'Raw',product_j,config$catchment_ID))
file.size(FileList)
# Initialize
time_stamp <- c()
for (i in seq(FileList)) {
# Read the files
FileData <- jsonlite::fromJSON((file.path(file.path(config$dir$DATA,'Raw',product_j,config$catchment_ID),FileList[i])))
if (i == 1){
coords_list_base <- FileData$pixels[,1:2]
# Initialize
product_data <- matrix(data = NA, nrow = length(coords_list_base[[1]]), ncol = 1, byrow = FALSE,dimnames = NULL)
} else {
if (length(suppressMessages(anti_join(coords_list_base, FileData$pixels[,1:2]))$e) != 0) {
print.noquote("The list of coordinates of this time stamp does not correspond to the list of the first time stamp.")
print.noquote("An error will be generated to stop the code.")
stop()
}
# Correction for precipitation data
#Only if we work with precipitation data (and starting at second time step)
if (category_var=='Precipitation'){
time_past <- as.POSIXct(time_stamp[length(time_stamp)], "%Y-%m-%dT%H:%M:%SZ",tz= 'UTC')
time_current <- as.POSIXct(gsub('.000Z','Z',FileData$event), "%Y-%m-%dT%H:%M:%SZ",tz= 'UTC')
while (as.numeric(difftime(time_current,time_past, units = "mins")) != 60){
product_data <- cbind(product_data,rep(0, length(coords_list_base[[1]])))
time_stamp <- append(time_stamp,format(time_past + 3600,"%Y-%m-%dT%H:%M:00Z"))
time_past <- time_past + 3600
print.noquote(paste0("0 mm/h considered for missing time step ",format(time_past,"%Y-%m-%dT%H:%M:00Z")))
}
}
}
if (i == 1){
product_data[,i] <- FileData$pixels[,3]
time_stamp[i] <- gsub('.000Z','Z',FileData$event)
}else{
product_data <- cbind(product_data,FileData$pixels[,3])
time_stamp <- append(time_stamp,gsub('.000Z','Z',FileData$event))
}
}
# Bring to Upper left corner centroids of pixels
if (grepl("catchmentId",request_base)) {
coords_list_base[,1] <- coords_list_base[,1] - sp_res/2
coords_list_base[,2] <- coords_list_base[,2] + sp_res/2
}
print.noquote('Finished loading RAW data.')
#########################################
# LOOP (i) OVER ALL PIXELS TO CREATE THE FILES
#########################################
# i <- 1
token <- get_token_server(config)
print.noquote('Start creating pixel data files.')
sensor_line <- c('Sensor', sensor_var)
category_line <- c('Category',category_var)
unit_line <- c('Unit',unit_var)
interpolation_line <- c('Interpolation', interpolation_var)
n_coordinates <- length(coords_list_base[[1]])
for (i in 1 :n_coordinates) {
# Filepath preparation
generic_filepath <- file.path(config$dir$DATA, paste0(period_start,"TO",period_end),sensor_var,product_j,
paste0(coords_list_base[i,1],'_',coords_list_base[i,2]),paste0(category_var,suffix, '.csv'))
if (!file.exists(generic_filepath)) {
### HEADER
lines <- c()
stations_names <- c('Station')
X <- c('X')
Y <- c('Y')
Z <- c('Z')
# Load Elevation data
dem_prod_id <- 3
dir.create(file.path(config$dir$ELEVATION,dem_prod_id),showWarnings = FALSE, recursive = TRUE)
output_file <- file.path(config$dir$ELEVATION,dem_prod_id,paste0('DEM_',coords_list_base[i,1],'_',coords_list_base[i,2],'.RData'))
if (file.exists(output_file)) {
load(output_file)
} else {
altitude_to_compute_k <- tryCatch({
altitude_to_compute_k <- download_altitude_of_pixel(config,token,c(coords_list_base[i,1],coords_list_base[i,2]))[[1]]$V1
},
error = function(e) {
print.noquote('No elevation value for this pixel !!')
altitude_to_compute_k <- NA
return(altitude_to_compute_k)
}
)
dir.create(file.path(config$dir$ELEVATION,dem_prod_id),showWarnings = FALSE, recursive = TRUE)
save(altitude_to_compute_k , file = file.path(config$dir$ELEVATION,dem_prod_id,paste0('DEM_',coords_list_base[i,1],'_',coords_list_base[i,2],'.RData')))
}
if (sensor_var == 'H_Snow'){
stations_names <- c(stations_names, 'Snow')
}else{
stations_names <- c(stations_names, paste0(coords_list_base[i,1],'_',coords_list_base[i,2]))
}
X <- c(X, coords_list_base[i,1])
Y <- c(Y, coords_list_base[i,2])
Z <- c(Z, altitude_to_compute_k)
### DATA
# product_data[i,][product_data[i,]>1000] <- 0 # To be reactivated in case of corrupeted files to be corrected in the database. (cf r1798 and r1814)
data_table <- cbind(format(as.POSIXct(time_stamp, "%Y-%m-%dT%H:%M:%SZ",tz= 'UTC'),"%d.%m.%Y %H:%M"), product_data[i,])
header_content <- rbind(stations_names, X, Y, Z, sensor_line, category_line, unit_line, interpolation_line)
## merge Header and data table
table_csv <- rbind(header_content, data_table)
# Write the file
dir.create(file.path(config$dir$DATA, paste0(period_start,"TO",period_end),sensor_var,product_j,
paste0(coords_list_base[i,1],'_',coords_list_base[i,2])),
showWarnings = FALSE, recursive = TRUE)
write.table(x = table_csv, file = generic_filepath, quote = F, sep = ';', row.names = F, col.names = F)
}
}
print.noquote('Finished creating pixel data files.')
}else{
sink(file = log_file, type = "output", append = TRUE)
print.noquote('Log file already exists. Nothing done.')
}
sink()
}
c4science · Help