0

I am trying to process a large number of HRRR files, about six months' worth of them, which is close to 4059 files. I am working in R, and this is the function I have, that I am trying to vectorize over all the files in list.files('path to HRRR files, full.names = T):

    hrrr_processing_fn <- function(hrrr_file,vars){
  file_date <- str_extract(tail(unlist(str_split(hrrr_file,"/")),1),"[^\\D$]{8}\\.[^\\D$]{2}")
  
  message(hrrr_file)
  
  # 1) Identify the layer number of each of the variables
  index.file <- varfiles[grepl(file_date,varfiles)]
  
  hrrr.idx <- fread(index.file, sep = ":", header = F)
  hrrr.idx[,V5 := str_replace_all(V5, " ","_")]
  var_num <- hrrr.idx[V5 %in% c("surface", "2_m_above_ground") & V4 %in% vars,V1]
  colnames <- hrrr.idx[var_num,paste0(V4,"_",V5)]
  
  # 2) Read the grib file as a RasterStack
  hrrr_stack <- stack(hrrr_file, bands = var_num)
  
  # 3) Reproject the rasterstack to the required CRS
  hrrr_stack_crs <- projectRaster(hrrr_stack, crs = crs_latlon$proj4string, method = "ngb")
  
  # 5) Mask the rasterstack to the shapefile
  hrrr_stack_masked <- mask(hrrr_stack_crs, mask = suffolk_county)
  
  # 6) Extracts the centroid and data of the masked raster as a data.table
  hrrr_mask_dt <- as.data.frame(hrrr_stack_masked,na.rm = T,xy = T)
  setDT(hrrr_mask_dt)
  setnames(hrrr_mask_dt,c("lon","lat",colnames))
  file_date <- str_replace(file_date, "\\.","_")
  hrrr_mask_dt[,date := ymd_h(file_date, tz = "America/New_York")]
  
}

where crs_latlon = Coordinate Reference System: EPSG: 4326 proj4string: "+proj=longlat +datum=WGS84 +no_defs"

and vars = 'TMP'

While most files are being processed, there are some files which throw an error message, such as this: Error in .local(.Object, ...) : /home/my_username/path_to_files/20160806.22.grib2 is a grib file, but no raster dataset was successfully identified. This is not the only error message I get, though. Some files are just not readable. I have ensured that the files are present and downloaded correctly, because I got a '200' message when I downloaded the file using download.file().

These errors halt the function from running. Given the large number of files, I do not want to stop at each bad file and start over again. I want to write a way into my function to not halt the processing, but log the files that are 'bad' in a different file, and keep processing the other files. One way I could think of was to use a conditional statement examining if the file was readable and read without problems, in which case to go ahead with the processing, or if the file could not be read, which would lead to logging the file name, but I am unsure how to operationalise it. I have tried the assertive package function is_executable_file, which does recognise the grib2 file in question to be an executable file, but the file is nevertheless not read by raster::stack. I have a vague idea that I need to use tryCatch for this, but my attempts so far have failed.

2
  • Take a look at try and tryCatch there is quite a bit of information out there on how to use these function for exactly what you describe. Commented May 19, 2020 at 0:28
  • try projecting the polygons to the raster instead of warping the grid to the polygons i.e. ` hrrr_stack_masked <- mask(hrrr_stack, mask = spTransform(suffolk_county, projection(hrrr_stack)))` it's virtually lossless compared to the drastic remodelling of a grid (mystifies me why this isn't commonly said given how fundamental here ...)
    – mdsumner
    Commented May 19, 2020 at 8:31

1 Answer 1

0

Thanks a lot @JeffreyEvans and @mdsumner for your comments. I was able to rework my function the following way to enable it to vectorize over all the files:

hrrr_processing_fn <- function(hrrr_file,vars){
  file_date <- str_extract(tail(unlist(str_split(hrrr_file,"/")),1),"[^\\D$]{8}\\.[^\\D$]{2}")

  message(hrrr_file)

  # 1) Identify the layer number of each of the variables
  index.file <- varfiles[grepl(file_date,varfiles)]

  hrrr.idx <- fread(index.file, sep = ":", header = F)
  hrrr.idx[,V5 := str_replace_all(V5, " ","_")]
  var_num <- hrrr.idx[V5 %in% c("surface", "2_m_above_ground") & V4 %in% vars,V1]
  colnames <- hrrr.idx[var_num,paste0(V4,"_",V5)]

  # 2) Read the grib file as a RasterStack
  hrrr_stack <- stack(hrrr_file, bands = var_num)

  # 3) Mask the rasterstack to the shapefile
  hrrr_stack_masked <- mask(hrrr_stack, mask = st_transform(suffolk_county,crs = crs(hrrr_stack)))

  # 4) Extracts the centroid and data of the masked raster as a data.table
  hrrr_mask_dt <- as.data.frame(hrrr_stack_masked,na.rm = T,xy = T)
  setDT(hrrr_mask_dt)
  setnames(hrrr_mask_dt,c("lon","lat",colnames))
  hrrr_mask_dt[,ID := .I]
  file_date <- str_replace(file_date, "\\.","_")
  hrrr_mask_dt[,date := ymd_h(file_date, tz = "America/New_York")]

}

hrrr_2017_dt <- lapply(grib.2017,function(hrrr_file, vars = "TMP"){
  tryCatch(hrrr_processing_fn(hrrr_file, vars = "TMP"),
           error = function(err){
             message( hrrr_file ," could not be processed. Error: ",err)
             write(hrrr_file,file = file.path(hrrr_path,"hrrr_error_files.txt"), 
                   append = TRUE)})})

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.