fixing elevation and adding DEM func

DESCRIPTION

@@ -22,6 +22,7 @@ Imports:
   cowplot,
   countrycode,
   dplyr,
+  elevatr,
   ggplot2,
   ggraph,
   ggrepel,
@@ -36,6 +37,7 @@ Imports:
   magrittr,
   mobility,
   propvacc,
+  raster,
   reshape2,
   rnaturalearth,
   sf,

R/download_country_DEM.R

@@ -0,0 +1,78 @@
+#' Download 1km DEM Raster for Multiple Countries and Save to Files
+#'
+#' This function downloads 1km resolution DEM (Digital Elevation Model) raster data for a given vector of ISO3 country codes, and writes each raster to a GeoTIFF file.
+#'
+#' @param PATHS A list containing paths where DEM data will be saved. Typically generated by the `get_paths()` function and should include:
+#' \itemize{
+#'   \item \strong{DATA_DEM}: Path to the directory where DEM rasters will be saved.
+#' }
+#' #' @param iso_codes A character vector of ISO3 country codes for which DEM data should be downloaded.
+#' @param zoom_level An integer representing the zoom level for the DEM data. Zoom levels 6-8 are recommended for 1km resolution (default = 6).
+#'
+#' @return The function does not return a value. It downloads the DEM data for each country and saves the results as GeoTIFF files in the specified directory.
+#'
+#' @importFrom elevatr get_elev_raster
+#' @importFrom sf st_bbox
+#' @importFrom rnaturalearth ne_countries
+#' @importFrom raster writeRaster
+#' @importFrom glue glue
+#' @examples
+#' \dontrun{
+#' # Define the ISO3 country codes
+#' iso_codes <- c("ZAF", "KEN", "NGA")
+#'
+#' # Get paths for data storage
+#' PATHS <- get_paths()
+#'
+#' # Download DEM rasters for the countries
+#' download_country_DEM(iso_codes, PATHS)
+#' }
+#' @export
+
+download_country_DEM <- function(PATHS, iso_codes, zoom_level = 6) {
+
+     requireNamespace('sf')
+     requireNamespace('raster')
+
+     # Ensure the output directory exists
+     if (!dir.exists(PATHS$DATA_DEM)) dir.create(PATHS$DATA_DEM, recursive = TRUE)
+
+
+     # Loop through each ISO3 country code
+     for (i in 1:length(iso_codes)) {
+
+          message(glue::glue("Downloading DEM for {iso_codes[i]}..."))
+
+          # Construct the path to the shapefile for each country
+          shapefile_path <- file.path(PATHS$DATA_SHAPEFILES, paste0(iso_codes[i], "_ADM0.shp"))
+
+          # Check if the shapefile exists
+          if (!file.exists(shapefile_path)) {
+               message(paste("Shapefile for", iso_codes[i], "not found. Skipping."))
+               next
+          }
+
+          # Load the country shapefile
+          country_shapefile <- sf::st_read(shapefile_path, quiet = TRUE)
+
+
+          if (is.null(country_shapefile) || nrow(country_shapefile) == 0) {
+               message(glue::glue("Skipping {iso_codes[i]}: shapefile not found."))
+               next
+          }
+
+          # Download the DEM data for the country's bounding box
+          dem_raster <- elevatr::get_elev_raster(locations = country_shapefile, z = zoom_level, clip = "bbox")
+
+          # Create a filename for the DEM file
+          dem_filename <- file.path(PATHS$DATA_DEM, glue::glue("{iso_codes[i]}_1km_DEM.tif"))
+
+          # Save the DEM raster as a GeoTIFF file
+          raster::writeRaster(dem_raster, filename = dem_filename, format = "GTiff", overwrite = TRUE)
+
+          message(glue::glue("DEM for {iso_codes[i]} saved to {dem_filename}."))
+
+     }
+
+     message(glue::glue("DEM data saved for all countries in {PATHS$DATA_DEM}."))
+}

R/est_seasonal_dynamics.R

