diff --git a/DESCRIPTION b/DESCRIPTION index 7a94566..60cd442 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: blackmarbler Title: Black Marble Data and Statistics -Version: 0.2.2 +Version: 0.2.4 Authors@R: c(person("Robert", "Marty", , "rmarty@worldbank.org", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-3164-3813")), @@ -15,7 +15,6 @@ URL: https://worldbank.github.io/blackmarbler/ BugReports: https://github.com/worldbank/blackmarbler/issues Imports: readr, - hdf5r, dplyr, purrr, lubridate, @@ -24,7 +23,6 @@ Imports: sf, exactextractr, stringr, - httr, httr2 Suggests: geodata, diff --git a/NAMESPACE b/NAMESPACE index 1ebb2c7..51c68c6 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -5,8 +5,6 @@ export(bm_raster) export(get_nasa_token) import(dplyr) import(exactextractr) -import(hdf5r) -import(httr) import(httr2) import(lubridate) import(purrr, except = c(flatten_df, values)) diff --git a/R/blackmarbler.R b/R/blackmarbler.R index abc2218..febaba4 100644 --- a/R/blackmarbler.R +++ b/R/blackmarbler.R @@ -182,6 +182,141 @@ apply_scaling_factor <- function(x, variable){ return(x) } +# file_to_raster <- function(h5_file, +# variable, +# quality_flag_rm){ +# # 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)]][,] +# +# # VNP46A1 does not have Mandatory_Quality_Flag +# if(h5_file %>% str_detect("VNP46A2")){ +# 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 +# } +# } +# } +# } +# +# #### 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) +# } + file_to_raster <- function(h5_file, variable, quality_flag_rm){ @@ -189,70 +324,42 @@ file_to_raster <- function(h5_file, # ARGS # --f: Filepath to h5 file - ## Data - h5_data <- h5file(h5_file, "r+") + # Load data - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + h5_data <- terra::rast(h5_file) - #### 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)]][,] - - # VNP46A1 does not have Mandatory_Quality_Flag + # Get extent - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + 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() + + # Load raster of select variable - - - - - - - - - - - - - - - - - - - - - - - + var_names <- names(h5_data) + + 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[[variable]] + + # Filter by quality - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + if(length(quality_flag_rm) > 0){ + #### Get quality raster if(h5_file %>% str_detect("VNP46A2")){ - qf <- h5_data[["HDFEOS/GRIDS/VNP_Grid_DNB/Data Fields/Mandatory_Quality_Flag"]][,] + qf <- h5_data$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 - } - } - } - } - - #### 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){ + } else{ variable_short <- variable %>% str_replace_all("_Num", "") %>% @@ -260,81 +367,55 @@ file_to_raster <- function(h5_file, 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 - } - - } - + qf <- h5_data[[qf_name]] } - if(class(out[1,1])[1] != "numeric"){ - out <- matrix(as.numeric(out), # Convert to numeric matrix - ncol = ncol(out)) + #### Filter + for(val in quality_flag_rm){ # out[qf %in% quality_flag_rm] doesn't work, so loop + out[qf == val] <- NA } - - 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" + # Add CRS and Extent - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + terra::crs(out) <- "EPSG:4326" + terra::ext(out) <- c(xMin,xMax,yMin,yMax) - ## 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)) + # Additional filtering - - - - - - - - - - - - - - - - - - - - - - - - - - - - #set fill values to NA - outr <- remove_fill_value(outr, variable) + out <- remove_fill_value(out, variable) #apply scaling factor - outr <- apply_scaling_factor(outr, variable) + out <- apply_scaling_factor(out, variable) - #h5closeAll() - h5_data$close_all() - - return(outr) + return(out) } read_bm_csv <- function(year, day, product_id){ - - - 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 + # + 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){ + #warning(paste0("Error with year: ", year, "; day: ", day)) + data.frame(NULL) + } + ) Sys.sleep(0.1) - return(df) + return(df_out) } create_dataset_name_df <- function(product_id, @@ -435,7 +516,7 @@ download_raster <- function(file_name, url <- paste0('https://ladsweb.modaps.eosdis.nasa.gov/archive/allData/5000/', product_id, '/', year, '/', day, '/', file_name) - headers <- c('Authorization' = paste('Bearer', bearer)) + #headers <- c('Authorization' = paste('Bearer', bearer)) if(is.null(h5_dir)){ download_path <- file.path(temp_dir, file_name) @@ -447,30 +528,23 @@ download_raster <- function(file_name, if(quiet == FALSE) message(paste0("Processing: ", file_name)) - if(quiet == TRUE){ - - response <- httr::GET(url, - httr::timeout(60), - httr::add_headers(headers), - httr::write_disk(download_path, overwrite = TRUE)) - - } else{ - response <- httr::GET(url, - httr::timeout(60), - httr::add_headers(headers), - httr::write_disk(download_path, overwrite = TRUE), - httr::progress()) - - } + response <- httr2::request(url) %>% + httr2::req_headers('Authorization' = paste('Bearer', bearer)) %>% + httr2::req_timeout(60) %>% + httr2::req_perform() if(response$status_code != 200){ - message("Error in downloading data") message(response) + stop("Error in downloading data") + } + + if(response$status_code == 200){ + if(length(response$body) < 10000){ + stop(paste0("Issue with bearer token. You may need to generate a new token. Ensure that select EULAs are accepted. Please see the instructions here: https://github.com/worldbank/blackmarbler?tab=readme-ov-file#bearer-token-")) + } } - #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-") - #} + writeBin(httr2::resp_body_raw(response), download_path) } @@ -578,7 +652,7 @@ get_nasa_token <- function(username, password) { #' * For `product_id` `:VNP46A1"`, uses `DNB_At_Sensor_Radiance_500m`. #' * For `product_id` `"VNP46A2"`, uses `Gap_Filled_DNB_BRDF-Corrected_NTL`. #' * For `product_id`s `"VNP46A3"` and `"VNP46A4"`, uses `NearNadir_Composite_Snow_Free`. -#' For information on other variable choices, see [here](https://ladsweb.modaps.eosdis.nasa.gov/api/v2/content/archives/Document%20Archive/Science%20Data%20Product%20Documentation/VIIRS_Black_Marble_UG_v1.2_April_2021.pdf); for `VNP46A1`, see Table 3; for `VNP46A2` see Table 6; for `VNP46A3` and `VNP46A4`, see Table 9. +#' To see all variable choices, set `variable = ""` (this will create an error message that lists all valid variables). For additional information on variable choices, see [here](https://ladsweb.modaps.eosdis.nasa.gov/api/v2/content/archives/Document%20Archive/Science%20Data%20Product%20Documentation/VIIRS_Black_Marble_UG_v1.2_April_2021.pdf); for `VNP46A1`, see Table 3; for `VNP46A2` see Table 6; for `VNP46A3` and `VNP46A4`, see Table 9. #' @param quality_flag_rm Quality flag values to use to set values to `NA`. Each pixel has a quality flag value, where low quality values can be removed. Values are set to `NA` for each value in the `quality_flag_rm` vector. Note that `quality_flag_rm` does not apply for `VNP46A1`. (Default: `NULL`). #' #' @@ -605,7 +679,7 @@ get_nasa_token <- function(username, password) { #' @param ... Additional arguments for `terra::approximate`, if `interpol_na = TRUE` #' #' @return Raster -#' +#' @author Robert Marty #' @examples #' \dontrun{ #' # Define bearer token @@ -760,131 +834,133 @@ bm_extract <- function(roi_sf, # Download data -------------------------------------------------------------- r_list <- lapply(date, function(date_i){ - out <- tryCatch( - { + # 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 name for raster based on date - date_name_i <- define_date_name(date_i, product_id) + #### Make raster + r <- bm_raster_i(roi_sf = roi_sf, + product_id = product_id, + date = date_i, + bearer = bearer, + variable = variable, + quality_flag_rm = quality_flag_rm, + check_all_tiles_exist = check_all_tiles_exist, + h5_dir = h5_dir, + quiet = quiet, + temp_dir = temp_dir) + names(r) <- date_name_i - #### 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, - date = date_i, - bearer = bearer, - variable = variable, - quality_flag_rm = quality_flag_rm, - check_all_tiles_exist = check_all_tiles_exist, - h5_dir = h5_dir, - quiet = quiet, - temp_dir = temp_dir) - names(r) <- date_name_i - - #### Extract - 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) - } else{ - 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$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 - + #### Extract + 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) } else{ - r_out <- bm_raster_i(roi_sf = roi_sf, - product_id = product_id, - date = date_i, - bearer = bearer, - variable = variable, - quality_flag_rm = quality_flag_rm, - check_all_tiles_exist = check_all_tiles_exist, - h5_dir = h5_dir, - 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 - } - - r_out <- exact_extract(x = r_out, y = roi_sf, fun = aggregation_fun, - progress = !quiet) + roi_df[[paste0("ntl_", aggregation_fun)]] <- r_agg + r_agg <- roi_df + } + + if(add_n_pixels){ - roi_df <- roi_sf - roi_df$geometry <- NULL + r_n_obs <- exact_extract(r, roi_sf, function(values, coverage_fraction) + sum(!is.na(values)), + progress = !quiet) - 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_n_obs_poss <- exact_extract(r, roi_sf, function(values, coverage_fraction) + length(values), + progress = !quiet) - r_out$date <- date_i + 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 } - return(r_out) + r_agg$date <- date_i + + #### Export + saveRDS(r_agg, out_path) - # HERE - }, - error=function(e) { - return(NULL) + } 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, + date = date_i, + bearer = bearer, + variable = variable, + quality_flag_rm = quality_flag_rm, + check_all_tiles_exist = check_all_tiles_exist, + h5_dir = h5_dir, + 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 + } + + 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) + + # TRY START + # + # }, + # error=function(e) { + # return(NULL) + # } + # ) + # TRY END }) @@ -945,7 +1021,7 @@ bm_extract <- function(roi_sf, #' * For `product_id` `:VNP46A1"`, uses `DNB_At_Sensor_Radiance_500m`. #' * For `product_id` `"VNP46A2"`, uses `Gap_Filled_DNB_BRDF-Corrected_NTL`. #' * For `product_id`s `"VNP46A3"` and `"VNP46A4"`, uses `NearNadir_Composite_Snow_Free`. -#' For information on other variable choices, see [here](https://ladsweb.modaps.eosdis.nasa.gov/api/v2/content/archives/Document%20Archive/Science%20Data%20Product%20Documentation/VIIRS_Black_Marble_UG_v1.2_April_2021.pdf); for `VNP46A1`, see Table 3; for `VNP46A2` see Table 6; for `VNP46A3` and `VNP46A4`, see Table 9. +#' To see all variable choices, set `variable = ""` (this will create an error message that lists all valid variables). For additional information on variable choices, see [here](https://ladsweb.modaps.eosdis.nasa.gov/api/v2/content/archives/Document%20Archive/Science%20Data%20Product%20Documentation/VIIRS_Black_Marble_UG_v1.2_April_2021.pdf); for `VNP46A1`, see Table 3; for `VNP46A2` see Table 6; for `VNP46A3` and `VNP46A4`, see Table 9. #' @param quality_flag_rm Quality flag values to use to set values to `NA`. Each pixel has a quality flag value, where low quality values can be removed. Values are set to `NA` for each value in the `quality_flag_rm` vector. Note that `quality_flag_rm` does not apply for `VNP46A1`. (Default: `NULL`). #' #' @@ -973,6 +1049,7 @@ bm_extract <- function(roi_sf, #' #' @return Raster #' +#' @author Robert Marty #' @examples #' \dontrun{ #' # Define bearer token @@ -1004,12 +1081,10 @@ bm_extract <- function(roi_sf, #' @export #' #' @import readr -#' @import hdf5r #' @import dplyr #' @import sf #' @import exactextractr #' @import stringr -#' @import httr #' @import httr2 #' @import lubridate #' @rawNamespace import(tidyr, except = c(extract)) @@ -1091,75 +1166,76 @@ bm_raster <- function(roi_sf, # 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) + #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){ - #### 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, - bearer = bearer, - variable = variable, - quality_flag_rm = quality_flag_rm, - check_all_tiles_exist = check_all_tiles_exist, - h5_dir = h5_dir, - 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, - bearer = bearer, - variable = variable, - quality_flag_rm = quality_flag_rm, - check_all_tiles_exist = check_all_tiles_exist, - h5_dir = h5_dir, - quiet = quiet, - temp_dir = temp_dir) - names(r_out) <- date_name_i - - } + r <- bm_raster_i(roi_sf = roi_sf, + product_id = product_id, + date = date_i, + bearer = bearer, + variable = variable, + quality_flag_rm = quality_flag_rm, + check_all_tiles_exist = check_all_tiles_exist, + h5_dir = h5_dir, + quiet = quiet, + temp_dir = temp_dir) + names(r) <- date_name_i - return(r_out) + terra::writeRaster(r, out_path) - }, - error=function(e) { - return(NULL) + } 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, + bearer = bearer, + variable = variable, + quality_flag_rm = quality_flag_rm, + check_all_tiles_exist = check_all_tiles_exist, + h5_dir = h5_dir, + quiet = quiet, + temp_dir = temp_dir) + names(r_out) <- date_name_i + + } + + return(r_out) + + # TRY START + # }, + # error=function(e) { + # return(NULL) + # } + # ) + # TRY END - #) }) @@ -1297,13 +1373,10 @@ bm_raster_i <- function(roi_sf, 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{ @@ -1319,3 +1392,7 @@ bm_raster_i <- function(roi_sf, return(r) } + + + + diff --git a/README.md b/README.md index d18b5ce..9b86128 100644 --- a/README.md +++ b/README.md @@ -45,43 +45,46 @@ devtools::install_github("worldbank/blackmarbler") ## Bearer Token -To obtain a bearer token, you'll need to have a registered [NASA Earth Data account](https://ladsweb.modaps.eosdis.nasa.gov/). On the webpage, click "login" and create a username/password if needed. +Follow the below steps to obtain a bearer token: -After an account is created, the NASA Bearer Token can be retrieved either using the `get_nasa_token` function or manually (see below). +1. Create a [NASA Earth Data account](https://ladsweb.modaps.eosdis.nasa.gov/) account. On the top right of the [webpage](https://ladsweb.modaps.eosdis.nasa.gov/), click "Login" then "Earthdata Login". Then click "register" (blue button). +2. Enter the information in the registration page. You __must__ include the following information; this information is not required to create an account, but the bearer token will not work without this information: -### Programmatically retrieve token - -The NASA Bearer Token can also be programmatically retrieved using the `get_nasa_token()` function. After making an account, the `get_nasa_token()` function uses your username and password to retrieve the Bearer token. - -```r -bearer <- get_nasa_token(username = "USERNAME-HERE", - password = "PASSWORD-HERE") -``` - -### Manually retrieve token - -The function requires using a **Earthdata Download Bearer Token**; to obtain a token, follow the below steps: + - Study Area + - User Type + - Organization + +3. Click "Register for EarthData Login" (green button at bottom). Check your email, and click the link in the email to activate the account. +4. Go to the [Earth Data Login](https://urs.earthdata.nasa.gov/users) page and login. +5. On the panel near the top, click "EULAs" then "Accept New EULAs". Accept: -1. Go to the [NASA LAADS Archive](https://ladsweb.modaps.eosdis.nasa.gov/missions-and-measurements/products/VNP46A4/) -2. Click "Login" (button on top right), and click "Earthdata Login"; create an account if needed. -3. Enter your username and password. -4. Click "Login", then "Generate Token" on the dropdown. -5. Click the "Generate Token" tab, then click the green button to generate a token -6. Click the blue "Show token" button; this is your bearer token. It will be a long string of text (over 500 characters). + - MERIS EULA + - Sentinel EULA + +6. On the "Profile Home" page, you should see something like below. Information should be filled in for each category, and "Agreed To Meris EULA" and "Agreed To Sentinel-3 EULA" should be True.

-NASA LAADS Bearer Token +NASA Profile Home Information

-

-NASA LAADS Bearer Token -

+7. Go to the [NASA LAADS Archive](https://ladsweb.modaps.eosdis.nasa.gov/missions-and-measurements/products/VNP46A4/) and login (login botton on top right). You will see a page to authorize use of Sentinel3 and Meris. Click the green "Authorize" button. + +8. To obtain the bearer token, go to the [Earth Data Login](https://urs.earthdata.nasa.gov/users) page and login. On the top panel, click "Generate token". On this page, click "Show Token" to see the bearer token.

NASA LAADS Bearer Token

+9. If the bearer token ever stops working, you make need to go to the "Generate token" page (see step 8), delete any existing tokens, and generate a new token. + +### Programmatically retrieve token + +After following the above steps, the bearer token can also be programmatically retrieved using the `get_nasa_token()` function and your usename and password. +```r +bearer <- get_nasa_token(username = "USERNAME-HERE", + password = "PASSWORD-HERE") +``` ## Usage @@ -289,7 +292,7 @@ Both functions take the following arguments: ### Optional arguments -* **variable:** Variable to used to create raster (default: `NULL`). For information on all variable choices, see [here](https://ladsweb.modaps.eosdis.nasa.gov/api/v2/content/archives/Document%20Archive/Science%20Data%20Product%20Documentation/VIIRS_Black_Marble_UG_v1.2_April_2021.pdf); for `VNP46A1`, see Table 3; for `VNP46A2` see Table 6; for `VNP46A3` and `VNP46A4`, see Table 9. If `NULL`, uses the following default variables: +* **variable:** Variable to used to create raster (default: `NULL`). To see all variable choices, set `variable = ""` (this will create an error message that lists all valid variables). For additional information on all variable choices, see [here](https://ladsweb.modaps.eosdis.nasa.gov/api/v2/content/archives/Document%20Archive/Science%20Data%20Product%20Documentation/VIIRS_Black_Marble_UG_v1.2_April_2021.pdf); for `VNP46A1`, see Table 3; for `VNP46A2` see Table 6; for `VNP46A3` and `VNP46A4`, see Table 9. If `NULL`, uses the following default variables: * For `product_id` `"VNP46A1"`, uses `DNB_At_Sensor_Radiance_500m`. * For `product_id` `"VNP46A2"`, uses `Gap_Filled_DNB_BRDF-Corrected_NTL`. diff --git a/man/bm_extract.Rd b/man/bm_extract.Rd index 08ff91d..720dbfb 100644 --- a/man/bm_extract.Rd +++ b/man/bm_extract.Rd @@ -128,3 +128,6 @@ ken_2021_r <- bm_extract(roi_sf = roi_sf, } } +\author{ +Robert Marty \href{mailto:rmarty@worldbank.org}{rmarty@worldbank.org} +} diff --git a/man/bm_raster.Rd b/man/bm_raster.Rd index 36406e1..2a4dc50 100644 --- a/man/bm_raster.Rd +++ b/man/bm_raster.Rd @@ -123,3 +123,6 @@ ken_2021_r <- bm_raster(roi_sf = roi_sf, } } +\author{ +Robert Marty \href{mailto:rmarty@worldbank.org}{rmarty@worldbank.org} +} diff --git a/man/figures/nasa_profile_info.png b/man/figures/nasa_profile_info.png new file mode 100644 index 0000000..db85554 Binary files /dev/null and b/man/figures/nasa_profile_info.png differ diff --git a/readme_figures/test_20250103.R b/readme_figures/test_20250103.R new file mode 100644 index 0000000..23d8865 --- /dev/null +++ b/readme_figures/test_20250103.R @@ -0,0 +1,47 @@ +# Test Package + +# INSTALL DEV VERSION FROM GITHUB +# devtools::install_github("worldbank/blackmarbler", auth_token = "") + +library(blackmarbler) +library(sf) +library(terra) +library(httr2) + +# ADD BEARER TOKEN HERE +bearer <- read.csv("~/Dropbox/bearer_bm.csv")$token + +# CHECK 1 ---------------------------------------------------------------------- +# With sessionInfo, under "other attached packages", make sure the blackmarbler +# version is: blackmarbler_0.2.4 +sessionInfo() + +# CHECK 2 ---------------------------------------------------------------------- +# What is the response? In particular, what is the: +# - Status +# - Content-Type +# - Body (how many bytes?) +url1 <- "https://ladsweb.modaps.eosdis.nasa.gov/archive/allData/5000/VNP46A3/2021/032/VNP46A3.A2021032.h21v09.001.2021144024822.h5" +response <- request(url1) %>% + req_headers('Authorization' = paste('Bearer', bearer)) %>% + req_timeout(60) %>% + req_perform() + +response + +length(response$body) < 10000 + +# CHECK 3 ---------------------------------------------------------------------- +roi_sf <- data.frame(lat = -1.943889, lon = 30.059444, id = 1) |> + st_as_sf(coords = c("lon", "lat"), + crs = 4326) |> + st_buffer(dist = 20000) + +r <- bm_raster(roi_sf = roi_sf, + product_id = "VNP46A2", + date = "2018-09-01", + bearer = bearer) + +plot(r) + + diff --git a/readme_figures/testing.R b/readme_figures/testing.R index 263fb88..9dbcd08 100644 --- a/readme_figures/testing.R +++ b/readme_figures/testing.R @@ -1,8 +1,9 @@ # Testing -#### Basic library(blackmarbler) library(sf) + +# Using package ---------------------------------------------------------------- bearer <- read.csv("~/Dropbox/bearer_bm.csv")$token roi_sf <- data.frame(lat = -1.943889, lon = 30.059444, id = 1) |> @@ -11,10 +12,32 @@ roi_sf <- data.frame(lat = -1.943889, lon = 30.059444, id = 1) |> st_buffer(dist = 20000) r_20210205 <- bm_raster(roi_sf = roi_sf, - product_id = "VNP46A4", - date = 2023:2024, + product_id = "VNP46A3", + date = "2018-09-01", bearer = bearer) + + + + +#### Basic + + +library(readr) +library(hdf5r) +library(dplyr) +library(purrr) +library(lubridate) +library(tidyr) +library(raster) +library(sf) +library(exactextractr) +library(stringr) +library(httr2) +library(httr) + + + library(readr) library(hdf5r) library(dplyr) diff --git a/readme_figures/testing_from_source.R b/readme_figures/testing_from_source.R new file mode 100644 index 0000000..ee746b0 --- /dev/null +++ b/readme_figures/testing_from_source.R @@ -0,0 +1,154 @@ +# Testing from Source + +#### Packages +library(readr) +library(hdf5r) +library(dplyr) +library(purrr) +library(lubridate) +library(tidyr) +#library(raster) +library(sf) +library(exactextractr) +library(stringr) +library(httr2) +#library(httr) +library(terra) + +source("~/Documents/Github/blackmarbler/R/blackmarbler.R") + +bearer <- read.csv("~/Dropbox/bearer_bm.csv")$token + +roi_sf <- data.frame(lat = -1.943889, lon = 30.059444, id = 1) |> + st_as_sf(coords = c("lon", "lat"), + crs = 4326) |> + st_buffer(dist = 20000) + +r_20210205 <- bm_raster(roi_sf = roi_sf, + product_id = "VNP46A3", + date = "2021-02-05", + bearer = bearer, + quiet = T) + +library(leaflet) +r <- r_20210205 +# Convert raster to a color palette +palette <- colorNumeric(palette = "viridis", domain = terra::values(r), na.color = "transparent") + +# Create Leaflet map and add raster +leaflet() %>% + addTiles() %>% # Add base map + addRasterImage(r, colors = palette, opacity = 0.8) %>% # Overlay raster + addLegend(pal = palette, values = terra::values(r), title = "Raster Values") + + +url1 <- "https://ladsweb.modaps.eosdis.nasa.gov/archive/allData/5000/VNP46A3/2021/032/VNP46A3.A2021032.h21v09.001.2021144024822.h5" +response <- request(url1) %>% + req_headers('Authorization' = paste('Bearer', bearer)) %>% + req_timeout(60) %>% + req_perform() + +# Assuming `response` is your HTTP response object +raw_data <- resp_body_raw(response) # Get the raw data from the response + +# Create a raw connection +raw_conn <- rawConnection(raw_data, open = "rb") + +result = base::unserialize(memDecompress(raw_conn)) + + +# Use `rast()` to load the data from the raw connection +h5_data <- rast(raw_data) +h5_data <- rast(raw_conn) +h5_data <- read_stars(raw_data) +h5_data <- read_stars(raw_conn) + +library(stars) + + + +close(raw_conn) + + +h5_data <- rast(resp_body_raw(response)) + +writeBin(resp_body_raw(response), "~/Desktop/VNP46A3.A2021036.h21v09.001.2021104201220.h5") + +h5_file <- "~/Desktop/VNP46A3.A2021036.h21v09.001.2021104201220.h5" +variable <- "AllAngle_Composite_Snow_Free" +quality_flag_rm <- NULL +library(terra) + +file_to_raster <- function(h5_file, + variable, + quality_flag_rm){ + # Converts h5 file to raster. + # ARGS + # --f: Filepath to h5 file + + # Load data - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + h5_data <- terra::rast(h5_file) + + # Get extent - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + 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() + + # Load raster of select variable - - - - - - - - - - - - - - - - - - - - - - - + var_names <- names(h5_data) + + 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[[variable]] + + # Filter by quality - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + if(length(quality_flag_rm) > 0){ + #### Get quality raster + if(h5_file %>% str_detect("VNP46A2")){ + qf <- h5_data$Mandatory_Quality_Flag + + } else{ + + variable_short <- variable %>% + str_replace_all("_Num", "") %>% + str_replace_all("_Std", "") + + qf_name <- paste0(variable_short, "_Quality") + + qf <- h5_data[[qf_name]] + } + + #### Filter + for(val in quality_flag_rm){ # out[qf %in% quality_flag_rm] doesn't work, so loop + out[qf == val] <- NA + } + } + + # Add CRS and Extent - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + crs(out) <- "EPSG:4326" + ext(out) <- c(xMin,xMax,yMin,yMax) + + # Additional filtering - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + #set fill values to NA + out <- remove_fill_value(out, variable) + + #apply scaling factor + out <- apply_scaling_factor(out, variable) + + return(out) +} diff --git a/readme_figures/testing_using_package.R b/readme_figures/testing_using_package.R new file mode 100644 index 0000000..ed711c5 --- /dev/null +++ b/readme_figures/testing_using_package.R @@ -0,0 +1,22 @@ +# Testing + +# devtools::install_github("worldbank/blackmarbler", auth_token = "") + +library(blackmarbler) +library(sf) +library(terra) + +# Using package ---------------------------------------------------------------- +bearer <- read.csv("~/Dropbox/bearer_bm.csv")$token + +roi_sf <- data.frame(lat = -1.943889, lon = 30.059444, id = 1) |> + st_as_sf(coords = c("lon", "lat"), + crs = 4326) |> + st_buffer(dist = 20000) + +r_20210205 <- bm_raster(roi_sf = roi_sf, + product_id = "VNP46A2", + date = "2018-09-01", + bearer = bearer) + +