inst/doc/appeears_vignette.R

In appeears: Interface to 'AppEEARS' NASA Web Services

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

# load the library
library(ncdf4)
library(terra)
library(sf)
library(appeears)
library(dplyr)
library(ggplot2)
library(patchwork)

# load demo data

r_polygon <- terra::rast(
  file.path(system.file(package = "appeears"),"extdata/polygon/MCD12Q2.006_Greenup_0_doy2010001_aid0001.tif")
  )

r_raster <- terra::rast(
  file.path(system.file(package = "appeears"),"extdata/raster/MCD12Q2.006_Greenup_0_doy2010001_aid0001.tif")
  )

time_series <- read.table(
  file.path(system.file(package = "appeears"),"extdata/time_series/time-series-MCD43A4-061-results.csv"),
  header = TRUE,
  sep = ","
  )

## ----eval = FALSE-------------------------------------------------------------
# library(appeears)
# 
# # set a key to the keychain
# rs_set_key(
#   user = "earth_data_user",
#   password = "XXXXXXXXXXXXXXXXXXXXXX"
#   )
# 
# # you can retrieve the password using
# rs_get_key(user = "earth_data_user")
# 
# # the output should be the key you provided
# # "XXXXXXXXXXXXXXXXXXXXXX"

## ----eval = FALSE-------------------------------------------------------------
# # request the current token
# token <- rs_login(user = "earth_data_user")
# 
# # invalidate the current session
# rs_logout(token)

## -----------------------------------------------------------------------------
# list all product information
products <- rs_products()

# print the start of all products with their versions
head(products$ProductAndVersion)

# list all layers for a particular
# product
layers <- rs_layers(
  product = "MCD12Q2.006"
)

head(layers)

## ----eval = FALSE-------------------------------------------------------------
# # Load the library
# library(appeears)
# 
# # list all products
# rs_products()
# 
# # list layers of the MOD11A2.061 product
# rs_layers("MOD11A2.061")
# 
# df <- data.frame(
#   task = "time_series",
#   subtask = "US-Ha1",
#   latitude = 42.5378,
#   longitude = -72.1715,
#   start = "2010-01-01",
#   end = "2010-12-31",
#   product = "MCD43A4.061",
#   layer = c("Nadir_Reflectance_Band3","Nadir_Reflectance_Band4")
# )
# 
# # build the area based request/task
# # rename the task name so data will
# # be saved in the "point" folder
# # as defined by the task name
# df$task <- "point"
# task <- rs_build_task(df = df)
# 
# # request the task to be executed
# rs_request(
#   request = task,
#   user = "earth_data_user",
#   transfer = TRUE,
#   path = "~/some_path",
#   verbose = TRUE
# )

## ----eval = FALSE-------------------------------------------------------------
# # read in data
# time_series <- read.table(
#   "~/some_path/time_series/time-series-MCD43A4-061-results.csv",
#   header = TRUE,
#   sep = ","
#   )

## ----warning=FALSE------------------------------------------------------------
# convert band 3 and 4 to NDVI
time_series <- time_series |>
  mutate(
    Date = as.Date(Date),
    NDVI = (MCD43A4_061_Nadir_Reflectance_Band4 - MCD43A4_061_Nadir_Reflectance_Band3)/
      (MCD43A4_061_Nadir_Reflectance_Band4 + MCD43A4_061_Nadir_Reflectance_Band3)
  )

# screen for quality control
time_series <- time_series |>
  mutate(
    NDVI = ifelse(MCD43A4_061_BRDF_Albedo_Band_Mandatory_Quality_Band4 == 255 |
                    MCD43A4_061_BRDF_Albedo_Band_Mandatory_Quality_Band3 == 255,
                  NA, NDVI)
  )

# plot the time series
ggplot(time_series) +
  geom_point(
    aes(
      Date,
      NDVI
    )
  ) +
  theme_bw()


## -----------------------------------------------------------------------------
# load the required libraries
library(appeears)
library(sf)
library(dplyr)
library(ggplot2)

df <- data.frame(
  task = "time_series",
  subtask = "subtask",
  latitude = 42.5378,
  longitude = -72.1715,
  start = "2010-01-01",
  end = "2010-12-31",
  product = "MCD12Q2.006",
  layer = c("Greenup")
)