@@ -0,0 +1,97 @@
+#' Get Seasonal Dynamics Data for Cholera and Precipitation
+#'
+#' This function retrieves historical precipitation data, processes cholera case data, fits seasonal dynamics models using Fourier series,
+#' and saves the results to a CSV file.
+#'
+#' @param PATHS A list containing paths where raw and processed data are stored. Typically generated by the `get_paths()` function.
+#'
+#' @return The function saves the parameter estimates, fitted values, and processed data as CSV files.
+#'
+#' @importFrom dplyr select mutate group_by summarize filter bind_rows
+#' @importFrom tidyr spread
+#' @importFrom lubridate week year
+#' @importFrom sf st_read st_coordinates
+#' @importFrom utils write.csv read.csv
+#' @importFrom glue glue
+#' @importFrom ggplot2 ggplot geom_sf labs theme_minimal
+#' @importFrom minpack.lm nlsLM
+#' @export
+
+est_seasonal_dynamics <- function(PATHS) {
+
+     iso_codes <- MOSAIC::iso_codes_mosaic
+
+     # Load cholera cases data
+     cholera_data <- utils::read.csv(file.path(PATHS$DATA_WHO_WEEKLY, "cholera_country_weekly_processed.csv"), stringsAsFactors = FALSE)
+
+     iso_codes_with_data <- sort(unique(cholera_data$iso_code))
+     iso_codes_no_data <- setdiff(iso_codes, unique(cholera_data$iso_code))
+
+     all_precip_data <- list()
+     all_fitted_values <- list()
+     all_param_values <- list()
+
+     for (country_iso_code in iso_codes) {
+
+          message(glue::glue("Processing country: {country_iso_code}"))
+
+
+          # instead load climate data
+
+          country_shp <- sf::st_read(file.path(PATHS$DATA_SHAPEFILES, paste0(country_iso_code, "_ADM0.shp")), quiet = TRUE)
+          grid_points <- MOSAIC::generate_country_grid_n(country_shp, n_points = 10)
+          coords <- sf::st_coordinates(grid_points)
+          coords <- as.data.frame(coords)
+
+          precipitation_data_list <- list()
+
+          for (i in seq_len(nrow(coords))) {
+               lat <- coords$Y[i]
+               lon <- coords$X[i]
+
+               precip_data <- get_historical_precip(lat, lon, as.Date("2014-09-01"), as.Date("2024-09-01"), api_key = "your_api_key_here")
+               precip_data$year <- lubridate::year(precip_data$date)
+               precip_data$week <- lubridate::week(precip_data$date)
+
+               weekly_precip <- precip_data %>%
+                    dplyr::group_by(year, week) %>%
+                    dplyr::summarize(weekly_precipitation_sum = sum(daily_precipitation_sum, na.rm = TRUE))
+
+               weekly_precip$iso_code <- country_iso_code
+               precipitation_data_list[[i]] <- weekly_precip
+          }
+
+          precip_data <- do.call(rbind, precipitation_data_list)
+          precip_data <- merge(precip_data, cholera_data, by = c("week", "iso_code"), all.x = TRUE)
+
+          # Scale precipitation and cholera data
+          precip_data <- precip_data %>%
+               dplyr::mutate(precip_scaled = scale(weekly_precipitation_sum),
+                             cases_scaled = scale(cases))
+
+          # Fit Fourier series model for precipitation
+          fit_precip <- minpack.lm::nlsLM(
+               precip_scaled ~ generalized_fourier(week, beta0, a1, b1, a2, b2, p = 52),
+               data = precip_data,
+               start = list(beta0 = 0, a1 = 1, b1 = 1, a2 = 0.5, b2 = 0.5)
+          )
+
+          param_precip <- coef(summary(fit_precip))
+          param_values_precip <- data.frame(
+               country_iso_code = country_iso_code,
+               response = "precipitation",
+               parameter = rownames(param_precip),
+               mean = param_precip[,"Estimate"],
+               se = param_precip[,"Std. Error"],
+               ci_lo = param_precip[,"Estimate"] - param_precip[,"Std. Error"] * 1.96,
+               ci_hi = param_precip[,"Estimate"] + param_precip[,"Std. Error"] * 1.96
+          )
+
+          all_param_values[[country_iso_code]] <- param_values_precip
+          all_precip_data[[country_iso_code]] <- precip_data
+     }
+
+     combined_param_values <- do.call(rbind, all_param_values)
+     utils::write.csv(combined_param_values, file = file.path(PATHS$MODEL_INPUT, "combined_param_values.csv"), row.names = FALSE)
+     message("Parameter values saved.")
+}

R/get_elevation.R

@@ -1,110 +1,110 @@
-#' Get Median Elevation Data for Multiple Countries Using Shapefiles
+#' Get Mean or Median Elevation from Downloaded DEM Data Using Country Boundaries
 #'
-#' This function retrieves elevation data for a specified number of points within the boundary of each country's shapefile from the `iso_codes_mosaic` list using the Open-Meteo API. It returns the median elevation for each country.
+#' This function calculates the mean or median elevation for a given set of ISO3 country codes by loading the downloaded DEM files and extracting raster values within the country's boundary.
 #'
