diff --git a/.gitignore b/.gitignore
index 38d97cdb..2b89cb3f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -74,4 +74,10 @@ code/mitchell_tests/
 
 # Insang's negative value exploration script
 tools/negative_exploration
-tools/metalearner_test
\ No newline at end of file
+tools/metalearner_test
+
+# NASA Earthdata login credentials
+.dodsrc
+
+# MODIS raw data
+input/data/modis
\ No newline at end of file
diff --git a/input/Rinput/download_functions/download_modis_data.R b/input/Rinput/download_functions/download_modis_data.R
index 214d7e3a..76cdc939 100644
--- a/input/Rinput/download_functions/download_modis_data.R
+++ b/input/Rinput/download_functions/download_modis_data.R
@@ -1,37 +1,55 @@
 ################################################################################
 # Date created: 2023-10-12
-# Packages required: stringr
+# Date modified: 2023-11-27
 ################################################################################
 
 ################################################################################
 #' download_modis_data:
-#' @description
+#' @description Need maintenance for the directory path change
+#' in NASA EOSDIS. This function first retrieves the all hdf download links
+#' on a certain day, then only selects the relevant tiles from the retrieved
+#' links. Download is only done at the queried horizontal-vertical tile number
+#' combinations. An exception is MOD06_L2 product, which is produced
+#' every five minutes every day.
+#' @note \code{date_start} and \code{date_end} should be in the same year.
+#'  Directory structure looks like
+#'  input/modis/raw/[version]/[product]/[year]/[day_of_year]
 #' @param date_start character(1). length of 10. Start date for downloading
 #' data. Format YYYY-MM-DD (ex. September 1, 2023 = "2023-09-01").
 #' @param date_end character(1). length of 10. End date for downloading data.
 #' Format YYYY-MM-DD (ex. September 1, 2023 = "2023-09-01").
-#' @param product character(1).
-#' @param version character(1).
-#' @param horizontal_tiles integer(2).
-#' @param vertical_tiles integer(2).
+#' @param product character(1). One of c("MOD09GA", "MOD11A1", "MOD06_L2",
+#'      "MCD19A2", "MOD13A2", "VNP46A2").
+#' @param version character(1). Default is "61", meaning v061.
+#' @param horizontal_tiles integer(2). Horizontal tile numbers
+#' c([start], [end]). Default is c(7, 13).
+#' @param vertical_tiles integer(2). Vertical tile numbers
+#' c([start], [end]). Default is c(3, 6).
 #' @param nasa_earth_data_token character(1).
+#'  Token for downloading data from NASA. Should be set before
+#'  trying running the function.
 #' @param directory_to_save character(1). Directory to save data.
 #' @param data_download_acknowledgement logical(1). By setting `= TRUE` the
 #' user acknowledge that the data downloaded using this function may be very
 #' large and use lots of machine storage and memory.
