Adding functions and scripts for downloading, extracting, and processing observations, initial conditions, land cover types, ERA5 drivers for anchor sites within NA.

modules/assim.sequential/inst/anchor/anchorSites_data_prep.R

@@ -0,0 +1,74 @@
+library(dplyr)
+library(xts)
+library(PEcAn.all)
+library(purrr)
+library(furrr)
+library(lubridate)
+library(nimble)
+library(ncdf4)
+library(PEcAnAssimSequential)
+library(dplyr)
+library(sp)
+library(raster)
+library(zoo)
+library(ggplot2)
+library(mnormt)
+library(sjmisc)
+library(stringr)
+library(doParallel)
+library(doSNOW)
+library(Kendall)
+library(lgarch)
+library(parallel)
+if (future::supportsMulticore()) {
+  future::plan(future::multicore)
+} else {
+  future::plan(future::multisession)
+}
+
+setwd("/projectnb/dietzelab/dongchen/anchorSites/")
+settings <- PEcAn.settings::read.settings("/projectnb/dietzelab/dongchen/anchorSites/SDA/pecan.xml")
+
+in.path <- "/projectnb/dietzelab/hamzed/ERA5/Data/Ensemble"
+out.path <- "/projectnb/dietzelab/dongchen/anchorSites//ERA5_2012_2021"
+
+# prepare ERA5
+paths <- ERA5_met_process(settings = settings, in.path = in.path, out.path = out.path, write.db = F, write = T)
+
+#prepare observations
+settings <- PEcAn.settings::read.settings("/projectnb/dietzelab/dongchen/anchorSites/SDA/pecan.xml")
+settings$state.data.assimilation$Obs_Prep$outdir <- "/projectnb/dietzelab/dongchen/anchorSites/Obs"
+obs <- PEcAnAssimSequential::SDA_OBS_Assembler(settings)
+
+#prepare LC and PFTs
+LC <- PEcAn.data.remote::MODIS_LC_prep(site_info = site_info, time_points = settings$state.data.assimilation$start.date, qc.filter = T)
+
+#tweak LC to actual pfts
+LC[which(LC[,2] %in% c("Deciduous Broadleaf Forests", 
+                       "Deciduous Needleleaf Forests", 
+                       "Mixed Forests"))] <- "temperate.deciduous.HPDA"
+LC[which(LC[,2] %in% c("Evergreen Broadleaf Forests", 
+                       "Evergreen Needleleaf Forests"))] <- "boreal.coniferous"
+LC[which(LC[,2] %in% c("Closed Shrublands", 
+                       "Open Shrublands", 
+                       "Woody Savannas", 
+                       "Savannas", 
+                       "Grasslands",
+                       "Permanent Wetlands",
+                       "Croplands",
+                       "Cropland/Natural Vegetation Mosaics",
+                       ))] <- "semiarid.grassland_HPDA"
+LC[which(LC[,2] %in% c("Urban and Built-up Lands", 
+                       "Permanent Snow and Ice", 
+                       "Barren",
+                       "Water Bodies"))] <- NA
+#delete settings and LC table with NAs.
+del_ind <- which(is.na(LC[,2]))
+LC <- LC[-del_ind,]
+settings <- settings[-del_ind]
+#write settings and PFT file.
+PEcAn.settings::write.settings(settings, outputfile = "pecan.xml")
+write.csv(LC, file = settings$run$settings.1$inputs$pft.site[[1]], row.names = F)
+# settings <- PEcAn.settings::read.settings("/projectnb/dietzelab/dongchen/anchorSites/SDA/pecan.xml")
+# settings <- PEcAn.settings::prepare.settings(settings)
+

modules/data.land/NAMESPACE

@@ -3,6 +3,7 @@
 export(BADM_IC_process)
 export(Clean_Tucson)
 export(EPA_ecoregion_finder)
+export(IC_ISCN_SOC)
 export(InventoryGrowthFusion)
 export(InventoryGrowthFusionDiagnostics)
 export(Read.IC.info.BADM)
@@ -11,11 +12,13 @@ export(Soilgrids_SoilC_prep)
 export(buildJAGSdata_InventoryRings)
 export(cohort2pool)
 export(dataone_download)
+export(download.SM_CDS)
 export(download_package_rm)
 export(ens_veg_module)
 export(extract.stringCode)
 export(extract_FIA)
 export(extract_NEON_veg)
+export(extract_SM_CDS)
 export(extract_soil_gssurgo)
 export(extract_soil_nc)
 export(extract_veg)

modules/data.land/R/IC_BADM_Utilities.R

