Add met2CF.ERA5.reanalysis: Convert ERA5 Reanalysis Data to CF-Compliant NetCDF

CHANGELOG.md

@@ -18,6 +18,7 @@ For more information about this file see also [Keep a Changelog](http://keepacha
   - You can compile all pkgdown pages locally at any time with `make pkgdocs`.
 - The ERA5 NC extraction function can now handle multi-site instead of one.
 - Added function for merging images from the same tiling system (MODIS, GLANCE, ICESat-2, HLS, etc.).
+- Added `met2CF.ERA5.reanalysis` function for converting ERA5 reanalysis NetCDF data to CF-compliant format
 
 ### Fixed
 - api to correctly use x_var from request in plotResults #3528

modules/data.atmosphere/NAMESPACE

@@ -72,6 +72,7 @@ export(met2CF.ALMA)
 export(met2CF.Ameriflux)
 export(met2CF.AmerifluxLBL)
 export(met2CF.ERA5)
+export(met2CF.ERA5.reanalysis)
 export(met2CF.FACE)
 export(met2CF.Geostreams)
 export(met2CF.ICOS)

modules/data.atmosphere/NEWS.md

@@ -1,8 +1,14 @@
 # PEcAn.data.atmosphere 1.9.1
 
+## Changed
+
 * `ERA5_met_process()` can now process ensemble data efficiently in parallel using new option `n_cores`
 * Dependency `ggplot2` is now suggested rather than required. It is used in two vignettes and for optional diagnostic plots from `debias_met_regression`.
 
+## Added
+
+* Added `met2CF.ERA5.reanalysis` function for converting ERA5 reanalysis NetCDF data to CF-compliant format
+
 
 # PEcAn.data.atmosphere 1.9.0
 

modules/data.atmosphere/R/met2cf.ERA5.reanalysis.R

@@ -0,0 +1,239 @@
+#' met2CF.ERA5.reanalysis - function for ERA5 reanalysis data processing
+#'
+#' @param slat latitude
+#' @param slon longitude
+#' @param in.path path to the directory containing the ERA5 reanalysis NetCDF files
+#' @param start_date start date
+#' @param end_date end date
+#' @param sitename the name of the site used for making the identifier
+#' @param outfolder path to directory where CF-compliant nc files need to be saved
+#' @param in.prefix initial portion of the filename that does not vary by date 
+#' @param vars variables to be extracted. If NULL all available variables will be returned
+#' @param overwrite logical if files need to be overwritten
+#' @param verbose logical flag defining if output of function be extra verbose
+#'
+#' @return dataframe with file information for processed NetCDF files
+#' @export
+#' @author Akash
+met2CF.ERA5.reanalysis <- function(slat, slon, in.path, start_date, end_date, sitename, outfolder,
+                                   in.prefix, vars = NULL, overwrite = FALSE, verbose = TRUE) {
+
+  if (is.na(slat) || is.na(slon)) {
+    PEcAn.logger::logger.severe("Invalid latitude or longitude provided")
+    return(NULL)
+  }
+  years <- seq(lubridate::year(start_date),
+               lubridate::year(end_date),
+               1)
+
+  # define variable mapping from ERA5 to CF standard names
+  cf_mapping <- c(
+    "t2m" = "air_temperature",
+    "sp" = "air_pressure", 
+    "tp" = "precipitation_flux",
+    "u10" = "eastward_wind",
+    "v10" = "northward_wind",
+    "ssrd" = "surface_downwelling_shortwave_flux_in_air",
+    "strd" = "surface_downwelling_longwave_flux_in_air"
+  )
+
+  # define units mapping for CF variables
+  cf_units_mapping <- c(
+    "air_temperature" = "K",
+    "air_pressure" = "Pa", 
+    "precipitation_flux" = "kg m-2 s-1",
+    "eastward_wind" = "m s-1",
+    "northward_wind" = "m s-1",
+    "surface_downwelling_shortwave_flux_in_air" = "W m-2",
+    "surface_downwelling_longwave_flux_in_air" = "W m-2",
+    "specific_humidity" = "1"
+  )
+
+  tryCatch({
+    # process each year
+    out.xts <- years %>%
+      purrr::map(function(year) {
+        ncfile <- file.path(in.path, paste0(in.prefix, year, ".nc"))
+        if (!file.exists(ncfile)) {
+          PEcAn.logger::logger.warn(paste0("File not found: ", ncfile))
+          return(NULL)
+        }
+        nc_data <- ncdf4::nc_open(ncfile)
+        on.exit(ncdf4::nc_close(nc_data), add = TRUE)
+
+        time_var <- if ("time" %in% names(nc_data$var)) "time" else "valid_time"
+        t <- ncdf4::ncvar_get(nc_data, time_var)
+        tunits <- ncdf4::ncatt_get(nc_data, time_var)
+        tustr <- strsplit(tunits$units, " ")
+
+        if(time_var == "time") {
+          # legacy format - "hours since YYYY-MM-DD HH:MM:SS"
+          timestamp <- as.POSIXct(t * 3600, tz = "UTC", origin = tustr[[1]][3])
+        } else {
+          # new format - "seconds since YYYY-MM-DD HH:MM:SS"
+          timestamp <- as.POSIXct(t, tz = "UTC", origin = tustr[[1]][3])
+        }
+
+        lon <- ncdf4::ncvar_get(nc_data, "longitude")
+        lat <- ncdf4::ncvar_get(nc_data, "latitude")
+        lon_ix <- which.min(abs(lon - slon))
+        lat_iy <- which.min(abs(lat - slat))
+
+        if(is.null(vars)) {
+          all_vars <- names(nc_data$var)
+          vars <- all_vars[sapply(all_vars, function(v) {
+            var_info <- nc_data$var[[v]]
+            var_info$ndims == 3 && 
+              var_info$prec %in% c("float", "double", "integer") &&
+              !v %in% c("longitude", "latitude", "time", "valid_time")
+          })]
+        } else {
+          vars <- vars
+        }
+
+        if(length(vars) == 0) {
+          PEcAn.logger::logger.warn(paste0("No valid variables found in ", ncfile))
+          return(NULL)
+        }
+
+        extracted_data <- vars %>%
+          purrr::set_names(vars) %>%
+          purrr::map_dfc(function(vname) {
+            vals <- ncdf4::ncvar_get(nc_data, vname)
+            if(!is.null(nc_data$var[[vname]]$misval)){
+              missval <- nc_data$var[[vname]]$missval
+            } else {
+              missval <- NA
+            }
+            if (length(dim(vals)) == 3) {
+              point_vals <- vals[lon_ix, lat_iy, ]
+            } else {
+              PEcAn.logger::logger.warn(paste0("Unexpected dimensions for variable ", vname))
+              point_vals <- rep(NA, length(timestamp))
+            }
+            point_vals[point_vals == missval] <- NA
+            return(point_vals)
+          })
+        # create xts object
+        ens <- xts::xts(extracted_data, order.by = timestamp)
+        return(ens)
+      }) %>%
+      purrr::discard(is.null)
+    if(length(out.xts) == 0) {
+      PEcAn.logger::logger.severe("No data successfully extracted from any files")
+      return(NULL)
+    }
+    rbind.xts <- do.call("::", list("xts", "rbind.xts"))
+    out.new <- do.call("rbind.xts", out.xts)
+    vars <- colnames(out.new)
+
+    time_diffs <- diff(as.numeric(zoo::index(out.new)))
+    if(length(time_diffs) > 0) {
+      timestep_seconds <- as.numeric(median(time_diffs))
+    } else {
+      timestep_seconds <- 3600 
+    }
+    timestep_hours <- timestep_seconds / 3600
+
+    # solar radiation conversions(J/m2 to W/m2)
+    if("ssrd" %in% vars) {
+      out.new[, "ssrd"] <- out.new[, "ssrd"] / timestep_seconds
+    }
+    if("strd" %in% vars) {
+      out.new[, "strd"] <- out.new[, "strd"] / timestep_seconds
+    }
+    # precipitation conversion(m to kg/m2/s)
+    if("tp" %in% vars) {
+      out.new[, "tp"] <- (out.new[, "tp"] * 1000) / timestep_seconds
+    }
+    specific_humidity <- NULL
+    if(all(c("t2m", "d2m", "sp") %in% vars)) {
+      t <- PEcAn.utils::ud_convert(out.new[, "t2m"] %>% as.numeric(), "K", "degC")
+      dewpoint <- PEcAn.utils::ud_convert(out.new[, "d2m"] %>% as.numeric(), "K", "degC")
+      beta <- (112 - (0.1 * t) + dewpoint) / (112 + (0.9 * t))
+      relative.humidity <- beta ^ 8
+
+      specific_humidity <- PEcAn.data.atmosphere::rh2qair(
+        relative.humidity,
+        out.new[, "t2m"] %>% as.numeric(),
+        out.new[, "sp"] %>% as.numeric()
+      )
+    }
+    era5_vars <- intersect(names(cf_mapping), vars)
+    if(length(era5_vars) == 0) {
+      PEcAn.logger::logger.severe("No mappable ERA5 variables found")
+      return(NULL)
+    }
+    cf_data <- out.new[, era5_vars, drop = FALSE]
+    colnames(cf_data) <- cf_mapping[era5_vars]
+    if(!is.null(specific_humidity)) {
+      colnames(specific_humidity) <- "specific_humidity"
+      cf_data <- xts::merge.xts(cf_data, specific_humidity)
+    }
+    cf_var_names <- colnames(cf_data)
+    cf_var_units <- cf_units_mapping[cf_var_names]
+    identifier <- paste("ERA5", sitename, "reanalysis", sep = "_")
+    reanalysis_folder <- file.path(outfolder, identifier)
+
+    if(!dir.exists(reanalysis_folder)) {
+      dir.create(reanalysis_folder, recursive = TRUE, showWarnings = FALSE)
+    }
+
+    start_date <- min(zoo::index(cf_data))
+    end_date <- max(zoo::index(cf_data))
+
+    results <- data.frame(
+      file = "",
+      host = PEcAn.remote::fqdn(),
+      mimetype = "application/x-netcdf",
+      formatname = "CF Meteorology",
+      startdate = paste0(format(start_date, "%Y-%m-%dT%H:%M:00 %z")),
+      enddate = paste0(format(end_date, "%Y-%m-%dT%H:%M:00 %z")),
+      dbfile.name = "ERA5.reanalysis",
+      stringsAsFactors = FALSE
+    )
+
+    years %>%
+      purrr::walk(function(year) {
+        year_data <- cf_data[year %>% as.character]
+        if(nrow(year_data) == 0) {
+          return(NULL)
+        }
+        identifier_file <- paste("ERA5", "reanalysis", year, sep = ".")
+        flname <- file.path(reanalysis_folder, paste(identifier_file, "nc", sep = "."))
+        if(!file.exists(flname) || overwrite) {
+          tryCatch({
+            time_vals <- as.numeric(zoo::index(year_data))
+            time_dim <- ncdf4::ncdim_def(
+              name = "time",
+              units = "seconds since 1970-01-01 00:00:00",
+              vals = time_vals,
+              create_dimvar = TRUE
+            )
+            lat_dim <- ncdf4::ncdim_def("latitude", "degree_north", slat, create_dimvar = TRUE)
+            lon_dim <- ncdf4::ncdim_def("longitude", "degree_east", slon, create_dimvar = TRUE)
+            nc_vars <- purrr::map2(cf_var_names, cf_var_units,
+                                    ~ ncdf4::ncvar_def(.x, .y, list(time_dim, lat_dim, lon_dim), 
+                                                      missval = NA_real_))
+            nc_flptr <- ncdf4::nc_create(flname, nc_vars, verbose = FALSE)
+
+            for(j in seq_along(cf_var_names)) {
+              ncdf4::ncvar_put(nc_flptr, nc_vars[[j]], 
+                               zoo::coredata(year_data)[, j])
+            }
+            ncdf4::nc_close(nc_flptr)
+          }, error = function(e) {
+            PEcAn.logger::logger.severe(paste0("Error writing NetCDF file for year ", year, ": ", conditionMessage(e)))
+          })
+        }
+      })
+
+    year <- years[1]
+    identifier.file <- paste("ERA5", "reanalysis", year, sep = ".")
+    results$file <- file.path(reanalysis_folder, paste(identifier.file, "nc", sep = "."))
+    return(results)
+  }, error = function(e) {
+    PEcAn.logger::logger.severe(paste0("Error in met2CF.ERA5.rea: ", conditionMessage(e)))
+    return(NULL)
+  })
+}