diff --git a/DESCRIPTION b/DESCRIPTION index 467ad0c..52ecc17 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -28,6 +28,7 @@ Imports: Suggests: geodata, ggplot2, + httr2, knitr, tidyterra, testthat (>= 3.0.0) diff --git a/R/blackmarbler.R b/R/blackmarbler.R index 2b07367..2f29ec0 100644 --- a/R/blackmarbler.R +++ b/R/blackmarbler.R @@ -1,43 +1,43 @@ # BlackMarblerR month_start_day_to_month <- function(x){ - + month <- NA - + if(x == "001") month <- "01" - + if(x == "032") month <- "02" - + if(x == "060") month <- "03" if(x == "061") month <- "03" - + if(x == "091") month <- "04" if(x == "092") month <- "04" - + if(x == "121") month <- "05" if(x == "122") month <- "05" - + if(x == "152") month <- "06" if(x == "153") month <- "06" - + if(x == "182") month <- "07" if(x == "183") month <- "07" - + if(x == "213") month <- "08" if(x == "214") month <- "08" - + if(x == "244") month <- "09" if(x == "245") month <- "09" - + if(x == "274") month <- "10" if(x == "275") month <- "10" - + if(x == "305") month <- "11" if(x == "306") month <- "11" - + if(x == "335") month <- "12" if(x == "336") month <- "12" - + return(month) } @@ -53,12 +53,12 @@ pad3 <- Vectorize(pad3) remove_fill_value <- function(x, variable){ # Remove fill values - + # https://viirsland.gsfc.nasa.gov/PDF/BlackMarbleUserGuide_v1.2_20220916.pdf # * Table 3 (page 12) # * Table 6 (page 16) # * Table 9 (page 18) - + #### 255 if(variable %in% c( "Granule", @@ -76,13 +76,13 @@ remove_fill_value <- function(x, variable){ )){ x[][x[] == 255] <- NA } - - #### -999.9 + + #### -999.9 if(variable %in% c("UTC_Time")){ x[][x[] == -999.9] <- NA } - - #### -32768 + + #### -32768 if(variable %in% c("Sensor_Azimuth", "Sensor_Zenith", "Solar_Azimuth", @@ -94,8 +94,8 @@ remove_fill_value <- function(x, variable){ "Moon_Phase_Angle")){ x[][x[] == -32768] <- NA } - - + + #### 65535 if(variable %in% c( "DNB_At_Sensor_Radiance_500m", @@ -138,28 +138,28 @@ remove_fill_value <- function(x, variable){ )){ x[][x[] == 65535] <- NA } - + return(x) } apply_scaling_factor <- function(x, variable){ # Apply scaling factor to variables according to Black Marble user guide - + # https://viirsland.gsfc.nasa.gov/PDF/BlackMarbleUserGuide_v1.2_20220916.pdf # * Table 3 (page 12) # * Table 6 (page 16) # * Table 9 (page 18) - + if(variable %in% c( - + # VNP46A1 "DNB_At_Sensor_Radiance", - + # VNP46A2 "DNB_BRDF-Corrected_NTL", "Gap_Filled_DNB_BRDF-Corrected_NTL", "DNB_Lunar_Irradiance", - + # VNP46A3/4 "AllAngle_Composite_Snow_Covered", "AllAngle_Composite_Snow_Covered_Std", @@ -174,11 +174,11 @@ apply_scaling_factor <- function(x, variable){ "OffNadir_Composite_Snow_Free", "OffNadir_Composite_Snow_Free_Std") ){ - + x <- x * 0.1 - + } - + return(x) } @@ -188,151 +188,151 @@ file_to_raster <- function(h5_file, # Converts h5 file to raster. # ARGS # --f: Filepath to h5 file - + ## Data h5_data <- h5file(h5_file, "r+") - + #### Daily if(h5_file %>% str_detect("VNP46A1|VNP46A2")){ - + tile_i <- h5_file %>% stringr::str_extract("h\\d{2}v\\d{2}") - + bm_tiles_sf <- read_sf("https://raw.githubusercontent.com/worldbank/blackmarbler/main/data/blackmarbletiles.geojson") grid_i_sf <- bm_tiles_sf[bm_tiles_sf$TileID %in% tile_i,] - + grid_i_sf_box <- grid_i_sf %>% st_bbox() - + xMin <- min(grid_i_sf_box$xmin) %>% round() yMin <- min(grid_i_sf_box$ymin) %>% round() xMax <- max(grid_i_sf_box$xmax) %>% round() yMax <- max(grid_i_sf_box$ymax) %>% round() - + var_names <- h5_data[["HDFEOS/GRIDS/VNP_Grid_DNB/Data Fields"]]$names - + if(!(variable %in% var_names)){ warning(paste0("'", variable, "'", " not a valid variable option. Valid options include:\n", paste(var_names, collapse = "\n") )) } - + out <- h5_data[[paste0("HDFEOS/GRIDS/VNP_Grid_DNB/Data Fields/", variable)]][,] qf <- h5_data[["HDFEOS/GRIDS/VNP_Grid_DNB/Data Fields/Mandatory_Quality_Flag"]][,] - + if(length(quality_flag_rm) > 0){ if(variable %in% c("DNB_BRDF-Corrected_NTL", "Gap_Filled_DNB_BRDF-Corrected_NTL", "Latest_High_Quality_Retrieval")){ - + for(val in quality_flag_rm){ # out[qf %in% quality_flag_rm] doesn't work, so loop out[qf == val] <- NA } } } - + # # Above doesn't fully capture # if(variable %in% "Latest_High_Quality_Retrieval"){ # out[out == 255] <- NA # } - + #### Monthly/Annually } else{ - + lat <- h5_data[["HDFEOS/GRIDS/VIIRS_Grid_DNB_2d/Data Fields/lat"]][] lon <- h5_data[["HDFEOS/GRIDS/VIIRS_Grid_DNB_2d/Data Fields/lon"]][] - + var_names <- h5_data[["HDFEOS/GRIDS/VIIRS_Grid_DNB_2d/Data Fields"]]$names - + if(!(variable %in% var_names)){ warning(paste0("'", variable, "'", " not a valid variable option. Valid options include:\n", paste(var_names, collapse = "\n") )) } - + out <- h5_data[[paste0("HDFEOS/GRIDS/VIIRS_Grid_DNB_2d/Data Fields/", variable)]][,] - + if(length(quality_flag_rm) > 0){ - + variable_short <- variable %>% str_replace_all("_Num", "") %>% str_replace_all("_Std", "") - + qf_name <- paste0(variable_short, "_Quality") - + if(qf_name %in% var_names){ - + qf <- h5_data[[paste0("HDFEOS/GRIDS/VIIRS_Grid_DNB_2d/Data Fields/", qf_name)]][,] - + for(val in quality_flag_rm){ # out[qf %in% quality_flag_rm] doesn't work, so loop out[qf == val] <- NA } - + } - + } - + if(class(out[1,1])[1] != "numeric"){ out <- matrix(as.numeric(out), # Convert to numeric matrix ncol = ncol(out)) } - + xMin <- min(lon) %>% round() yMin <- min(lat) %>% round() xMax <- max(lon) %>% round() yMax <- max(lat) %>% round() - + } - + ## Metadata nRows <- nrow(out) nCols <- ncol(out) res <- nRows #nodata_val <- NA myCrs <- "EPSG:4326" - + ## Make Raster - + #transpose data to fix flipped row and column order #depending upon how your data are formatted you might not have to perform this out <- t(out) - + #assign data ignore values to NA #out[out == nodata_val] <- NA - + #turn the out object into a raster outr <- terra::rast(out, crs = myCrs, extent = c(xMin,xMax,yMin,yMax)) - + #set fill values to NA outr <- remove_fill_value(outr, variable) - + #apply scaling factor outr <- apply_scaling_factor(outr, variable) - + #h5closeAll() h5_data$close_all() - + return(outr) } read_bm_csv <- function(year, day, product_id){ - - - # + + + # df_out <- tryCatch( { df <- readr::read_csv(paste0("https://ladsweb.modaps.eosdis.nasa.gov/archive/allData/5000/",product_id,"/",year,"/",day,".csv"), show_col_types = F) - - + + df$year <- year df$day <- day - + df }, error = function(e){ @@ -340,9 +340,9 @@ read_bm_csv <- function(year, data.frame(NULL) } ) - + Sys.sleep(0.1) - + return(df_out) } @@ -351,7 +351,7 @@ create_dataset_name_df <- function(product_id, years = NULL, months = NULL, days = NULL){ - + #### Prep dates if(product_id %in% c("VNP46A1", "VNP46A2")) months <- NULL if(product_id %in% c("VNP46A3")) days <- NULL @@ -359,71 +359,71 @@ create_dataset_name_df <- function(product_id, days <- NULL months <- NULL } - + #### Determine end year year_end <- Sys.Date() %>% substring(1,4) %>% as.numeric() - + #### Make parameter dataframe if(product_id %in% c("VNP46A1", "VNP46A2")){ param_df <- cross_df(list(year = 2012:year_end, day = pad3(1:366))) } - + if(product_id == "VNP46A3"){ param_df <- cross_df(list(year = 2012:year_end, day = c("001", "032", "061", "092", "122", "153", "183", "214", "245", "275", "306", "336", "060", "091", "121", "152", "182", "213", "244", "274", "305", "335"))) } - + if(product_id == "VNP46A4"){ param_df <- cross_df(list(year = 2012:year_end, day = "001")) } - + #### Add month if daily or monthly data if(product_id %in% c("VNP46A1", "VNP46A2", "VNP46A3")){ - + param_df <- param_df %>% dplyr::mutate(month = day %>% month_start_day_to_month() %>% as.numeric()) - + } - + #### Subset time period ## Year if(!is.null(years)){ param_df <- param_df[param_df$year %in% years,] } - + ## Month if(product_id %in% c("VNP46A1", "VNP46A2", "VNP46A3")){ - + if(!is.null(months)){ param_df <- param_df[as.numeric(param_df$month) %in% as.numeric(months),] } - + if(!is.null(days)){ param_df <- param_df[as.numeric(param_df$day) %in% as.numeric(days),] } - + } - + #### Create data # files_df <- purrr::map2_dfr(param_df$year, # param_df$day, # read_bm_csv, # product_id) - + files_df <- purrr::map2(param_df$year, param_df$day, read_bm_csv, product_id) %>% bind_rows() - - + + return(files_df) } @@ -434,97 +434,97 @@ download_raster <- function(file_name, quality_flag_rm, h5_dir, quiet){ - + year <- file_name %>% substring(10,13) day <- file_name %>% substring(14,16) product_id <- file_name %>% substring(1,7) - + url <- paste0('https://ladsweb.modaps.eosdis.nasa.gov/archive/allData/5000/', product_id, '/', year, '/', day, '/', file_name) - + headers <- c('Authorization' = paste('Bearer', bearer)) - + if(is.null(h5_dir)){ download_path <- file.path(temp_dir, file_name) } else{ download_path <- file.path(h5_dir, file_name) } - + if(!file.exists(download_path)){ - + if(quiet == FALSE) message(paste0("Processing: ", file_name)) - + if(quiet == TRUE){ - - response <- httr::GET(url, + + response <- httr::GET(url, httr::timeout(60), - httr::add_headers(headers), + httr::add_headers(headers), httr::write_disk(download_path, overwrite = TRUE)) - + } else{ - response <- httr::GET(url, + response <- httr::GET(url, httr::timeout(60), - httr::add_headers(headers), + httr::add_headers(headers), httr::write_disk(download_path, overwrite = TRUE), httr::progress()) - + } - - - + + + if(response$status_code != 200){ message("Error in downloading data") message(response) } - + if(response$all_headers[[1]]$status != 200){ message("**Error in downloading data; bearer token likely invalid.** Try regenerating the bearer token; please see this link for instructions to obtain a bearer token: https://github.com/worldbank/blackmarbler?tab=readme-ov-file#bearer-token-") } - + } - + r <- file_to_raster(download_path, variable, quality_flag_rm) - + return(r) } define_variable <- function(variable, product_id){ - + if(is.null(variable)){ if(product_id == "VNP46A1") variable <- "DNB_At_Sensor_Radiance_500m" if(product_id == "VNP46A2") variable <- "Gap_Filled_DNB_BRDF-Corrected_NTL" if(product_id %in% c("VNP46A3", "VNP46A4")) variable <- "NearNadir_Composite_Snow_Free" } - + return(variable) } define_date_name <- function(date_i, product_id){ - + #### Make name for raster based on date if(product_id %in% c("VNP46A1", "VNP46A2")){ date_name_i <- paste0("t", date_i %>% str_replace_all("-", "_")) } - + if(product_id %in% c("VNP46A3")){ date_name_i <- paste0("t", date_i %>% str_replace_all("-", "_") %>% substring(1,7)) } - + if(product_id %in% c("VNP46A4")){ date_name_i <- paste0("t", date_i %>% str_replace_all("-", "_") %>% substring(1,4)) } - + return(date_name_i) } count_n_obs <- function(values, coverage_fraction) { ## Function to count observations, for exact_extract - + orig_vars <- names(values) - + values %>% dplyr::mutate(across(orig_vars, ~ as.numeric(!is.na(.)) )) %>% dplyr::summarise(across(orig_vars, sum, .names = "n_non_na_pixels.{.col}"), @@ -627,29 +627,29 @@ bm_extract <- function(roi_sf, h5_dir = NULL, quiet = FALSE, ...){ - + # Errors & Warnings ---------------------------------------------------------- if( (interpol_na == T) & (length(date) == 1) ){ stop("If interpol_na = TRUE, then must have more than one date") } - + if( (interpol_na == T) & (output_location_type == "file") ){ interpol_na <- F warning("interpol_na ignored. Interpolation only occurs when output_location_type = 'memory'") } - + if(class(roi_sf)[1] == "SpatVector") roi_sf <- roi_sf %>% st_as_sf() if(!("sf" %in% class(roi_sf))){ stop("roi must be an sf object") } - - # Required parameters used in try statement, so error not generated when used, + + # Required parameters used in try statement, so error not generated when used, # so use them here roi_sf <- roi_sf product_id <- product_id date <- date bearer <- bearer - + # Assign interpolation variables --------------------------------------------- if(interpol_na == T){ if(!exists("method")) method <- "linear" @@ -659,32 +659,32 @@ bm_extract <- function(roi_sf, if(!exists("z")) z <- NULL if(!exists("NArule")) NArule <- 1 } - + # Define Tempdir ------------------------------------------------------------- temp_main_dir = tempdir() - + current_time_millis = as.character(as.numeric(Sys.time())*1000) %>% str_replace_all("[:punct:]", "") temp_dir = file.path(temp_main_dir, paste0("bm_raster_temp_", current_time_millis)) - + dir.create(temp_dir, showWarnings = F) - + # NTL Variable --------------------------------------------------------------- variable <- define_variable(variable, product_id) - + # Filename root -------------------------------------------------------------- # Define outside of lapply, as use this later to aggregate rasters if(output_location_type == "file"){ - out_name_begin <- paste0(file_prefix, - product_id, "_", + out_name_begin <- paste0(file_prefix, + product_id, "_", variable, "_", - "qflag", + "qflag", quality_flag_rm %>% paste0(collapse="_"), "_", aggregation_fun %>% paste0(collapse="_")) } - + if(interpol_na == T){ - + #### Create raster bm_r <- bm_raster(roi_sf = roi_sf, product_id = product_id, @@ -697,7 +697,7 @@ bm_extract <- function(roi_sf, h5_dir = h5_dir, quiet = quiet, temp_dir = temp_dir) - + bm_r <- terra::approximate(bm_r, method = method, rule = rule, @@ -705,53 +705,53 @@ bm_extract <- function(roi_sf, ties = ties, z = z, NArule = NArule) - + #### Extract roi_df <- roi_sf %>% st_drop_geometry() roi_df$date <- NULL - + n_obs_df <- exact_extract(bm_r, roi_sf, count_n_obs, progress = !quiet) %>% bind_cols(roi_df) %>% tidyr::pivot_longer(cols = -c(names(roi_df)), names_to = c(".value", "date"), names_sep = "\\.t") %>% dplyr::mutate(prop_non_na_pixels = .data$n_non_na_pixels / .data$n_pixels) - + ntl_df <- exact_extract(bm_r, roi_sf, aggregation_fun, progress = !quiet) %>% tidyr::pivot_longer(cols = everything(), names_to = c(".value", "date"), names_sep = "\\.t") - - names(ntl_df)[names(ntl_df) != "date"] <- + + names(ntl_df)[names(ntl_df) != "date"] <- paste0("ntl_", names(ntl_df)[names(ntl_df) != "date"]) - + ntl_df$date <- NULL r <- bind_cols(n_obs_df, ntl_df) - + # Apply through each date, extract, then append } else{ - + # Download data -------------------------------------------------------------- r_list <- lapply(date, function(date_i){ - + #out <- tryCatch( # { - + #### Make name for raster based on date date_name_i <- define_date_name(date_i, product_id) - + #### If save to file if(output_location_type == "file"){ - + out_name_end <- paste0("_", date_name_i, ".Rds") out_name <- paste0(out_name_begin, out_name_end) out_path <- file.path(file_dir, out_name) - + make_raster <- TRUE if(file_skip_if_exists & file.exists(out_path)) make_raster <- FALSE - + if(make_raster){ - + #### Make raster r <- bm_raster_i(roi_sf = roi_sf, product_id = product_id, @@ -764,13 +764,13 @@ bm_extract <- function(roi_sf, quiet = quiet, temp_dir = temp_dir) names(r) <- date_name_i - + #### Extract - r_agg <- exact_extract(x = r, y = roi_sf, fun = aggregation_fun, + r_agg <- exact_extract(x = r, y = roi_sf, fun = aggregation_fun, progress = !quiet) roi_df <- roi_sf roi_df$geometry <- NULL - + if(length(aggregation_fun) > 1){ names(r_agg) <- paste0("ntl_", names(r_agg)) r_agg <- bind_cols(r_agg, roi_df) @@ -778,33 +778,33 @@ bm_extract <- function(roi_sf, roi_df[[paste0("ntl_", aggregation_fun)]] <- r_agg r_agg <- roi_df } - + if(add_n_pixels){ - + r_n_obs <- exact_extract(r, roi_sf, function(values, coverage_fraction) sum(!is.na(values)), progress = !quiet) - + r_n_obs_poss <- exact_extract(r, roi_sf, function(values, coverage_fraction) length(values), progress = !quiet) - + r_agg$n_pixels <- r_n_obs_poss r_agg$n_non_na_pixels <- r_n_obs - r_agg$prop_non_na_pixels <- r_agg$n_non_na_pixels / r_agg$n_pixels + r_agg$prop_non_na_pixels <- r_agg$n_non_na_pixels / r_agg$n_pixels } - + r_agg$date <- date_i - + #### Export saveRDS(r_agg, out_path) - + } else{ warning(paste0('"', out_path, '" already exists; skipping.\n')) } - + r_out <- NULL # Saving as file, so output from function should be NULL - + } else{ r_out <- bm_raster_i(roi_sf = roi_sf, product_id = product_id, @@ -817,80 +817,80 @@ bm_extract <- function(roi_sf, quiet = quiet, temp_dir = temp_dir) names(r_out) <- date_name_i - + if(add_n_pixels){ - + r_n_obs <- exact_extract(r_out, roi_sf, function(values, coverage_fraction) sum(!is.na(values)), progress = !quiet) - + r_n_obs_poss <- exact_extract(r_out, roi_sf, function(values, coverage_fraction) length(values), progress = !quiet) - + roi_sf$n_pixels <- r_n_obs_poss roi_sf$n_non_na_pixels <- r_n_obs - roi_sf$prop_non_na_pixels <- roi_sf$n_non_na_pixels / roi_sf$n_pixels + roi_sf$prop_non_na_pixels <- roi_sf$n_non_na_pixels / roi_sf$n_pixels } - + r_out <- exact_extract(x = r_out, y = roi_sf, fun = aggregation_fun, progress = !quiet) - + roi_df <- roi_sf roi_df$geometry <- NULL - + if(length(aggregation_fun) > 1){ names(r_out) <- paste0("ntl_", names(r_out)) r_out <- bind_cols(r_out, roi_df) } else{ - + roi_df[[paste0("ntl_", aggregation_fun)]] <- r_out r_out <- roi_df } - + r_out$date <- date_i } - + return(r_out) - + ## HERE # }, # error=function(e) { # return(NULL) # } #) - + }) - + # Clean output --------------------------------------------------------------- # Remove NULLs r_list <- r_list[!sapply(r_list,is.null)] - + r <- r_list %>% bind_rows() - + } - + # Output dataframe when output_location_type = "file" ------------------------ if(output_location_type == "file"){ if(!file_return_null){ - + ## Output path date_names <- define_date_name(date, product_id) - + out_name_end <- paste0("_", date_names, ".Rds") out_name <- paste0(out_name_begin, out_name_end) - + ## Restrict to files that exist, in case user enters date with no data ## (eg, year after current year) all_files <- list.files(file_dir) out_name <- out_name[out_name %in% all_files] - + r <- file.path(file_dir, out_name) %>% map_df(readRDS) - + # r <- file_dir %>% # list.files(full.names = T, # pattern = paste0("*.Rds")) %>% @@ -900,7 +900,7 @@ bm_extract <- function(roi_sf, r <- NULL } } - + unlink(temp_dir, recursive = T) return(r) } @@ -993,7 +993,7 @@ bm_extract <- function(roi_sf, #' @rawNamespace import(tidyr, except = c(extract)) #' @rawNamespace import(purrr, except = c(flatten_df, values)) #' @rawNamespace import(terra, except = c(intersect, values, origin, union)) -#' +#' # @rawNamespace import(utils, except = c(stack, unstack)) bm_raster <- function(roi_sf, product_id, @@ -1011,29 +1011,29 @@ bm_raster <- function(roi_sf, h5_dir = NULL, quiet = FALSE, ...){ - + # Errors & Warnings ---------------------------------------------------------- if( (interpol_na == T) & (length(date) == 1) ){ stop("If interpol_na = TRUE, then must have more than one date") } - + if( (interpol_na == T) & (output_location_type == "file") ){ interpol_na <- F stop("interpol_na ignored. Interpolation only occurs when output_location_type = 'memory'") } - + if(class(roi_sf)[1] == "SpatVector") roi_sf <- roi_sf %>% st_as_sf() if(!("sf" %in% class(roi_sf))){ stop("roi must be an sf object") } - - # Required parameters used in try statement, so error not generated when used, + + # Required parameters used in try statement, so error not generated when used, # so use them here roi_sf <- roi_sf product_id <- product_id date <- date bearer <- bearer - + # Assign interpolation variables --------------------------------------------- if(interpol_na == T){ if(!exists("method")) method <- "linear" @@ -1043,54 +1043,54 @@ bm_raster <- function(roi_sf, if(!exists("z")) z <- NULL if(!exists("NArule")) NArule <- 1 } - + # Define Tempdir ------------------------------------------------------------- temp_main_dir = tempdir() - + current_time_millis = as.character(as.numeric(Sys.time())*1000) %>% str_replace_all("[:punct:]", "") temp_dir = file.path(temp_main_dir, paste0("bm_raster_temp_", current_time_millis)) - + dir.create(temp_dir, showWarnings = F) - + # NTL Variable --------------------------------------------------------------- variable <- define_variable(variable, product_id) - + # Filename root -------------------------------------------------------------- # Define outside of lapply, as use this later to aggregate rasters if(output_location_type == "file"){ - out_name_begin <- paste0(file_prefix, - product_id, "_", + out_name_begin <- paste0(file_prefix, + product_id, "_", variable, "_", - "qflag", + "qflag", quality_flag_rm %>% paste0(collapse="_")) } - + # Download data -------------------------------------------------------------- r_list <- lapply(date, function(date_i){ - + out <- tryCatch( { - + #### Make name for raster based on date date_name_i <- define_date_name(date_i, product_id) - + #### If save as tif format if(output_location_type == "file"){ - + ## Output path out_name_end <- paste0("_", date_name_i, ".tif") out_name <- paste0(out_name_begin, out_name_end) - + out_path <- file.path(file_dir, out_name) - + make_raster <- TRUE if(file_skip_if_exists & file.exists(out_path)) make_raster <- FALSE - + if(make_raster){ - + r <- bm_raster_i(roi_sf = roi_sf, product_id = product_id, date = date_i, @@ -1102,17 +1102,17 @@ bm_raster <- function(roi_sf, quiet = quiet, temp_dir = temp_dir) names(r) <- date_name_i - + writeRaster(r, out_path) - + } else{ message(paste0('"', out_path, '" already exists; skipping.\n')) } - + r_out <- NULL # Saving as tif file, so output from function should be NULL - + } else{ - + r_out <- bm_raster_i(roi_sf = roi_sf, product_id = product_id, date = date_i, @@ -1124,23 +1124,23 @@ bm_raster <- function(roi_sf, quiet = quiet, temp_dir = temp_dir) names(r_out) <- date_name_i - + } - + return(r_out) - + }, error=function(e) { return(NULL) } ) - + }) - + # Clean output --------------------------------------------------------------- # Remove NULLs r_list <- r_list[!sapply(r_list,is.null)] - + if(length(r_list) == 1){ r <- r_list[[1]] } else if (length(r_list) > 1){ @@ -1148,7 +1148,7 @@ bm_raster <- function(roi_sf, } else{ r <- NULL } - + # Interpolate ---------------------------------------------------------------- if(interpol_na %in% T){ r <- terra::approximate(r, @@ -1159,21 +1159,21 @@ bm_raster <- function(roi_sf, z = z, NArule = NArule) } - + unlink(temp_dir, recursive = T) - + # Output raster when output_location_type = "file" --------------------------- if(output_location_type == "file"){ if(!file_return_null){ - + ## Output path date_names <- define_date_name(date, product_id) - + out_name_end <- paste0("_", date_names, ".tif") out_name <- paste0(out_name_begin, out_name_end) - + ## Restrict to files that exist, in case user enters date with no data ## (eg, year after current year) all_files <- list.files(file_dir) @@ -1186,7 +1186,7 @@ bm_raster <- function(roi_sf, r <- NULL } } - + return(r) } @@ -1200,29 +1200,29 @@ bm_raster_i <- function(roi_sf, h5_dir, quiet, temp_dir){ - + # Black marble grid ---------------------------------------------------------- bm_tiles_sf <- read_sf("https://raw.githubusercontent.com/worldbank/blackmarbler/main/data/blackmarbletiles.geojson") - + # Prep dates ----------------------------------------------------------------- ## For monthly, allow both yyyy-mm and yyyy-mm-dd (where -dd is ignored) if(product_id == "VNP46A3"){ - + if(nchar(date) %in% 7){ date <- paste0(date, "-01") } - + } - + ## For year, allow both yyyy and yyyy-mm-dd (where -mm-dd is ignored) if(product_id == "VNP46A4"){ - + if(nchar(date) %in% 4){ date <- paste0(date, "-01-01") } - + } - + # Grab tile dataframe -------------------------------------------------------- year <- date %>% year() month <- date %>% month() @@ -1233,19 +1233,19 @@ bm_raster_i <- function(roi_sf, years = year, months = month, days = day) - - + + # Intersecting tiles --------------------------------------------------------- # Remove grid along edges, which causes st_intersects to fail bm_tiles_sf <- bm_tiles_sf[!(bm_tiles_sf$TileID %>% str_detect("h00")),] bm_tiles_sf <- bm_tiles_sf[!(bm_tiles_sf$TileID %>% str_detect("v00")),] - - + + inter <- tryCatch( { inter <- st_intersects(bm_tiles_sf, roi_sf, sparse = F) %>% apply(1, sum) - + inter }, error = function(e){ @@ -1253,44 +1253,79 @@ bm_raster_i <- function(roi_sf, stop("Issue with `roi_sf` intersecting with blackmarble tiles; try buffering by a width of 0: eg, st_buffer(roi_sf, 0)") } ) - + grid_use_sf <- bm_tiles_sf[inter > 0,] - + # Make Raster ---------------------------------------------------------------- tile_ids_rx <- grid_use_sf$TileID %>% paste(collapse = "|") bm_files_df <- bm_files_df[bm_files_df$name %>% str_detect(tile_ids_rx),] - + if( (nrow(bm_files_df) < nrow(grid_use_sf)) & check_all_tiles_exist){ warning("Not all satellite imagery tiles for this location exist, so skipping. To ignore this error and process anyway, set check_all_tiles_exist = FALSE") stop("Not all satellite imagery tiles for this location exist, so skipping. To ignore this error and process anyway, set check_all_tiles_exist = FALSE") } - + unlink(file.path(temp_dir, product_id), recursive = T) - + if(quiet == F){ message(paste0("Processing ", nrow(bm_files_df), " nighttime light tiles")) } - - + + r_list <- lapply(bm_files_df$name, function(name_i){ download_raster(name_i, temp_dir, variable, bearer, quality_flag_rm, h5_dir, quiet) }) - - - + + + if(length(r_list) == 1){ r <- r_list[[1]] } else{ - + r <- do.call(terra::mosaic, c(r_list, fun = "max")) } - + ## Crop r <- r %>% terra::crop(roi_sf) - + unlink(file.path(temp_dir, product_id), recursive = T) - + return(r) } +#' @title Get a NASA Earthdata bearer token +#' @description Fetch a NASA Earthdata bearer token from using the Earthdata API. +#' If none exist, this will create one, or if one already exists it will fetch +#' that one instead. +#' +#' @param username character. NASA Earthdata username +#' @param username character. NASA Earthdata password +#' +#' @return character +#' @author Simon E. H. Smart +#' @export +get_nasa_token <- function(username, password) { + + stopifnot("username must be a character string" = inherits(username, "character")) + stopifnot("password must be a character string" = inherits(password, "character")) + + token_url <- "https://urs.earthdata.nasa.gov/api/users/find_or_create_token" + req <- httr2::request(token_url) + + response <- tryCatch( + req %>% + httr2::req_auth_basic(username, password) %>% + httr2::req_method("POST") %>% + httr2::req_perform(), + httr2_http_401 = function(cnd){stop("Incorrect username or password", call. = FALSE)} + ) + + if (httr2::resp_status(response) == 200){ + body <- response %>% httr2::resp_body_json() + token <- body$access_token + return(token) + } else { + stop("Something went wrong when trying to retrieve the token") + } +} diff --git a/tests/testthat/test-bm_raster.R b/tests/testthat/test-bm_raster.R index ea5a2b2..3e09e71 100644 --- a/tests/testthat/test-bm_raster.R +++ b/tests/testthat/test-bm_raster.R @@ -1,88 +1,104 @@ test_that("Query VNP46A1", { - + skip() - + # Define bearer token bearer <- Sys.getenv("BEARER_NASA_TOKEN") - + # sf polygon of Switzerland (covers 2 black marble tiles) roi_sf <- geodata::gadm(country = "CHE", level = 0, path = tempdir()) - + # Daily data r <- bm_raster(roi_sf = roi_sf, product_id = "VNP46A1", date = "2021-10-03", bearer = bearer) - + expect_true(class(r) == "SpatRaster", info = "r is not a SpatRaster object" ) - + }) test_that("Query VNP46A2", { - + skip() - + # Define bearer token bearer <- Sys.getenv("BEARER_NASA_TOKEN") - + # sf polygon of Switzerland (covers 2 black marble tiles) roi_sf <- geodata::gadm(country = "CHE", level = 0, path = tempdir()) - + # Daily data r <- bm_raster(roi_sf = roi_sf, product_id = "VNP46A2", date = "2021-10-03", bearer = bearer) - + expect_true(class(r) == "SpatRaster", info = "r is not a SpatRaster object" ) - + }) test_that("Query VNP46A3", { - + skip() - + # Define bearer token bearer <- Sys.getenv("BEARER_NASA_TOKEN") - + # sf polygon of Switzerland (covers 2 black marble tiles) roi_sf <- geodata::gadm(country = "CHE", level = 0, path = tempdir()) - + # Daily data r <- bm_raster(roi_sf = roi_sf, product_id = "VNP46A3", date = "2021-10", bearer = bearer) - + expect_true(class(r) == "SpatRaster", info = "r is not a SpatRaster object" ) - + }) test_that("Query VNP46A4", { - + skip() - + # Define bearer token bearer <- Sys.getenv("BEARER_NASA_TOKEN") - + # sf polygon of Switzerland (covers 2 black marble tiles) roi_sf <- geodata::gadm(country = "CHE", level = 0, path = tempdir()) - + # Daily data r <- bm_raster(roi_sf = roi_sf, product_id = "VNP46A4", date = 2021, bearer = bearer) - + expect_true(class(r) == "SpatRaster", info = "r is not a SpatRaster object" ) - + }) +test_that("Get NASA token", { + + expect_error(get_nasa_token(123, "black"), "username must be a character string") + expect_error(get_nasa_token("black", 123), "password must be a character string") + expect_error(get_nasa_token("black", "marble"), "Incorrect username or password") + + username <- Sys.getenv("NASA_USERNAME") + password <- Sys.getenv("NASA_PASSWORD") + + token <- get_nasa_token(username, password) + expect_type(token, "character") + expect_gt(nchar(token), 100) + +}) + +