-#' @param PATHS A list containing the paths to where the elevation data should be saved. Typically generated by the `get_paths()` function and should include:
+#' @param iso_codes A character vector of ISO3 country codes for which the elevation data will be calculated.
+#' @param PATHS A list containing paths where DEM data is stored. Typically generated by the `get_paths()` function and should include:
 #' \itemize{
-#'   \item \strong{DATA_SHAPEFILES}: Path to the directory containing country shapefiles.
-#'   \item \strong{DATA_ELEVATION}: Path to the directory where the elevation data will be saved.
+#'   \item \strong{DATA_DEM}: Path to the directory where DEM rasters are saved (GeoTIFF format).
+#'   \item \strong{DATA_SHAPEFILES}: Path to the directory where shapefiles are saved.
 #' }
-#' @param n_points An integer specifying the number of points to generate within each country's boundary.
-#' @param api_key A character string representing the API key for the Open-Meteo API.
+#' @param method A character string specifying the method to use for calculating elevation statistics. Can be `"mean"` or `"median"` (default is `"median"`).
 #'
-#' @return A data frame with the median elevation for each country in `iso_codes_mosaic`.
-#'
-#' @details The function loads the shapefiles for all countries in `iso_codes_mosaic` from the specified directory, generates `n_points` random points within each country's boundary, and queries the Open-Meteo API for the elevation data at each point. It returns the median elevation for each country and saves the data as a CSV file.
+#' @return A data frame with ISO3 codes, country names, and the calculated mean or median elevation for each country.
 #'
+#' @importFrom exactextractr exact_extract
+#' @importFrom raster raster
+#' @importFrom sf st_read
+#' @importFrom dplyr bind_rows
+#' @importFrom utils write.csv
+#' @importFrom glue glue
 #' @examples
 #' \dontrun{
-#' # Example usage to get median elevation data for all countries in iso_codes_mosaic
-#' n_points <- 10
-#' api_key <- "your_api_key_here"
+#' # Define the ISO3 country codes
+#' iso_codes <- c("ZAF", "KEN", "NGA")
+#'
+#' # Get paths for data storage
 #' PATHS <- get_paths()
-#' median_elevations <- get_elevation(n_points, api_key, PATHS)
-#' print(median_elevations)
-#'}
-#' @importFrom httr GET content status_code
-#' @importFrom jsonlite fromJSON
-#' @importFrom utils download.file write.csv txtProgressBar setTxtProgressBar
-#' @importFrom sf st_read st_sample st_coordinates
+#'
+#' # Calculate mean elevation for the countries
+#' mean_elevations <- get_elevation(iso_codes, PATHS, method = "mean")
+#' print(mean_elevations)
+#' }
 #' @export
 