# load the north carolina demo data
# included in the {sf} package
# and only retain Camden county
roi <- st_read(system.file("gpkg/nc.gpkg", package="sf"), quiet = TRUE) |>
  filter(
    NAME == "Camden"
  )

## ----eval = FALSE-------------------------------------------------------------
# # build the area based request/task
# # rename the task name so data will
# # be saved in the "polygon" folder
# # as defined by the task name
# df$task <- "polygon"
# task <- rs_build_task(
#   df = df,
#   roi = roi,
#   format = "geotiff"
# )
# 
# # request the task to be executed
# rs_request(
#   request = task,
#   user = "earth_data_user",
#   transfer = TRUE,
#   path = "~/some_path",
#   verbose = TRUE
# )

## ----eval=FALSE---------------------------------------------------------------
# library(terra)
# r_polygon <- terra::rast(
#   file.path("~/some_path","polygon/MCD12Q2.006_Greenup_0_doy2010001_aid0001.tif")
#   )

## -----------------------------------------------------------------------------
# convert to data frame for plotting
# with ggplot2, otherwise use the
# tidyterra package and geom_spatrast()
df_polygon <- r_polygon |>
  as.data.frame(xy=TRUE)

# convert incremental values (days since Jan 1, 1970
# to DOY)
df_polygon <- df_polygon |>
  mutate(
    DOY =  as.numeric(
      format(
        as.Date("1970-01-01") + MCD12Q2.006_Greenup_0_doy2010001_aid0001, "%j"
        )
      )
  ) |>
  filter(
    DOY < 180
  )

head(df_polygon)

ggplot() +
  geom_raster(
    data = df_polygon,
    aes(
      x,
      y,
      fill = DOY
    )
  ) +
  labs(
    x = "",
    y = ""
  ) +
  scale_fill_viridis_c() +
  geom_sf(
    data = roi,
    fill = NA,
    colour = "red",
    lwd = 2
    ) +
  theme_bw()
  

## -----------------------------------------------------------------------------
# load the required libraries
library(terra)
library(ggplot2)
library(patchwork)

# create a SpatRaster ROI from the terra demo file
f <- system.file("ex/elev.tif", package="terra")
roi <- terra::rast(f)

## ----eval = FALSE-------------------------------------------------------------
# 
# # build the area based request/task
# # rename the task name so data will
# # be saved in the "raster" folder
# # as defined by the task name
# df$task <- "raster"
# task <- rs_build_task(
#   df = df,
#   roi = roi,
#   format = "geotiff"
# )
# 
# # request the task to be executed
# rs_request(
#   request = task,
#   user = "earth_data_user",
#   transfer = TRUE,
#   path = "~/some_path",
#   verbose = TRUE
# )

## ----eval=FALSE---------------------------------------------------------------
# r_raster <- terra::rast(
#   file.path("~/some_path","raster/MCD12Q2.006_Greenup_0_doy2010001_aid0001.tif")
#   )

## -----------------------------------------------------------------------------
# convert to data frame for plotting
# with ggplot2, otherwise use the
# tidyterra package and geom_spatrast()

df_raster <- r_raster |>
  as.data.frame(xy=TRUE)

# convert incremental values (days since Jan 1, 1970
# to DOY)
df_raster <- df_raster |>
  mutate(
    DOY =  as.numeric(
      format(
        as.Date("1970-01-01") + MCD12Q2.006_Greenup_0_doy2010001_aid0001, "%j"
        )
      )
  ) |>
  filter(
    DOY < 180
  )

df_source <- roi |>
  as.data.frame(xy=TRUE)

p <- ggplot() +
  geom_raster(
    data = df_source,
    aes(
      x,
      y,
      fill = elevation
    )
  ) +
  labs(
    x = "",
    y = ""
  ) +
  scale_fill_viridis_c() +
  theme_bw()

p2 <- ggplot() +
  geom_raster(
    data = df_raster,
    aes(
      x,
      y,
      fill = DOY
    )
  ) +
  labs(
    x = "",
    y = ""
  ) +
  scale_fill_viridis_c() +
  theme_bw()

# patchwork side by side plot
p | p2