-#' @author Mitchell Manware
-#' @return NULL;
+#' @param write_command_only logical(1). Only return a text file
+#' with wget commands instead of running it
+#' @author Mitchell Manware, Insang Song
+#' @import rvest
+#' @returns NULL;
 #' @export
 download_modis_data <- function(
     date_start = "2023-09-01",
     date_end = "2023-09-01",
-    product = NULL,
+    product = c("MOD09GA", "MOD11A1", "MOD06_L2",
+                "MCD19A2", "MOD13A2", "VNP46A2"),
     version = "61",
     horizontal_tiles = c(7, 13),
     vertical_tiles = c(3, 6),
     nasa_earth_data_token = NULL,
     directory_to_save = "./input/modis/raw/",
-    data_download_acknowledgement = FALSE) {
+    data_download_acknowledgement = FALSE,
+    write_command_only = FALSE) {
   #### 1. directory setup
   chars_dir_save <- nchar(directory_to_save)
   if (substr(directory_to_save, chars_dir_save, chars_dir_save) != "/") {
@@ -39,71 +57,147 @@ download_modis_data <- function(
   }
   #### 2. check for data download acknowledgement
   if (data_download_acknowledgement == FALSE) {
-    cat(paste0(
-      "Data download acknowledgement is set to FALSE. ",
-      "Please acknowledge that the data downloaded using this ",
-      "function may be very large and use lots of machine storage ",
-      "and memory.\n"
-    ))
-    stop()
+    stop(
+         paste0(
+                "Data download acknowledgement is set to FALSE. ",
+                "Please acknowledge that the data downloaded using this ",
+                "function may be very large and use lots of machine storage ",
+                "and memory.\n"))
   }
   #### 3. check for NASA earth data token
-  if (nasa_earth_data_token == FALSE) {
-    cat(paste0("Please provide NASA EarthData Login token.\n"))
-    stop()
+  if (is.null(nasa_earth_data_token)) {
+    stop("Please provide NASA EarthData Login token.\n")
   }
   #### 4. check for product
-  if (is.null(product) == TRUE) {
-    cat(paste0("Please select a MODIS product.\n"))
-    stop()
+  product <- match.arg(product)
+  ismod13 <- startsWith(product, "MOD13")
+  ismod06 <- startsWith(product, "MOD06")
+
+  if (substr(date_start, 1, 4) != substr(date_end, 1, 4)) {
+    stop("date_start and date_end should be in the same year.\n")
   }
+
   #### 5. check for version
   if (is.null(version) == TRUE) {
-    cat(paste0("Please select a data version.\n"))
-    stop()
+    stop("Please select a data version.\n")
   }
   #### 6. check for valid horizontal tiles
-  for (h in seq_along(horizontal_tiles)) {
-    if (horizontal_tiles[h] < 0 || horizontal_tiles[h] > 35) {
-      cat(paste0("Horizontal tiles invalid.\n"))
-      stop()
-    }
-  }
-  #### 7. check for valid vertical tiles
-  for (v in seq_along(vertical_tiles)) {
-    if (vertical_tiles[v] < 0 || vertical_tiles[v] > 17) {
-      cat(paste0("Vertical tiles invalid.\n"))
-      stop()
-    }
-  }
-  #### 8. define date sequence
+  # revised to reduce cyclomatic complexity
+  stopifnot(all(horizontal_tiles %in% seq(0, 35)))
+  stopifnot(all(vertical_tiles %in% seq(0, 17)))
+
+  #### 8. Reuse ladsweb home url
+  ladsurl <- "https://ladsweb.modaps.eosdis.nasa.gov/"
+  version <- ifelse(startsWith(product, "VNP"), "5000", version)
+
+  #### 9. define date sequence
   date_start_date_format <- as.Date(date_start, format = "%Y-%m-%d")
   date_end_date_format <- as.Date(date_end, format = "%Y-%m-%d")
   date_sequence <- seq(date_start_date_format, date_end_date_format, "day")
-  #### 9. initiate "..._wget_commands.txt" file
+  if (ismod13) {
+    date_start_yearday <- as.numeric(strftime(date_start_date_format, "%j"))
+    date_end_yearday <- as.numeric(strftime(date_end_date_format, "%j"))
+    year_mod13 <- strftime(date_start_date_format, "%Y")
+    date_sequence <- seq(1, 366, 16)
+    date_sequence <-
+      date_sequence[date_sequence >= date_start_yearday &
+        date_sequence <= date_end_yearday]
+  }
+
+  # In a certain year, list all available dates
+  year <- ifelse(ismod13, year_mod13, as.character(substr(date_start, 1, 4)))
+  filedir_year_url <-
+    paste0(
+          ladsurl,
+          "archive/allData/",
+          version,
+          "/",
+          product,
+          "/",
+          year)
+
+  list_available_d <-
+    rvest::read_html(filedir_year_url) |>
+    rvest::html_elements("tr") |>
+    rvest::html_attr("data-name")
+  # no conditional assignment at this moment.
+  # 11/28/2023 I. Song
+  # remove NAs
+  date_sequence <- list_available_d[!is.na(list_available_d)]
+
+  #### 10. define horizontal tiles
+  tiles_horizontal <- seq(horizontal_tiles[1],
+    horizontal_tiles[2],
+    by = 1
+  )
+  tiles_horizontal <- paste0(
+    "h",
+    sprintf("%02d", tiles_horizontal)
+  )
+
+  #### 11. define vertical tiles
+  tiles_vertical <- seq(vertical_tiles[1],
+    vertical_tiles[2],
+    by = 1
+  )
+  tiles_vertical <- paste0(
+    "v",
+    sprintf("%02d", tiles_vertical)
+  )
+  #### 12. define requested tiles
+  tiles_df <- expand.grid(
+    h = tiles_horizontal,
+    v = tiles_vertical
+  )
+  tiles_requested <-
+    paste0(tiles_df$h, tiles_df$v)
+
+
+  #### 13. initiate "..._wget_commands.txt" file
   commands_txt <- paste0(
     directory_to_save,
     product,
+    "_",
+    date_start,
+    "_",
+    date_end,
+    "_",
+    product,
     "_wget_commands.txt"
   )
+
+  # avoid any possible errors by removing existing command files
+  if (file.exists(commands_txt)) {
+    unlink(commands_txt)
+  }
+
   sink(commands_txt)
-  #### 10. concatenate and print download commands to "..._wget_commands.txt"
+
+  #### 14. append download commands to text file
   for (d in seq_along(date_sequence)) {
     date <- date_sequence[d]
-    year <- as.character(substr(date, 1, 4))
-    day <- strftime(date, "%j")
-    download_url <- paste0(
-      "https://ladsweb.modaps.eosdis.nasa.gov/",
-      "archive/allData/",
-      version,
-      "/",
-      product,
-      "/",
-      year,
-      "/",
-      day,
-      "/"
-    )
+    day <- date
+    filedir_url <-
+      paste0(
+             filedir_year_url,
+             "/",
+             day)
+
+    filelist <-
+      rvest::read_html(filedir_url) |>
+      rvest::html_elements("tr") |>
+      rvest::html_attr("data-path")
+
+    filelist_sub <-
+      grep(
+           paste("(", paste(tiles_requested, collapse = "|"), ")"),
+           filelist, value = TRUE)
+    if (ismod06) {
+      filelist_sub <- filelist
+    }
+    download_url <- sprintf("%s%s", ladsurl, filelist_sub)
+
+    # Main wget run
     download_command <- paste0(
       "wget -e robots=off -m -np -R .html,.tmp ",
       "-nH --cut-dirs=3 \"",
@@ -114,84 +208,24 @@ download_modis_data <- function(
       directory_to_save,
       "\n"
     )
+    #### 15. concatenate and print download commands to "..._wget_commands.txt"
     cat(download_command)
   }
-  #### 11. finish "..._wget_commands.txt"
-  sink()
-  #### 12. build system command
-  system_command <- paste0(
-    ". ",
-    commands_txt,
-    "\n"
-  )
-  #### 13. download data
-  system(command = system_command)
-  #### 14. remove "..._wget_commands.txt"
-  file.remove(commands_txt)
-  #### 15. define horizontal tiles
-  tiles_horizontal <- seq(horizontal_tiles[1],
-    horizontal_tiles[2],
-    by = 1
-  )
-  tiles_horizontal <- paste0(
-    "h",
-    stringr::str_pad(tiles_horizontal,
-      side = "left",
-      width = 2,
-      pad = "0"
-    )
-  )
-  #### 16. define vertical tiles
-  tiles_vertical <- seq(vertical_tiles[1],
-    vertical_tiles[2],
-    by = 1
-  )
-  tiles_vertical <- paste0(
-    "v",
-    stringr::str_pad(tiles_vertical,
-      wide = "left",
-      width = 2,
-      pad = "0"
-    )
-  )
-  #### 17. define requested tiles
-  tiles_requested <- as.vector(NULL)
-  for (t in seq_along(tiles_horizontal)) {
-    tile <- paste0(
-      tiles_horizontal[t],
-      tiles_vertical
-    )
-    tiles_requested <- c(tiles_requested, tile)
-  }
-  #### 18. remove data outside of requested tiles
-  for (s in seq_along(date_sequence)) {
-    date <- date_sequence[s]
-    year <- as.character(substr(date, 1, 4))
-    day <- strftime(date, "%j")
-    directory_with_data <- paste0(
-      directory_to_save,
-      product,
-      "/",
-      year,
-      "/",
-      day,
-      "/"
-    )
-    data_paths <- list.files(
-      path = directory_with_data,
-      full.names = TRUE
-    )
-    path_splitter <- paste0(
-      "A",
-      year,
-      day
+
+  #### 16. finish "..._wget_commands.txt"
+  sink(file = NULL)
+
+  if (!write_command_only) {
+    #### 17. build system command
+    system_command <- paste0(
+      ". ",
+      commands_txt,
+      "\n"
     )
-    for (p in seq_along(data_paths)) {
-      path_tiles <- as.data.frame(strsplit(data_paths[p], path_splitter))[2, ]
-      path_tiles <- substr(path_tiles, 2, 7)
-      if (!(path_tiles %in% tiles_requested)) {
-        file.remove(data_paths[p])
-      }
-    }
+    #### 18. download data
+    system(command = system_command)
   }
+
+  message("Requests were processed.\n")
+
 }
diff --git a/input/Rinput/filter_minimum_poc.R b/input/Rinput/processing_functions/filter_minimum_poc.R
similarity index 100%
rename from input/Rinput/filter_minimum_poc.R
rename to input/Rinput/processing_functions/filter_minimum_poc.R
diff --git a/input/Rinput/processing_functions/filter_unique_sites.R b/input/Rinput/processing_functions/filter_unique_sites.R
new file mode 100644
index 00000000..3c0f5a9c
--- /dev/null
+++ b/input/Rinput/processing_functions/filter_unique_sites.R
@@ -0,0 +1,86 @@
+## Date: 2023-11-24
+## Description: Get the unique site locations and the WGS84 coordinates of those
+check_installed_load <- function(lib) {
+  if (!require(lib, character.only = TRUE)) {
+    install.packages(lib)
+    library(lib, character.only = TRUE)
+  }
+}
+
+pkgs <- c("data.table", "sf", "terra")
+invisible(suppressMessages(sapply(pkgs, check_installed_load)))
+options(sf_use_s2 = FALSE)
+
+
+
+#' Filter unique sites with or without temporal information
+#' @param path_measurement character(1). Path to daily measurement data.
+#' @param include_time logical(1). Should the output include
+#'  temporal information?
+#' @param date_start character(1). Start date.
+#'  Should be in "YYYY-MM-DD" format.
+#' @param date_end character(1). End date.
+#'  Should be in "YYYY-MM-DD" format.
+#' @returns data.table object with three or four fields.
+#' - "site_id"
+#' - "lon": in WGS 1984 (EPSG:4326)
+#' - "lat": in WGS 1984 (EPSG:4326)
+#' - "date"
+#' "date" field will be attached only when include_time == TRUE.
+#' @import data.table
+#' @note \code{include_time = TRUE} will return a massive data.table
+#' object. Please choose proper \code{date_start} and \code{date_end} values.
+#' @export
+filter_unique_sites <-
+  function(
+    path_measurement = "./tests/testdata/daily_88101_2018-2022.rds",
+    include_time = FALSE,
+    date_start = "2018-01-01",
+    date_end = "2022-12-31"
+  ) {
+    sites <- readRDS(path_measurement)
+    # data manipulation
+    source("./R/manipulate_spacetime_data.R")
+
+    ## get unique sites
+    ##
+    sites$site_id <-
+      sprintf("%02d%03d%04d%05d",
+              sites$State.Code,
+              sites$County.Code,
+              sites$Site.Num,
+              sites$Parameter.Code)
+
+    # select relevant fields only
+    sites_v <- unique(sites[, c("site_id", "Longitude", "Latitude", "Datum")])
+    names(sites_v)[2:3] <- c("lon", "lat")
+    sites_v <- as.data.table(sites_v)
+
+    # subset mainland
+    sites_v <- sites_v[!grepl("^(02|15|72|78|6)", site_id), ]
+
+    # NAD83 to WGS84
+    sites_v_nad <-
+      project_dt(sites_v[Datum == "NAD83"],
+                 "EPSG:4269", "EPSG:4326")
+    sites_v_nad <- sites_v_nad[, c(3, 6, 5)]
+    sites_v_wgs <- sites_v[Datum == "WGS84"][, -4]
+    final_sites <- rbind(sites_v_wgs, sites_v_nad)
+
+    if (include_time) {
+      date_start <- as.Date(date_start)
+      date_end <- as.Date(date_end)
+      date_sequence <- seq(date_start, date_end, "day")
+      final_sites <-
+        split(date_sequence, date_sequence) |>
+        lapply(function(x) {
+          fs_time <- final_sites
+          fs_time$date <- x
+          return(fs_time)
+        })
+      final_sites <- Reduce(rbind, final_sites)
+    }
+
+    return(final_sites)
+}
+# File ends
\ No newline at end of file
diff --git a/input/Rinput/processing_functions/process_modis.R b/input/Rinput/processing_functions/process_modis.R
new file mode 100644
index 00000000..95a13824
--- /dev/null
+++ b/input/Rinput/processing_functions/process_modis.R
@@ -0,0 +1,144 @@
+# description: process modis data with
+# buffer radii and the list of target covariates
+# Insang Song
+# Updated: 11/24/2023
+# set base directory to ddn (provided that users ssh at triton)
+rootdir <- "/ddn/gs1"
+basedir <- "group/set/Projects/NRT-AP-model/"
+userdir <- "home/songi2/projects/NRTAPModel"
+inputdir <- "input/data"
+modisdir <- "modis/raw"
+rlibdir <- "home/songi2/r-libs/"
+download_dir <- file.path(rootdir, basedir, inputdir, modisdir)
+download_dir_user <- file.path(rootdir, userdir, inputdir, modisdir)
+lib_dir <- file.path(rootdir, rlibdir)
+
+if (!dir.exists(download_dir)) {
+  dir.create(download_dir)
+}
+if (!dir.exists(download_dir_user)) {
+  dir.create(download_dir_user)
+}
+
+# register custom library path for debugging on Triton
+.libPaths(lib_dir)
+
+# filter unique sites
+source("./input/Rinput/processing_functions/filter_unique_sites.R")
+
+unique_sites <- filter_unique_sites()
+unique_sites_t <- filter_unique_sites(include_time = TRUE)
+
+
+# check if remotes package is available
+if (!require(remotes)) {
+  install.packages("remotes")
+  library(remotes)
+}
+
+# install scalable_gis
+if (!require(scomps)) {
+  remotes::install_github("Spatiotemporal-Exposures-and-Toxicology/Scalable_GIS")
+}
+
+pkgs <- c("terra", "sf", "future.apply", "scomps", "exactextractr", "foreach", "data.table", "tigris")
+suppressMessages(invisible(sapply(pkgs, library, character.only = TRUE, quietly = TRUE)))
+tigris_cache_dir("~/tigris_cache")
+options(sf_use_s2 = FALSE, tigris_use_cache = TRUE)
+# rr <- terra::rast("...hdf")
+# summary(rr["1 km 16 days NDVI"])
+
+# https://opendap.cr.usgs.gov/opendap/hyrax/MOD09GA.061/h10v07.ncml.dap.nc4?dap4.ce=/MODIS_Grid_1km_2D_eos_cf_projection;/Latitude_1[0:1:2399][0:1:2399];/Longitude_1[0:1:2399][0:1:2399];/MODIS_Grid_500m_2D_eos_cf_projection;/sur_refl_b03_1[6574:1:6574][0:1:2399][0:1:2399];/sur_refl_b07_1[6574:1:6574][0:1:2399][0:1:2399];/sur_refl_b02_1[6574:1:6574][0:1:2399][0:1:2399];/sur_refl_b05_1[6574:1:6574][0:1:2399][0:1:2399];/sur_refl_b04_1[6574:1:6574][0:1:2399][0:1:2399];/sur_refl_b06_1[6574:1:6574][0:1:2399][0:1:2399];/sur_refl_b01_1[6574:1:6574][0:1:2399][0:1:2399];/time[6574:1:6574]
+# https://opendap.cr.usgs.gov/opendap/hyrax/MOD09GA.061/h10v07.ncml.dap.nc?dap4.ce=/MODIS_Grid_1km_2D_eos_cf_projection
+# https://opendap.cr.usgs.gov/opendap/hyrax/MOD09GA.061/h10v07.ncml.dap.nc4?dap4.ce=/Latitude_1%5B0:1:2399%5D%5B0:1:2399%5D;/Longitude_1%5B0:1:2399%5D%5B0:1:2399%5D;/MODIS_Grid_500m_2D_eos_cf_projection;/sur_refl_b03_1%5B6574:1:6574%5D%5B0:1:2399%5D%5B0:1:2399%5D;/sur_refl_b07_1%5B6574:1:6574%5D%5B0:1:2399%5D%5B0:1:2399%5D;/sur_refl_b02_1%5B6574:1:6574%5D%5B0:1:2399%5D%5B0:1:2399%5D;/sur_refl_b05_1%5B6574:1:6574%5D%5B0:1:2399%5D%5B0:1:2399%5D;/sur_refl_b04_1%5B6574:1:6574%5D%5B0:1:2399%5D%5B0:1:2399%5D;/sur_refl_b06_1%5B6574:1:6574%5D%5B0:1:2399%5D%5B0:1:2399%5D;/sur_refl_b01_1%5B6574:1:6574%5D%5B0:1:2399%5D%5B0:1:2399%5D;/time%5B6574:1:6574%5D
+## main part
+## hierarchy (arbitrary subsets of census regions)
+states <- tigris::states(year = 2020)
+states <- vect(states)
+states_main <- states[!states$GEOID %in% c("02", "15", "60", "66", "68", "69", "72", "78"), ]
+mod_cenreg <- read.csv("./input/data/regional_divisions.csv")
+mod_cenreg <- mod_cenreg |>
+  tidyr::separate_wider_delim(
+    col = "States",
+    delim = "-",
+    names = paste0("state_", sprintf("%02d", 1:6)),
+    too_few = "align_start") |>
+  tidyr::pivot_longer(cols = 2:7) |>
+  dplyr::filter(!is.na(value)) |>
+  dplyr::select(-name)
+states_m <- merge(states_main, mod_cenreg, by.x = "STUSPS", by.y = "value")
+mod_cr <- terra::aggregate(states_m, by = "Region") |>
+  terra::project("EPSG:5070")
+
+
+## path setting
+modis_mod13 <-
+  list.files(path = "~/Documents/input_test/61/MOD13A2/2018/001/",
+             pattern = "*.hdf$",
+             full.names = TRUE)
+
+# virtual raster on the same date
+# should run terra::describe(..., sds = TRUE)
+# when you do not know about the exact layer to use
+terra::describe(modis_mod13[1], sds = TRUE)
+# $var
+# get names
+names(terra::rast(modis_mod13[1], lyrs = 1))
+
+# mosaic all in daily data
+modis_mod13vrt <- lapply(
+  modis_mod13, \(x) terra::rast(x, lyrs = 1)) |>
+  do.call(what = terra::mosaic, args = _)
+
+
+
+plan(multicore, workers = 6L)
+parallel::makeCluster(6L)
+
+
+sites_v <- final_sites |>
+  terra::vect(crs = "EPSG:4326") |>
+  terra::project("EPSG:5070")
+
+# mod13_10k <- scomps::distribute_process_hierarchy(
+#       regions = mod_cr,
+#       split_level = mod_cr$Region,
+#       fun_dist = scomps::extract_with_buffer,
+#       points = sites_v,
+#       surf = modis_mod13vrt,
+#       radius = 1e4L,
+#       id = "Site.ID",
+#       func = "mean"
+#     )
+
+
+# the data values in MODIS HDF files seem
+# to have 0.00000001 scale, not like 0.00001 in the documentation.
+# strange thing is that the value scale is restored after
+# running gdal_translate or h4toh5.
+# should consult NASA for details.
+scomps::extract_with_buffer(
+  points = sites_v,
+  surf = modis_mod13vrt,
+  radius = 1e4L,
+  id = "Site.ID"
+)
+
+mod13_10k <- 
+  foreach(image = modis_mod13,
+          .combine = dplyr::bind_rows,
+          .export = c("mod_cr", "sites_v"),
+          .packages = c("terra", "scomps", "exactextractr", "foreach", "future", "dplyr")) %do%
+    {
+    modis_in <- terra::rast(image)
+    scomps::distribute_process_hierarchy(
+      regions = mod_cr,
+      split_level = mod_cr$region,
+      fun_dist = scomps::extract_with_buffer,
+      points = final_sites,
+      surf = modis_in,
+      radius = 1e4L,
+      id = "Site.ID",
+      func = "mean"
+    )
+    }
diff --git a/input/Rinput/process_noaa_hms_smoke_data.R b/input/Rinput/processing_functions/process_noaa_hms_smoke_data.R
similarity index 100%
rename from input/Rinput/process_noaa_hms_smoke_data.R
rename to input/Rinput/processing_functions/process_noaa_hms_smoke_data.R
diff --git a/input/Rinput/run_modis_download.R b/input/Rinput/run_modis_download.R
new file mode 100644
index 00000000..93acb7bf
--- /dev/null
+++ b/input/Rinput/run_modis_download.R
@@ -0,0 +1,61 @@
+# assuming that the working directory is the parent directory of the repository
+rootdir <- "/ddn/gs1"
+basedir <- "group/set/Projects/NRT-AP-model/"
+userdir <- "home/songi2/projects/NRTAPModel"
+inputdir <- "input/data"
+modisdir <- "modis/raw"
+rlibdir <- "home/songi2/r-libs/"
+download_dir <- file.path(rootdir, basedir, inputdir, modisdir)
+download_dir_user <- file.path(rootdir, userdir, inputdir, modisdir)
+lib_dir <- file.path(rootdir, rlibdir)
+
+
+# sourcing modis data download function
+source("./input/Rinput/download_functions/download_modis_data.R")
+
+
+products <-
+  c(
+    "MOD09GA",
+    "MOD11A1",
+    "MOD06_L2",
+    "MCD19A2",
+    "MOD13A2",
+    "VNP46A2"
+  )
+# MOD06_L2 are not tiled;
+# VNP46A2 is based on h5 formats
+
+
+# Earthdata token (My Profile > Generate Token)
+# Valid for 60 days; should reauthorize
+edtoken <- readLines("~/.edtoken")[1]
+
+
+args <- commandArgs(trailingOnly = TRUE)
+product_query <- as.character(args[1])
+dstart <- as.character(args[2])
+dend <- as.character(args[3])
+
+# 2023-11-27: The shell argument should describe
+# the full dates of start and end in an appropriate format
+
+#
+# 1826 days
+# MOD09GA 2.0-2.5 GB/day (4.2 TB)
+# MOD11A1 80 MB/day
+# MCD19A2 160 MB/day
+# MOD13A2 300 MB/day
+# VNP46A2 110 MB/day
+# 1 TB for other four;
+download_modis_data(
+  date_start = dstart,
+  date_end = dend,
+  product = product_query,
+  version = "61",
+  horizontal_tiles = c(7, 13),
+  vertical_tiles = c(3, 6),
+  nasa_earth_data_token = edtoken,
+  directory_to_save = download_dir_user,
+  data_download_acknowledgement = TRUE
+)
diff --git a/input/data/regional_divisions.csv b/input/data/regional_divisions.csv
new file mode 100644
index 00000000..e684128e
--- /dev/null
+++ b/input/data/regional_divisions.csv
@@ -0,0 +1,14 @@
+Region,States
+Northwest,WA-OR-MT-ID-WY
+Southwest Pacific,CA-NV
+Southwest Mountain,AZ-NM-CO-UT
+West North Central North,ND-SD-MN
+West North Central South,NE-KS-IA-MO
+East North Central,WI-IL-IN-MI-OH
+West South Central West,OK-TX
+West South Central East,AR-LA-MS
+East South Central,KY-TN-AL
+South Atlantic South,SC-GA-FL
+South Atlantic North,NC-VA-WV-DC-MD-DE
+Middle Atlantic,NY-PA-NJ
+New England,RI-CT-MA-VT-NH-ME