refactor ecmwf transform function to convert units, calculate relativ…

…e humidity, and remove unused data

R/calculate_forecasts_anomalies.R

@@ -0,0 +1,25 @@
+#' .. content for \description{} (no empty lines) ..
+#'
+#' .. content for \details{} ..
+#'
+#' @title
+#' @param ecmwf_forecasts_transformed
+#' @param ecmwf_forecasts_transformed_directory
+#' @param weather_historical_means
+#' @param forecast_anomalies_directory
+#' @param model_dates
+#' @param model_dates_selected
+#' @param overwrite
+#' @return
+#' @author Emma Mendelsohn
+#' @export
+calculate_forecasts_anomalies <- function(ecmwf_forecasts_transformed,
+                                          ecmwf_forecasts_transformed_directory,
+                                          weather_historical_means,
+                                          forecast_anomalies_directory,
+                                          model_dates, model_dates_selected,
+                                          overwrite = FALSE) {
+
+  NULL
+
+}

R/download_ecmwf_forecasts.R

@@ -7,15 +7,15 @@
 #' @return
 #' @author Emma Mendelsohn
 #' @export
-download_ecmwf_forecasts <- function(ecmwf_api_parameters,
+download_ecmwf_forecasts <- function(ecmwf_forecasts_api_parameters,
                                      download_directory,
                                      overwrite = FALSE){
 
   existing_files <- list.files(download_directory)
 
-  system <- ecmwf_api_parameters$system
-  year <- ecmwf_api_parameters$year
-  month <- unlist(ecmwf_api_parameters$month)
+  system <- ecmwf_forecasts_api_parameters$system
+  year <- ecmwf_forecasts_api_parameters$year
+  month <- unlist(ecmwf_forecasts_api_parameters$month)
 
   filename <- paste("ecmwf", "seasonal_forecast", paste0("sys", system), year, sep = "_")
   filename <- paste0(filename, ".grib")
@@ -29,12 +29,12 @@ download_ecmwf_forecasts <- function(ecmwf_api_parameters,
   request <- list(
     originating_centre = "ecmwf",
     system = system,
-    variable = unlist(ecmwf_api_parameters$variables),
-    product_type = unlist(ecmwf_api_parameters$product_types),
+    variable = unlist(ecmwf_forecasts_api_parameters$variables),
+    product_type = unlist(ecmwf_forecasts_api_parameters$product_types),
     year = year,
     month = month,
-    leadtime_month = unlist(ecmwf_api_parameters$leadtime_months),
-    area = round(unlist(ecmwf_api_parameters$spatial_bounds), 1),
+    leadtime_month = unlist(ecmwf_forecasts_api_parameters$leadtime_months),
+    area = round(unlist(ecmwf_forecasts_api_parameters$spatial_bounds), 1),
     format = "grib",
     dataset_short_name = "seasonal-monthly-single-levels",
     target = filename

R/transform_ecmwf_forecasts.R

@@ -47,30 +47,92 @@ transform_ecmwf_forecasts <- function(ecmwf_forecasts_downloaded,
     mutate(data_date = ymd(data_date))  |> 
     select(-grid_type, -packing_type, -level, -type_of_level, -centre, -edition)
 
-  # Calculate per-pixel mean and std for each unique combination, across all models
+  # Rename layers with the metadata names
+  names(grib) <- meta$variable_id
+
+  # Select only the means columns
+  grib_subset <- subset(grib, which(str_detect(names(grib), "fcmean"))) 
+
+  # Units conversions
+  # 2d = "2m_dewpoint_temperature" = 	2 metre dewpoint temperature K
+  # 2t = "2m_temperature" = 	2 metre temperature K
+  # tprate = "total_precipitation" =  Total precipitation m/s
+  # https://cds.climate.copernicus.eu/cdsapp#!/dataset/seasonal-monthly-single-levels?tab=overview  
+
+  # Convert kelvin to celsius (note these columns are mislabeled as C)
+  temp_cols <- which(str_ends(names(grib_subset), "2d|2t"))
+  grib_temp <- subset(grib_subset, temp_cols)
+  grib_temp <- app(grib_temp, function(i) i - 273.15)
+
+  # Convert precipitation from m/second to mm/day
+  precip_cols <- which(str_ends(names(grib_subset), "tprate"))
+  grib_precip <- subset(grib_subset, precip_cols)
+  grib_precip <- ifel(grib_precip < 0, 0, grib_precip) # some very small negative values
+  grib_precip <- app(grib_precip, function(i) i * 8.64e+7)
+
+  # Calculate per-pixel mean for each unique combination, across all models
   # transform to template
-  grib_index <- as.integer(factor(meta$variable_id, levels = unique(meta$variable_id)))
-  grib_means <- terra::tapp(grib, grib_index, "mean", cores = n_workers) |> 
-    setNames(unique(meta$variable_id)) |> 
+  grib_converted <- c(grib_precip, grib_temp)
+  grib_index <- as.integer(factor(names(grib_converted), levels = unique(names(grib_converted))))
+  grib_means <- terra::tapp(grib_converted, grib_index, "mean", cores = n_workers) |> 
+    setNames(unique(names(grib_converted))) |> 
     transform_raster(template = continent_raster_template) |> 
-    as.data.frame(xy = TRUE) |> 
+    as.data.frame(xy = TRUE) 
+
+  # Using the means, calculate relative humidity from dew point and temp
+  dp_cols <- names(grib_means)[str_ends(names(grib_means), "2d")] |> sort() # dew point
+  t_cols <- names(grib_means)[str_ends(names(grib_means), "2t")] |> sort() # temp
+  assertthat::assert_that(identical(str_remove(dp_cols, "_2d"), str_remove(t_cols, "_2t")))
+
+  rel_humidity <- map2_dfc(dp_cols, t_cols, function(dp, t){
+    dp <- grib_means[,dp]
+    t <- grib_means[, t]
+    rh <- 100 * exp((17.625 * dp)/(243.04+dp))/exp((17.625 * t)/(243.04+t))
+    assertthat::assert_that(all(rh <= 100 & rh >=0))
+    return(rh)
+  }) |> set_names(str_replace(dp_cols, "_2d", "_rh"))
+
+  # Remove dewpoint temperature and add relative humidity
+  grib_means <- grib_means |> 
+    select(-dp_cols) |> 
+    bind_cols(rel_humidity)|> 
     reshape2::melt(id.vars = c("x", "y"), variable.name = "variable_id", value.name = "mean") 
-  grib_sds <- terra::tapp(grib, grib_index, "sd", cores = n_workers) |> 
-    setNames(unique(meta$variable_id)) |> 
-    transform_raster(template = continent_raster_template) |> 
-    as.data.frame(xy = TRUE) |> 
-    reshape2::melt(id.vars = c("x", "y"), variable.name = "variable_id", value.name = "std") |> 
-    pull(std)
-
-  # Link with metadata for export
-  dat_out <- grib_means |> 
-    mutate(std = grib_sds) |> 
-    left_join(distinct(meta), by = "variable_id") |> 
+
+  # Not using Sds at the moment
+  # grib_sds <- terra::tapp(grib_converted, grib_index, "sd", cores = n_workers) |> 
+  #   setNames(unique(names(grib_converted))) |> 
+  #   transform_raster(template = continent_raster_template) |> 
+  #   as.data.frame(xy = TRUE) |> 
+  #   reshape2::melt(id.vars = c("x", "y"), variable.name = "variable_id", value.name = "std") |> 
+  #   pull(std)
+
+  # Add relative humidity to metafile 
+  meta_out <- meta |> 
+    distinct() |> 
+    mutate(short_name = if_else(short_name == "2d", "rh", short_name)) |> 
+    mutate(variable_id = str_replace(variable_id, "_2d", "_rh"))
+
+  # Link means with metadata for export
+  grib_means <- grib_means |> 
+    left_join(meta_out, by = "variable_id") |> 
     arrange(variable_id, x, y) |> 
-    select(x,y, mean, std, data_date, step_range, data_type, short_name)
+    select(x, y, mean, data_date, step_range, short_name)
+
+  # Standardize leads
+  n_leads <- grib_means |> 
+    group_by(data_date) |> 
+    summarize(n = n_distinct(step_range)) |> 
+    distinct(n)
+  assertthat::are_equal(6, n_leads$n)
 
+  grib_means <- grib_means |> 
+    group_by(data_date) |> 
+    mutate(lead_month = factor(rank(step_range), labels = 1:6)) |> 
+    select(-step_range) |> 
+    ungroup()
+
   # Save as parquet 
-  write_parquet(dat_out, here::here(ecmwf_forecasts_transformed_directory, save_filename), compression = "gzip", compression_level = 5)
+  write_parquet(grib_means, here::here(ecmwf_forecasts_transformed_directory, save_filename), compression = "gzip", compression_level = 5)
 
   return(file.path(ecmwf_forecasts_transformed_directory, save_filename))
 

_targets.R

@@ -298,7 +298,7 @@ data_targets <- tar_plan(
   tar_target(model_dates_selected, model_dates |> filter(select_date) |> pull(date)),
 
 
-  # weather anomalies --------------------------------------------------
+  # recorded weather anomalies --------------------------------------------------
   tar_target(weather_historical_means_directory, 
              create_data_directory(directory_path = "data/weather_historical_means")),
 
@@ -335,6 +335,22 @@ data_targets <- tar_plan(
              pattern = weather_anomalies,
              cue = tar_cue("never")), # only run this if you need to upload new data  
 
+
+  # forecast weather anomalies ----------------------------------------------------------------------
+  tar_target(forecasts_anomalies_directory, 
+             create_data_directory(directory_path = "data/forecast_anomalies")),
+
+  tar_target(forecasts_anomalies, calculate_forecasts_anomalies(ecmwf_forecasts_transformed,
+                                                                ecmwf_forecasts_transformed_directory,
+                                                                weather_historical_means,
+                                                                forecast_anomalies_directory,
+                                                                model_dates,
+                                                                model_dates_selected,
+                                                                overwrite = FALSE),
+             pattern = model_dates_selected,
+             format = "file", 
+             repository = "local"),    
+
   # ndvi anomalies --------------------------------------------------
   tar_target(ndvi_date_lookup, 
              create_ndvi_date_lookup(sentinel_ndvi_transformed,