diff --git a/book_source/03_topical_pages/06_data/01_meteorology.Rmd b/book_source/03_topical_pages/06_data/01_meteorology.Rmd index e37b8fe39eb..d55a1ed2142 100644 --- a/book_source/03_topical_pages/06_data/01_meteorology.Rmd +++ b/book_source/03_topical_pages/06_data/01_meteorology.Rmd @@ -132,3 +132,13 @@ Resolution: 30 min Availability: Varies by [site](https://meta.icos-cp.eu/collections/q4V7P1VLZevIrnlsW6SJO1Rz) Notes: To use this option, set `source` as `ICOS` and a `product` tag containing `etc` in `pecan.xml` + +## ECMWF Open Data + +Scale: Global + +Resolution: 6 hour, 0.4 degree + +Availability: Varies by [day](https://data.ecmwf.int/) + +Notes: 15 day forecast is provided at two hours 00 and 12, but current implementation only uses forecast hour 00 when download function is called from within the PEcAn workflow diff --git a/modules/data.atmosphere/DESCRIPTION b/modules/data.atmosphere/DESCRIPTION index a29287049c5..3a3b61ecccb 100644 --- a/modules/data.atmosphere/DESCRIPTION +++ b/modules/data.atmosphere/DESCRIPTION @@ -49,6 +49,7 @@ Imports: raster, REddyProc, reshape2, + reticulate, rgdal, rlang (>= 0.2.0), sp, @@ -67,7 +68,6 @@ Suggests: foreach, parallel, progress, - reticulate Remotes: github::chuhousen/amerifluxr, github::ropensci/geonames, diff --git a/modules/data.atmosphere/NAMESPACE b/modules/data.atmosphere/NAMESPACE index 7183782887f..f94aaf01808 100644 --- a/modules/data.atmosphere/NAMESPACE +++ b/modules/data.atmosphere/NAMESPACE @@ -20,6 +20,7 @@ export(debias.met.regression) export(download.Ameriflux) export(download.AmerifluxLBL) export(download.CRUNCEP) +export(download.ECMWF) export(download.ERA5.old) export(download.FACE) export(download.Fluxnet2015) diff --git a/modules/data.atmosphere/R/ECMWF_helper_functions.R b/modules/data.atmosphere/R/ECMWF_helper_functions.R new file mode 100644 index 00000000000..6bc1b58a994 --- /dev/null +++ b/modules/data.atmosphere/R/ECMWF_helper_functions.R @@ -0,0 +1,290 @@ +#' Helper functions for met2CF.ECMWF +#' +#' mergenc_ECMWF writes 1 CF and 50 PF netCDF files by extracting data from GRIB2 files and filling NaNs at 3h position wherever data is unavailable. +#' +#' @param lat latitude +#' @param lon longitude +#' @param outfolder Path to directory where nc files need to be saved. +#' @param overwrite Logical if files needs to be overwritten. +#' @param verbose Logical flag defining if output of function be extra verbose. +#' +#' @return list of dataframes +#' @export +#' +#' \dontrun{ +#' results <- mergenc_ECMWF(lat.in, +#' lon.in, +#' outfolder, +#' overwrite = overwrite) +#' +#' @author Swarnalee Mazumder +#' +nanfill <- function(variable) { # nanfill function fills NaN at every 3h period from 144h to 360h (original step 6h) + + len6 <- variable + len3nan <- rep(NA,length(len6)*2-1) + len3nan[seq(1, length(len6)*2-1, 2)] <- len6 + + return(len3nan) +} + +combine2nc_ECMWF <- function(lat.in, + lon.in, + outfolder, + overwrite = FALSE, + verbose = TRUE) { + + if (!file.exists(outfolder)) { + dir.create(outfolder, showWarnings = FALSE, recursive = TRUE) + } + + + # Python script "download_ecmwf.py" to get latest forecast date and download forecast datasets + script.path = file.path(system.file("ECMWF/download_ecmwf.py", package = "PEcAn.data.atmosphere")) + reticulate::source_python(script.path) + + time <- 0 + stream <- "enfo" + type <- c("cf", "pf") + step_15day <- 360 + + all_parameters <- c("10u", "10v", "2t", "sp", "tp") + + date_latestdata <- ecmwflatest(time= time, step15= step_15day, stream= stream, type= type, params= all_parameters) + + latest_filedate <- strtoi(gsub("-", "", substr(as.character.Date(date_latestdata), 1, 10))) + + current_date <- strtoi(gsub("-", "", Sys.Date())) + + in_fname <- paste(latest_filedate, time, stream, sep = "_") + + script.path = file.path(system.file("ECMWF/met2CFutils.py", package = "PEcAn.data.atmosphere")) + reticulate::source_python(script.path) + + fnames_3h <- all_filenames_3h(in_fname) + fnames_6h <- all_filenames_6h(in_fname) + + ######### CF ######### + + ##### 3h CF -> 0 - 141h at 3h timestep + cf_3h <- combine_files_cf(fnames_3h, lat_in, lon_in) + names(cf_3h) <- c("valid_time", "lon", "lat", "u10", "v10", "t2m", "sp", "tp") + + days_since <- paste("days since", cf_3h$valid_time$time$data, "00:00:00") + cf_3h$valid_time <- cf_3h$valid_time$values + + ##### 3h CF -> 144h - 360h at 6h timestep + cf_6h_nc <- combine_files_cf(fnames_6h, lat_in, lon_in) + names(cf_6h) <- c("valid_time", "lon", "lat", "u10", "v10", "t2m", "sp", "tp") + + cf_6h$valid_time <- cf_6h$valid_time$values + + ### Final CF netCDF filename + cf_nc_fname <- paste(latest_filedate, time, stream, "cf.nc", sep = "_") + + # values from cf_3h + lat <- cf_3h$lat + lon <- cf_3h$lon + validtime <- seq(0, 360, 3) + validtimeunits <- cf_3h$valid_time + + # create and write the netCDF file -- ncdf4 version + # define dimensions + londim <- ncdf4::ncdim_def("longitude","degrees_east",as.double(lon)) + latdim <- ncdf4::ncdim_def("latitude","degrees_north",as.double(lat)) + timedim <- ncdf4::ncdim_def("valid_time",validtimeunits,as.double(validtime)) + time <- ncdf4::ncdim_def("time", days_since, as.double(validtime)) + + # define variables + fillvalue <- NA + + cf_u10_3 <- cf_3h$u10 + cf_u10_6_nan <- nanfill(cf_6h$u10) + cf_u10 <- c(cf_u10_3, cf_u10_6_nan) + + cf_v10_3 <- cf_3h$v10 + cf_v10_6_nan <- nanfill(cf_6h$v10) + cf_v10 <- c(cf_v10_3, cf_v10_6_nan) + + cf_t2m_3 <- cf_3h$t2m + cf_t2m_6_nan <- nanfill(cf_6h$t2m) + cf_t2m <- c(cf_t2m_3, cf_t2m_6_nan) + + cf_sp_3 <- cf_3h$sp + cf_sp_6_nan <- nanfill(cf_6h$tp) + cf_sp <- c(cf_sp_3, cf_sp_6_nan) + + cf_tp_3 <- cf_3h$tp + cf_tp_6_nan <- nanfill(cf_6h$tp) + cf_tp <- c(cf_tp_3, cf_tp_6_nan) + + dlname <- "Eastward Component of Wind" + u10.def <- ncdf4::ncvar_def("eastward_wind","m s**-1",list(timedim),fillvalue,dlname,prec="double") + + dlname <- "Northward Component of Wind" + v10.def <- ncdf4::ncvar_def("northward_wind","m s**-1",list(timedim),fillvalue,dlname,prec="double") + + dlname <- "Near surface air temperature" + t2m.def <- ncdf4::ncvar_def("air_temperature","K",list(timedim),fillvalue,dlname,prec="double") + + dlname <- "Surface pressure" + sp.def <- ncdf4::ncvar_def("air_pressure","Pa",list(timedim),fillvalue,dlname,prec="double") + + dlname <- "Rainfall rate" + tp.def <- ncdf4::ncvar_def("precipitation_flux","kg m**-2 s**-1",list(timedim),fillvalue,dlname,prec="double") + + if (!file.exists(cf_nc_fname) || overwrite) { + + tryCatch({ + + # create netCDF file and put arrays + ncout_cf <- ncdf4::nc_create(cf_nc_fname, list(u10.def, v10.def, t2m.def, sp.def, tp.def),force_v4=TRUE) + + # put variables in netCDF file + ncdf4::ncvar_put(ncout_cf,u10.def, cf_u10, verbose = TRUE) + ncdf4::ncvar_put(ncout_cf,v10.def, cf_v10, verbose = TRUE) + ncdf4::ncvar_put(ncout_cf,t2m.def, cf_t2m, verbose = TRUE) + ncdf4::ncvar_put(ncout_cf,sp.def, cf_sp, verbose = TRUE) + ncdf4::ncvar_put(ncout_cf,tp.def, cf_tp, verbose = TRUE) + + # add global attributes + ncdf4::ncatt_put(ncout_cf,0,"GRIB_edition",2) + ncdf4::ncatt_put(ncout_cf,0,"Institution","ecmwf") + ncdf4::ncatt_put(ncout_cf,0,"Source","European Centre for Medium-Range Weather Forecasts") + history <- paste("PEcAn Project", date(Sys.Date()), sep=", ") + ncdf4::ncatt_put(ncout_cf,0,"Conventions","CF-1.7") + + ncdf4::nc_close(ncout_cf) # writes ensemble wise pf files to disk + + }, + + error = function(e) { + PEcAn.logger::logger.severe("Something went wrong during the writing of the nc file.", + conditionMessage(e)) + }) + + } else { + PEcAn.logger::logger.info(paste0( + "The file ", + flname, + " already exists. It was not overwritten." + )) + } + + rows <- 51 + results <- data.frame( + file = character(rows), + host = character(rows), + mimetype = character(rows), + formatname = character(rows), + startdate = character(rows), + enddate = character(rows), + dbfile.name = "ECMWF", + stringsAsFactors = FALSE + ) + + start_date <- latest_filedate + end_date <- strtoi(gsub("-", "", lubridate::ymd(start_date) + days(14))) + + results$file[1] <- + file.path(outfolder, cf_nc_fname) + results$host[1] <- PEcAn.remote::fqdn() + results$startdate[1] <- start_date + results$enddate[1] <- end_date + results$mimetype[1] <- "application/x-netcdf" + results$formatname[1] <- "CF Meteorology" + + ######### Perturbed Forecast ######### + + ##### 3h PF -> 0 - 141h at 3h timestep + pf_3h_out_fname <- paste(latest_filedate, time, stream, "pf_3h", sep = "_") + pf_3h <- combine_files_pf(fnames_3h, lat_in, lon_in) + names(pf_3h) <- c("valid_time", "lon", "lat", "u10", "v10", "t2m", "sp", "tp") + + days_since <- paste("days since", pf_3h$valid_time$time$data, "00:00:00") + pf_3h$valid_time <- pf_3h$valid_time$values + + ##### 6h PF -> 144h - 360h at 6h timestep + pf_6h_out_fname <- paste(latest_filedate, time, stream, "pf_6h", sep = "_") + pf_6h <- combine_files_pf(fnames_6h, lat_in, lon_in) + names(pf_6h) <- c("valid_time", "lon", "lat", "u10", "v10", "t2m", "sp", "tp") + + pf_6h$valid_time <- pf_6h$valid_time$values + + + ##### Creating ensemble wise netCDF files + + pf_u10_3 <- t(pf_3h$u10) + pf_u10_6 <- t(pf_6h$u10) + + pf_v10_3 <- t(pf_3h$v10) + pf_v10_6 <- t(pf_6h$v10) + + pf_t2m_3 <- t(pf_3h$t2m) + pf_t2m_6 <- t(pf_6h$t2m) + + pf_sp_3 <- t(pf_3h$sp) + pf_sp_6 <- t(pf_6h$sp) + + pf_tp_3 <- t(pf_3h$tp) + pf_tp_6 <- t(pf_6h$tp) + + for (ens in 1:50){ + + ncfname_pf <- paste(latest_filedate, time, stream, "pf_ens", ens, sep = "_") + + if (!file.exists(ncfname_pf) || overwrite) { + + tryCatch({ + + ncfname_pf <- paste0(latest_filedate, time, stream, "pf_ens", ens, sep = "_") + + pf_u10 <- c(pf_u10_3[, ens], nanfill(u10_6[, ens])) + pf_v10 <- c(pf_v10_3[, ens], nanfill(v10_6[, ens])) + pf_t2m <- c(pf_t2m_3[, ens], nanfill(pf_t2m_6[, ens])) + pf_sp <- c(pf_sp_3[, ens], nanfill(pf_sp_6[, ens])) + pf_tp <- c(pf_tp_3[, ens], nanfill(pf_tp_6[, ens])) + + ncout_pf <- ncdf4::nc_create(ncfname_pf, list(u10.def, v10.def, t2m.def, sp.def, tp.def),force_v4=TRUE) + + ncdf4::ncvar_put(ncout_pf,u10.def, pf_u10, verbose = TRUE) + ncdf4::ncvar_put(ncout_pf,v10.def, pf_v10, verbose = TRUE) + ncdf4::ncvar_put(ncout_pf,t2m.def, pf_t2m, verbose = TRUE) + ncdf4::ncvar_put(ncout_pf,sp.def, pf_sp, verbose = TRUE) + ncdf4::ncvar_put(ncout_pf,tp.def, pf_tp, verbose = TRUE) + + # add global attributes + ncdf4::ncatt_put(ncout_pf,0,"GRIB_edition",2) + ncdf4::ncatt_put(ncout_pf,0,"Institution","ecmwf") + ncdf4::ncatt_put(ncout_pf,0,"Source","European Centre for Medium-Range Weather Forecasts") + history <- paste("PEcAn Project", date(Sys.Date()), sep=", ") + ncdf4::ncatt_put(ncout_pf,0,"Conventions","CF-1.7") + + ncdf4::nc_close(ncout_pf) # writes ensemble wise pf files to disk + + row <- ens # final data frame has 1+50 rows pertaining to 1cf + 50pf files + results$file[row+1] <- file.path(outfolder, ncout_pf) + results$host[row+1] <- "PEcAn.remote::fqdn()" + results$startdate[row+1] <- start_date + results$enddate[row+1] <- end_date + results$mimetype[row+1] <- "application/x-netcdf" + results$formatname[row+1] <- "CF Meteorology" + + }, + + error = function(e) { + PEcAn.logger::logger.severe("Something went wrong during the writing of the nc file.", + conditionMessage(e)) + }) + + } else { + PEcAn.logger::logger.info(paste0( + "The file ", + flname, + " already exists. It was not overwritten." + )) + } + } + + return(results) +} \ No newline at end of file diff --git a/modules/data.atmosphere/R/download.ECMWF.R b/modules/data.atmosphere/R/download.ECMWF.R new file mode 100644 index 00000000000..f1078eba62a --- /dev/null +++ b/modules/data.atmosphere/R/download.ECMWF.R @@ -0,0 +1,242 @@ +#' Downloads ECMWF Open Data 15 day forecast data products +#' https://confluence.ecmwf.int/display/DAC/ECMWF+open+data%3A+real-time+forecasts +#' +#' https://github.com/ecmwf/ecmwf-opendata +#' Under the hood, this function uses the Python `ecmwf-opendata` module. +#' The module and dependencies can be accessed via the `reticulate` package. +#' +#' `reticulate::conda_install(c("scipy", "xarray", "eccodes", "cfgrib", "ecmwf-opendata", "dask", "netCDF4"), envname = "r-reticulate", pip = TRUE)` +#' `reticulate::use_condaenv("r-reticulate")` +#' +#' +#' @param outfolder location on disk where outputs will be stored +#' +#' @param lat.in,lon.in site coordinates, decimal degrees (numeric) +#' +#' @param fchour reference time of the forecasts. Values are 00, 12 for 15 day forecast. +#' `fchour` is not passed via upstream functions and we are using 00 by default when the function is called from within the pecan workflow. +#' +#' @param type the type of forecast data. Current values: `cf` (controlled forecast), `pf` (perturbed forecast) +#' +#' @param parameters character vector of product types, or `"all"`. +#' +#' @param reticulate_python Path to Python binary for `reticulate` +#' (passed to [reticulate::use_python()]). If `NULL` (default), use +#' the system default. +#' +#' @param overwrite Logical. If `FALSE` (default), skip any files with +#' the same target name (i.e. same variable) that already exist in +#' `outfolder`. If `TRUE`, silently overwrite existing files. +#' +#' Forecast hour is reference time of the forecasts. Values are 00, 12 for 15 day forecast. +# `fchour` is not passed via upstream functions currently as we are using 00 by default (time <- 0) for function calls from within the pecan workflow +#' +#' @return information about the output file +#' +#' @export +#' +#' @examples +#' +#' \dontrun{ +#' files <- download.ECMWF( +#' "ECMWF_output", +#' lat.in = 0, +#' lon.in = 180, +#' fchour = 00, +#' product = "NULL, +#' type = c("cf", "pf"), +#' parameters = "all", +#' ) +#' } +#' @author Swarnalee Mazumder +#' + +download.ECMWF <- function(outfolder, + lat.in, + lon.in, + fchour = 00, + product = NULL, + type = c("cf", "pf"), + parameters = "all", + overwrite = FALSE, + reticulate_python = NULL,...) + +{ + if (!file.exists(outfolder)) { + dir.create(outfolder, showWarnings = FALSE, recursive = TRUE) + } + + ############## Forecast Hour ################ + + # 15 day forecast hours can be 00h or 12h. Here forecast hour is assigned to 00 + # hence the time variable required by the download function becomes 0 + + # all_hours <- c(00, 12) + # if (any(!hour %in% all_hours)) { + # bad_hours <- setdiff(hour, all_hours) + # PEcAn.logger::logger.severe(sprintf( + # "Invalid forecast hours %s. 15 day forecast hours must be one of the following: %s", + # paste0("`", bad_hours, "`", collapse = ", "), + # paste0("`", all_hours, "`", collapse = ", ") + # )) + # } + + # for forecast hour = 00 + if (fchour == 00){ + time <- 0 + } else { + PEcAn.logger::logger.severe(sprintf( + "Invalid forecast hour. Currently, forecast hour 00 is the default." + )) + } + + ############## Product (Forecasting system) ################ + + # For product = NULL,the stream is set to "enfo" + # stream is forecasting system that produces the data. Current value is `"enfo"` + + if(is.null(product)){ + stream <- "enfo" + } else{ + PEcAn.logger::logger.severe(sprintf( + "Invalid data stream. Currently, data stream `enfo` is supported." + )) + } + + if (stream == "enfo"){ + time <- time + step_15day <- 360 + type <- type # type of forecast here c("cf", "pf") + } else{ + PEcAn.logger::logger.severe(sprintf( + "Invalid data stream. Currently, data stream `enfo` is supported." + )) + } + + + ############## Product Types ################ + + # Currently supported parameters are + # "2t" - 2 metre temperature + # "tp" - total precipitation + # "10u" - 10 metre U wind component + # "10v" - 10 metre V wind component + # "sp" - Surface Pressure + + + all_parameters <- c("2t", "tp", "10u", "10v", "sp") + + if (tolower(parameters) == "all") { + parameter_types <- all_parameters + } + + if (any(!parameter_types %in% all_parameters)) { + bad_parameters <- setdiff(parameter_types, all_parameters) + PEcAn.logger::logger.severe(sprintf( + "Invalid parameter types %s. Products must be one of the following: %s", + paste0("`", bad_parameters, "`", collapse = ", "), + paste0("`", all_parameters, "`", collapse = ", ") + )) + } + + + # Python script "download_ecmwf.py" to get latest forecast date and download forecast datasets + script.path = file.path(system.file("ECMWF/download_ecmwf.py", package = "PEcAn.data.atmosphere")) + reticulate::source_python(script.path) + + date_latestdata <- ecmwflatest(time, step, stream, type, all_parameters) + latest_filedate <- strtoi(gsub("-", "", substr(as.character.Date(date_latestdata), 1, 10))) + + current_date <- strtoi(gsub("-", "", Sys.Date())) + + date <- if (latest_filedate == current_date){ + # today + 0 + } else if ((latest_filedate - current_date) == -1){ + # yesterday + -1 + } else if ((latest_filedate - current_date) == -2){ + # the day before yesterday + -2 + } else { PEcAn.logger::logger.severe(sprintf( + "Invalid file download date." + )) + } + + # Skip download if there is any already existing file with base name same as `latest_filedate`` + if ((length(list.files(path = outfolder, pattern = as.character(latest_filedate)) ) > 0) == TRUE){ + PEcAn.logger::logger.severe(sprintf( + "Files already exist for %s", + latest_filedate + )) + } else { + + in_fname <- paste(latest_filedate, time, stream, sep = "_") + + data_download <- ecmwfdownload(date, time, stream, type, all_parameters, in_fname) + + rows <- 85 # 48 files at 3h timestep and 36 files at 6h timestep + results <- data.frame( + file = character(rows), + host = character(rows), + mimetype = character(rows), + formatname = character(rows), + startdate = character(rows), + enddate = character(rows), + dbfile.name = "ECMWF", + stringsAsFactors = FALSE + ) + + out_3h <- list() + for (h3 in seq(0, 141, 3)){ + out_3h <- append(out_3h, paste0(in_fname,"_",h3,".grib2")) + } + + out_6h <- list() + for (h6 in seq(144, 360, 6)){ + out_6h <- append(out_6h, paste0(in_fname,"_",h6,".grib2")) + } + + for (row in 1:48){ + results$file[row] <- + file.path(outfolder, out_3h[[row]]) + results$host[row] <- PEcAn.remote::fqdn() + results$startdate[row] <- "start_date" + results$enddate[row] <- "end_date" + results$mimetype[row] <- "application/x-netcd" + results$formatname[row] <- "CF Meteorology" + } + + for (row in 1:37){ + results$file[row+48] <- + file.path(outfolder, out_6h[[row]]) + results$host[row+48] <- PEcAn.remote::fqdn() + results$startdate[row+48] <- "start_date" + results$enddate[row+48] <- "end_date" + results$mimetype[row+48] <- "application/x-netcd" + results$formatname[row+48] <- "CF Meteorology" + } + return(results) + } + +} + +######### code to reproduce +# library(reticulate) +# reticulate::conda_install(c("scipy", "xarray", "eccodes", "cfgrib", "ecmwf-opendata", "dask", "netCDF4"), envname = "r-reticulate", pip = TRUE) +# reticulate::use_condaenv("r-reticulate") +# +# date <- -1 +# time <- 0 +# stream <- "enfo" +# type <- c("cf", "pf") +# step_15day <- 360 +# +# all_parameters <- c("10u", "10v", "2t", "sp", "tp") +# fname <- "15day_ecmwfxxx" +# +# reticulate::source_python("downloadecmwf.py") +# +# latest <- ecmwflatest(time, step_15day, stream, type, all_parameters) +# +# data_download <- ecmwfdownload(date= date, time= time, stream= stream, type= type, params= all_parameters, filename= fname) \ No newline at end of file diff --git a/modules/data.atmosphere/R/met2CF.ECMWF.R b/modules/data.atmosphere/R/met2CF.ECMWF.R new file mode 100644 index 00000000000..ecaa7932aa5 --- /dev/null +++ b/modules/data.atmosphere/R/met2CF.ECMWF.R @@ -0,0 +1,36 @@ +#' met2CF.ECMWF +#' +#' Converts ECMWF Open Data variables to CF format. +#' +#' Variables present in the output netCDF files: +#' air_temperature, air_temperature, air_pressure, +#' precipitation_flux, eastward_wind, northward_wind +#' +#' @param lat.in latitude +#' @param lon.in longitude +#' @param outfolder Path to directory where nc files need to be saved. +#' @param overwrite Logical if files needs to be overwritten. +#' @param verbose Logical flag defining if ouput of function be extra verbose. +#' +#' @return list of files in a dataframe +#' @export +#' +#' @author Swarnalee Mazumder +#' +met2CF.ECMWF <- function(lat.in, + lon.in, + outfolder, + overwrite = FALSE, + verbose = TRUE) { + + if (!file.exists(outfolder)) { + dir.create(outfolder, showWarnings = FALSE, recursive = TRUE) + } + + results <- + PEcAn.data.atmosphere::combine2nc_ECMWF(lat.in, + lon.in, + outfolder, + overwrite = overwrite) + return(results) +} \ No newline at end of file diff --git a/modules/data.atmosphere/inst/ECMWF/download_ecmwf.py b/modules/data.atmosphere/inst/ECMWF/download_ecmwf.py new file mode 100644 index 00000000000..18f6b4c6c01 --- /dev/null +++ b/modules/data.atmosphere/inst/ECMWF/download_ecmwf.py @@ -0,0 +1,104 @@ +#This script helps in viewing the latest date for which 15 day forecast data is available on https://data.ecmwf.int/ and helps in downloading the same +#Author: Swarnalee Mazumder + +def ecmwflatest(time, step15, stream, type, params): + + from ecmwf.opendata import Client + + """ + Fetches latest ECMWF Open Data forecast based on specified forecast time, step, stream, parameters + + Parameters + ---------- + time (numeric) -- reference time of the forecasts. Values are 0 (00), 12. + + step (numeric) -- forecast time step. Current value: 360 + + stream (str) -- forecasting system that produces the data. Current value: `enfo` + + lat_in (numeric) -- Site coordinates, decimal degrees + + lon_in (numeric) -- Site coordinates, decimal degrees + + Returns + ------- + Latest date for which data is available. + """ + from ecmwf.opendata import Client + + client = Client("ecmwf", beta=True) + + if stream == "mmsf": + latest= client.latest( + time = int(time), + step= int(step15), + stream= stream, + type= type, + param= params, + ) + + elif stream == "enfo": + latest= client.latest( + time = int(time), + step= int(step15), + stream= stream, + type= type, + param= params, + ) + + return latest + + +def ecmwfdownload(date, time, stream, type, params, filename): + + """ + Downloads latest ECMWF Open Data forecast based on specified forecast time, step, stream, parameters + + Parameters + ---------- + date (numeric) -- reference date of the forecasts. Values are 0 (today), -1 (yesterday), -2 (day before yesterday) + + time (numeric) -- reference time of the forecasts. Values are 0 (00), 12. + + stream (str) -- forecasting system that produces the data. Current value: `enfo` + + type (str) -- the type of data. Current values: `cf`, `pf` + + params (str) -- the meteorological parameters, such as wind, pressure or humidity. Current parameters are "2t", "tp", "10u", "10v", "q", "r", "sp". + + filename (str) -- output filename for the Grib2 file + + Returns + ------- + Downloads ECMWF Open Data forecast for the latest date for which forecast data is available. + """ + + from ecmwf.opendata import Client + + client = Client("ecmwf", beta=True) + + for i in range(0, 142, 3): + file_3h = client.retrieve( + date= int(date), + time= int(time), + step= int(i), + stream= stream, + type= type, + param= params, + target= filename+str(i)+".grib2" + ) + + + for j in range(144, 361, 6): + file_6h = client.retrieve( + date= int(date), + time= int(time), + step= int(j), + stream= stream, + type= type, + param= params, + target= filename+str(j)+".grib2" + ) + + return file_3h, file_6h + diff --git a/modules/data.atmosphere/inst/ECMWF/met2CFutils.py b/modules/data.atmosphere/inst/ECMWF/met2CFutils.py new file mode 100644 index 00000000000..8957559dc52 --- /dev/null +++ b/modules/data.atmosphere/inst/ECMWF/met2CFutils.py @@ -0,0 +1,140 @@ +#The functions contained in this script help in extraction of CF and PF data from GRIB2 files and conversion to CF standard +#Author: Swarnalee Mazumder + +def all_filenames_3h(filename): + files_3 = [] + + for i in range(0, 142,3): + files_3.append(filename+"_"+str(i)+".grib2") + + return files_3 + +def all_filenames_6h(filename): + files_6 = [] + for j in range(144, 366, 6): + files_6.append(filename+"_"+str(j)+".grib2") + + return files_6 + +def combine_files_cf(files, lat_in, lon_in): + + """ + Converts ECMWF Open Data 15 day controlled forecast data from Grib2 into PEcAn standard NetCDF + + Parameters + ---------- + in_filename (str) -- Grib2 file consisting of forecast data + + lat_in (numeric) -- Site coordinates, decimal degrees + + lon_in (numeric) -- Site coordinates, decimal degrees + + Returns + ------- + Arrays of variables + """ + + import os + import xarray + import cfgrib + + lat_in = float(lat_in) + lon_in = float(lon_in) + + + cf_u10 = xarray.open_mfdataset(files, engine = 'cfgrib', concat_dim="valid_time",combine="nested", + backend_kwargs={'filter_by_keys': {'dataType': 'cf', 'shortName':'10u'}}) + + cf_v10 = xarray.open_mfdataset(files, engine = 'cfgrib', concat_dim="valid_time",combine="nested", + backend_kwargs={'filter_by_keys': {'dataType': 'cf', 'shortName':'10v'}}) + + cf_t2m = xarray.open_mfdataset(files, engine = 'cfgrib', concat_dim="valid_time",combine="nested", + backend_kwargs={'filter_by_keys': {'dataType': 'cf', 'shortName':'2t'}}) + + cf_sp = xarray.open_mfdataset(files, engine = 'cfgrib', concat_dim="valid_time",combine="nested", + backend_kwargs={'filter_by_keys': {'dataType': 'cf', 'shortName':'sp'}}) + + cf_tp = xarray.open_mfdataset(files, engine = 'cfgrib', concat_dim="valid_time",combine="nested", + backend_kwargs={'filter_by_keys': {'dataType': 'cf', 'shortName':'tp'}}) + + cf = xarray.merge([cf_u10, cf_v10, cf_t2m, cf_sp, cf_tp], + compat = 'override') + + # total precipitation unit as per CF standard + # all other variables already in CF standard units + cf.tp.attrs.update({'GRIB_units': 'kg m**-2 s**-1', 'units': 'kg m**-2 s**-1'}) + attrs_tp_cf = cf.tp.attrs + + + # 1 m = 1/86.4 kg/m2/s + cf['tp'] = cf.tp / 86.4 + cf.tp.attrs = attrs_tp_cf + + + cf_lat_lon = cf.sel(latitude = lat_in, longitude = lon_in, method = "nearest") + + return cf_lat_lon["valid_time"], cf_lat_lon["longitude"].values, cf_lat_lon["latitude"].values, cf_lat_lon['u10'].values, cf_lat_lon["v10"].values, cf_lat_lon["t2m"].values, cf_lat_lon["sp"].values, cf_lat_lon["tp"].values + + +def combine_files_pf(files, lat_in, lon_in): + + """ + Converts ECMWF Open Data 15 day perturbed forecast from Grib2 into PEcAn standard NetCDF + + Parameters + ---------- + in_filename (str) -- Grib2 file consisting of forecast data + + lat_in (numeric) -- Site coordinates, decimal degrees + + lon_in (numeric) -- Site coordinates, decimal degrees + + Returns + ------- + Arrays of variables + """ + import os + import xarray + import cfgrib + + lat_in = float(lat_in) + lon_in = float(lon_in) + + ######### Perturbed Forecast + # 10 metre U wind component + pf_u10 = xarray.open_mfdataset(files, engine='cfgrib', concat_dim="valid_time",combine="nested", + backend_kwargs={'filter_by_keys': {'dataType': 'pf', 'shortName':'10u'}}) + + # 10 metre V wind component + pf_v10 = xarray.open_mfdataset(files, engine='cfgrib', concat_dim="valid_time",combine="nested", + backend_kwargs={'filter_by_keys': {'dataType': 'pf', 'shortName':'10v'}}) + + # 2 metre temperature + pf_t2m = xarray.open_mfdataset(files, engine='cfgrib', concat_dim="valid_time",combine="nested", + backend_kwargs={'filter_by_keys': {'dataType': 'pf', 'shortName':'2t'}}) + + # Surface pressure + pf_sp = xarray.open_mfdataset(files, engine='cfgrib', concat_dim="valid_time",combine="nested", + backend_kwargs={'filter_by_keys': {'dataType': 'pf', 'shortName':'sp'}}) + + # total precipitation + pf_tp = xarray.open_mfdataset(files, engine='cfgrib', concat_dim="valid_time",combine="nested", + backend_kwargs={'filter_by_keys': {'dataType': 'pf', 'shortName':'tp'}}) + + + pf = xarray.merge([pf_u10, pf_v10, pf_t2m, pf_sp, pf_tp], + compat = 'override') + + + pf.tp.attrs.update({'GRIB_units': 'kg m**-2 s**-1', 'units': 'kg m**-2 s**-1'}) + attrs_tp_pf = pf.tp.attrs + + # 1 m = 1/86.4 kg/m2/s + pf['tp'] = pf.tp / 86.4 + pf.tp.attrs = attrs_tp_pf + + + pf_lat_lon = pf.sel(latitude = lat_in, longitude = lon_in, method = "nearest") + + return pf_lat_lon["valid_time"], pf_lat_lon["longitude"].values, pf_lat_lon["latitude"].values, pf_lat_lon['u10'].values, pf_lat_lon["v10"].values, pf_lat_lon["t2m"].values, pf_lat_lon["sp"].values, pf_lat_lon["tp"].values + diff --git a/modules/data.atmosphere/inst/registration/register.ECMWF.xml b/modules/data.atmosphere/inst/registration/register.ECMWF.xml new file mode 100644 index 00000000000..726f4dcdb9e --- /dev/null +++ b/modules/data.atmosphere/inst/registration/register.ECMWF.xml @@ -0,0 +1,12 @@ + + + regional + 51 + TRUE + + 33 + CF Meteorology + application/x-netcdf + nc + + \ No newline at end of file diff --git a/modules/data.atmosphere/man/download.ECMWF.Rd b/modules/data.atmosphere/man/download.ECMWF.Rd new file mode 100644 index 00000000000..fa999e4b8c2 --- /dev/null +++ b/modules/data.atmosphere/man/download.ECMWF.Rd @@ -0,0 +1,57 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/download.ECMWF.R +\name{download.ECMWF} +\alias{download.ECMWF} +\title{Downloads ECMWF Open Data 15 day forecast data products +https://confluence.ecmwf.int/display/DAC/ECMWF+open+data%3A+real-time+forecasts} +\usage{ +download.ECMWF( + outfolder, + lat.in, + lon.in, + fchour = 0, + product = NULL, + type = c("cf", "pf"), + parameters = "all", + overwrite = FALSE, + reticulate_python = NULL, + ... +) +} +\arguments{ +\item{outfolder}{location on disk where outputs will be stored} + +\item{lat.in, lon.in}{site coordinates, decimal degrees (numeric)} + +\item{fchour}{reference time of the forecasts. Values are 00, 12 for 15 day forecast. +`fchour` is not passed via upstream functions and we are using 00 by default when the function is called from within the pecan workflow.} + +\item{type}{the type of forecast data. Current values: `cf` (controlled forecast), `pf` (perturbed forecast)} + +\item{parameters}{character vector of product types, or `"all"`.} + +\item{overwrite}{Logical. If `FALSE` (default), skip any files with + the same target name (i.e. same variable) that already exist in + `outfolder`. If `TRUE`, silently overwrite existing files. + +Forecast hour is reference time of the forecasts. Values are 00, 12 for 15 day forecast.} + +\item{reticulate_python}{Path to Python binary for `reticulate` +(passed to [reticulate::use_python()]). If `NULL` (default), use +the system default.} +} +\value{ +information about the output file +} +\description{ +https://github.com/ecmwf/ecmwf-opendata +Under the hood, this function uses the Python `ecmwf-opendata` module. +The module and dependencies can be accessed via the `reticulate` package. +} +\details{ +`reticulate::conda_install(c("scipy", "xarray", "eccodes", "cfgrib", "ecmwf-opendata"), envname = "r-reticulate", pip = TRUE)` +`reticulate::use_condaenv("r-reticulate")` +} +\author{ +Swarnalee Mazumder +}