-get_elevation <- function(PATHS, n_points, api_key = NULL) {
+get_elevation <- function(PATHS, iso_codes, method = "median") {
 
-     requireNamespace('sf')
-     requireNamespace('utils')
-
-     # Load the list of ISO3 country codes for the MOSAIC model
-     iso_codes <- MOSAIC::iso_codes_africa
-
-     # Ensure the elevation directory exists
-     if (!dir.exists(PATHS$DATA_ELEVATION)) dir.create(PATHS$DATA_ELEVATION, recursive = TRUE)
+     # Ensure the method is either "mean" or "median"
+     if (!method %in% c("mean", "median")) {
+          stop("Invalid method. Choose 'mean' or 'median'.")
+     }
 
+     # Initialize an empty list to store results
      results_list <- list()
 
-     pb <- utils::txtProgressBar(min = 0, max = length(iso_codes), style = 3)
+     # Loop through each ISO3 country code
+     for (i in 1:length(iso_codes)) {
 
-     for (i in seq_along(iso_codes)) {
+          iso_code <- iso_codes[i]
 
-          country_results_list <- list()
+          message(glue::glue("Calculating {method} elevation for {iso_code}..."))
 
-          # Construct the path to the shapefile for each country
-          shapefile_path <- file.path(PATHS$DATA_SHAPEFILES, paste0(iso_codes[i], "_ADM0.shp"))
+          # Load the shapefile for the country
+          shapefile_path <- file.path(PATHS$DATA_SHAPEFILES, paste0(iso_code, "_ADM0.shp"))
 
-          # Check if the shapefile exists
           if (!file.exists(shapefile_path)) {
-               message(paste("Shapefile for", iso_codes[i], "not found. Skipping."))
+               message(glue::glue("Shapefile for {iso_code} not found. Skipping."))
                next
           }
 
-          # Load the country shapefile
           country_shapefile <- sf::st_read(shapefile_path, quiet = TRUE)
 
-          # Generate n_points random points within the country's shapefile boundary
-          pts <- MOSAIC::generate_country_grid_n(country_shapefile, n_points)
-          coords <- sf::st_coordinates(pts)
-          rm(pts)
-          rm(country_shapefile)
+          # Load the DEM raster
+          dem_filename <- file.path(PATHS$DATA_DEM, glue::glue("{iso_code}_1km_DEM.tif"))
 
-          for (j in 1:nrow(coords)) {
-
-               lat <- coords[j, "Y"]
-               lon <- coords[j, "X"]
-
-               url <- paste0(
-                    "https://customer-api.open-meteo.com/v1/elevation?",
-                    "latitude=", lat,
-                    "&longitude=", lon,
-                    "&apikey=", api_key
-               )
+          if (!file.exists(dem_filename)) {
+               message(glue::glue("DEM file for {iso_code} not found. Skipping."))
+               next
+          }
 
-               temp_file <- tempfile(fileext = ".json")
-               utils::download.file(url, temp_file, mode = "wb", quiet = TRUE)
-               data <- jsonlite::fromJSON(temp_file)
+          dem_raster <- raster::raster(dem_filename)
 
-               country_df <- data.frame(country = iso_codes[i],
-                                        latitude = lat,
-                                        longitude = lon,
-                                        elevation = data$elevation)
+          # Use exactextractr::exact_extract() for faster extraction within the country boundary
+          extracted_values <- exactextractr::exact_extract(dem_raster, country_shapefile)[[1]]
 
-               country_results_list <- c(country_results_list, list(country_df))
+          if (method == 'mean') {
+               extracted_stat <- mean(extracted_values$value, na.rm=TRUE)
+          } else {
+               extracted_stat <- median(extracted_values$value, na.rm=TRUE)
+          }
 
+          # Check if any values were extracted
+          if (is.null(extracted_values) || length(extracted_values) == 0) {
+               message(glue::glue("No elevation data available for {iso_code}. Skipping."))
+               next
           }
 
-          result_country <- do.call(rbind, country_results_list)
-          result_country <- aggregate(elevation ~ country, data = result_country, FUN = median, na.rm = TRUE)
-          results_list <- c(results_list, list(result_country))
+          # Get the country name
+          country_name <- convert_iso_to_country(iso_code)
 
-          utils::setTxtProgressBar(pb, i)
+          # Store the result
+          result_df <- data.frame(
+               country = country_name,
+               iso_code = iso_code,
+               elevation = extracted_stat
+          )
 
+          results_list <- c(results_list, list(result_df))
      }
 
-     close(pb)
+     # Combine results into a single data frame
+     results_df <- dplyr::bind_rows(results_list)
 
-     out <- do.call(rbind, results_list)
+     # Define the output file path
+     file_path <- file.path(PATHS$DATA_ELEVATION, glue::glue("country_elevation_{method}.csv"))
 
-     file_path <- file.path(PATHS$DATA_ELEVATION, "median_elevation_data.csv")
-     utils::write.csv(out, file = file_path, row.names = FALSE)
-     message(paste("Elevation data saved to:", file_path))
+     # Save the results as a CSV file
+     utils::write.csv(results_df, file = file_path, row.names = FALSE)
+     message(glue::glue("Elevation data saved to: {file_path}"))
 
 }

R/get_paths.R

@@ -42,6 +42,7 @@ get_paths <- function(root=NULL) {
      PATHS$DATA_RAW <- file.path(root, "MOSAIC-data/raw")
      PATHS$DATA_PROCESSED <- file.path(root, "MOSAIC-data/processed")
      PATHS$DATA_SHAPEFILES <- file.path(root, "MOSAIC-data/processed/shapefiles")
+     PATHS$DATA_DEM <- file.path(root, "MOSAIC-data/raw/DEM")
      PATHS$DATA_ELEVATION <- file.path(root, "MOSAIC-data/processed/elevation")
      PATHS$DATA_CLIMATE <- file.path(root, "MOSAIC-data/processed/climate")
      PATHS$DATA_ENSO <- file.path(root, "MOSAIC-data/processed/ENSO")

docs/_pkgdown.yml

@@ -50,6 +50,7 @@ reference:
       - get_climate_historical
       - get_climate_forecast
       - get_climate_future
+      - download_country_DEM
       - get_elevation
 
   - title: "Downloading and Processing Data: Cholera Incidence, Disease Dynamics, and Population Characteristics"
@@ -68,6 +69,7 @@ reference:
 
   - title: "Estimating model quantities"
     contents:
+      - est_seasonal_dynamics
       - est_mobility
       - est_WASH_coverage
       - est_symptomatic_prop