@@ -156,7 +156,7 @@ Read.IC.info.BADM <-function(lat, long){
             Site = Gdf$SITE_ID %>% unique(),
             Var = Gdf$VARIABLE[1],
             Date = Date.in,
-            Organ = Organ.in,
+            # Organ = Organ.in,
             AGB = PlantWoodIni,
             soil_organic_carbon_content = SoilIni,
             litter_carbon_content = litterIni
@@ -239,15 +239,19 @@ netcdf.writer.BADM <- function(lat, long, siteid, outdir, ens){
 #' @export
 #'
 BADM_IC_process <- function(settings, dir, overwrite=TRUE){
-
-
-  new.site <-
-    data.frame(
-      id = settings$run$site$id %>% as.numeric(),
-      lat = settings$run$site$lat ,
-      lon = settings$run$site$lon %>% as.numeric()
-    )
-
+  # create site info.
+  new.site <- 
+    settings %>% 
+    purrr::map(~.x[['run']] ) %>% 
+    purrr::map('site')%>% 
+    purrr::map(function(site.list){
+      #conversion from string to number
+      site.list$lat <- as.numeric(site.list$lat)
+      site.list$lon <- as.numeric(site.list$lon)
+      list(id=site.list$id, lat=site.list$lat, lon=site.list$lon)
+    })%>% 
+    dplyr::bind_rows() %>% 
+    as.list()
 
   out.ense <- seq_len(settings$ensemble$size) %>%
     purrr::map(~ netcdf.writer.BADM(new.site$lat,
@@ -273,36 +277,49 @@ BADM_IC_process <- function(settings, dir, overwrite=TRUE){
 #' @return a dataframe with codes corresponding to level1 and level2 codes as two columns
 #' @export
 #'
-EPA_ecoregion_finder <- function(Lat, Lon){
+EPA_ecoregion_finder <- function(Lat, Lon, filenames = c("NA_CEC_Eco_Level1.shp", "NA_CEC_Eco_Level2.shp")){
+  #try eco-region map for just CONUS US.
+  if (is.null(filenames)) {
+    filenames = c("eco-region.json", "eco-regionl2.json")
+  }
   #lat long to spatial point
   U.S.SB.sp <-
     data.frame(Lati = Lat %>% as.numeric(),
                Long = Lon %>% as.numeric())
 
   sp::coordinates(U.S.SB.sp) <- ~ Long + Lati
 
-
   # L1 layer
-  L1 <-
-    sf::read_sf(system.file("extdata","eco-region.json", package = "PEcAn.data.land")) %>%
-    sf::st_set_crs(
-      "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"
-    ) %>%
-    sf::st_transform("+proj=longlat +datum=WGS84")
-  # L2 layer
-  L2 <-
-    sf::read_sf(system.file("extdata","eco-regionl2.json", package = "PEcAn.data.land")) %>%
-    sf::st_set_crs(
-      "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"
-    ) %>%
-    sf::st_transform("+proj=longlat +datum=WGS84")
-
-  sp::proj4string(U.S.SB.sp) <- sp::proj4string(sf::as_Spatial(L1))
-  # finding the code for each site
-  over.out.L1 <- sp::over(U.S.SB.sp, sf::as_Spatial(L1))
-  over.out.L2 <- sp::over(U.S.SB.sp, sf::as_Spatial(L2))
+  if (length(filenames) >= 1) {
+    L1 <-
+      sf::read_sf(system.file("extdata",filenames[1], package = "PEcAn.data.land")) %>%
+      sf::st_set_crs(
+        "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"
+      ) %>%
+      sf::st_transform("+proj=longlat +datum=WGS84")
+    #reproject points based on projection of the eco-region map.
+    sp::proj4string(U.S.SB.sp) <- sp::proj4string(sf::as_Spatial(L1))
+    # finding the code for each site
+    over.out.L1 <- sp::over(U.S.SB.sp, sf::as_Spatial(L1))
+  } else {
+    PEcAn.logger::logger.info("Please provide file name to the Level 1 eco-region map!")
+    return(0)
+  }
 
-  return(data.frame(L1 = over.out.L1$NA_L1CODE, L2 = over.out.L2$NA_L2CODE))
+  # L2 layer
+  if (length(filenames) > 1) {
+    L2 <-
+      sf::read_sf(system.file("extdata",filenames[2], package = "PEcAn.data.land")) %>%
+      sf::st_set_crs(
+        "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"
+      ) %>%
+      sf::st_transform("+proj=longlat +datum=WGS84")
+    # finding the code for each site
+    over.out.L2 <- sp::over(U.S.SB.sp, sf::as_Spatial(L2))
+    return(data.frame(L1 = over.out.L1$NA_L1CODE, L2 = over.out.L2$NA_L2CODE))
+  } else {
+    return(data.frame(L1 = over.out.L1$NA_L1CODE))
+  }
 }
 
 

modules/data.land/R/ISCN_extract.R

@@ -0,0 +1,31 @@
+#' Extract ISCN SOC initial conditions from existing ISCN database.
+#'
+#' @param site_info Bety list of site info including site_id, lon, and lat.
+#' @param ens ensemble number.
+#'
+#' @return A data frame containing sampled SOC, each row represent each site.
+#' @export
+#' 
+#' @examples
+#' @author Dongchen Zhang
+#' @importFrom magrittr %>%
+IC_ISCN_SOC <- function(site_info, ens = 100) {
+  iscn_soc <- PEcAn.data.land::iscn_soc
+  site_eco <- PEcAn.data.land::EPA_ecoregion_finder(site_info$lat, site_info$lon)
+  soc <- iscn_soc %>%
+    dplyr::filter(
+      .data$NA_L2CODE %in% site_eco$L2
+    ) %>%
+    dplyr::select("soc (g cm-2)", "NA_L2CODE") %>%
+    na.omit()
+  ic_sample_soc <- data.frame(matrix(NA, ens, length(site_info$site_id))) %>%
+    `colnames<-`(site_info$site_id)
+  for (i in seq_along(site_eco$L2)) {
+    if (is.na(site_eco$L2[i])) {
+      next
+    } else {
+      ic_sample_soc[,i] <- sample(soc$`soc (g cm-2)`[which(soc$NA_L2CODE == site_eco$L2[i])], ens, replace = T)
+    }
+  }
+  return(ic_sample_soc)
+}

modules/data.land/R/download.SM_CDS.R

@@ -0,0 +1,91 @@
+#' Download CDS soil moisture data for the SDA workflow.
+#'
+#' @param outfolder physical paths to where the unziped soil moisture files are downloaded.
+#' @param time_points A vector contains each time point within the start and end date.
+#' @param overwrite flag determine if we want to overwrite existing files when downloading.
+#'
+#' @return A vector containing file paths to the downloaded files.
+#' @export
+#' 
+#' @examples
+#' @author Dongchen Zhang
+#' @importFrom magrittr %>%
+download.SM_CDS <- function(outfolder, time.points, overwrite = FALSE) {
+  #load cdsapi from python environment.
+  tryCatch({
+    cdsapi <- reticulate::import("cdsapi")
+  }, error = function(e) {
+    PEcAn.logger::logger.severe(
+      "Failed to load `cdsapi` Python library. ",
+      "Please make sure it is installed to a location accessible to `reticulate`.",
+      "You should be able to install it with the following command: ",
+      "`pip install --user cdsapi`.",
+      "The following error was thrown by `reticulate::import(\"cdsapi\")`: ",
+      conditionMessage(e)
+    )
+  })
+  #check if the token exists for the cdsapi.
+  if (!file.exists(file.path(Sys.getenv("HOME"), ".cdsapirc")))
+    PEcAn.logger::logger.severe(
+      "Please create a `${HOME}/.cdsapirc` file as described here:",
+      "https://cds.climate.copernicus.eu/api-how-to#install-the-cds-api-key."
+    )
+  #grab the client object.
+  tryCatch({
+    cclient <- cdsapi$Client()
+  }, error = function(e) {
+    PEcAn.logger::logger.severe(
+      "The following error was thrown by `cdsapi$Client()`: ",
+      conditionMessage(e)
+    )
+  })
+  #loop over each time point.
+  file.names <- c()
+  #setup progress bar.
+  pb <- utils::txtProgressBar(min = 0, max = length(time.points), style = 3)
+  for (i in seq_along(time.points)) {
+    #name file.
+    fname <- file.path(outfolder, paste('surface_soil_moisture', time.points[i], "nc", sep = "."))
+    fname.zip <- gsub(".nc", ".zip", fname, fixed = T)
+    #add new extracted file into vector.
+    file.names <- c(file.names, fname)
+    #if we have already downloaded this file.
+    if (file.exists(fname) && !overwrite) {
+      PEcAn.logger::logger.warn(glue::glue(
+        "File `{fname}` already exists, and `overwrite` is FALSE. ",
+        "Skipping to next variable."
+      ))
+      next
+    }
+    #prepare file through cds server.
+    while ("try-error" %in% class(try(do_next <- cclient$retrieve(
+      'satellite-soil-moisture',
+      list(
+        'variable'= 'surface_soil_moisture',
+        'type_of_sensor'= 'active',
+        'time_aggregation'= 'day_average',
+        'year'= sprintf("%04d", lubridate::year(time.points[i])),
+        'month'= sprintf("%02d", lubridate::month(time.points[i])),
+        'day'= sprintf("%02d", lubridate::day(time.points[i])),
+        'type_of_record'= 'cdr',
+        'version'= 'v202212'
+      ),
+      'download.zip'
+    )))) {
+      Sys.sleep(10)
+      PEcAn.logger::logger.info("Encounter error! Will try download in 10 seconds.")
+    }
+    #download file to local.
+    utils::download.file(do_next$reply$location, destfile = fname.zip)
+    #unzip file.
+    unzipPath <- utils::unzip(zipfile = fname.zip, exdir = outfolder)
+    #rename unziped file.
+    base::file.rename(unzipPath, fname)
+    #remove zip file.
+    base::file.remove(fname.zip)
+    #update progress bar
+    pbi <- i
+    utils::setTxtProgressBar(pb, pbi)
+  }
+  file.names
+}