packages_to_check <- c("stars", "httr", "jsonlite", "tmap")
+packages_to_check <- c("stars", "httr", "jsonlite", "tmap")
# Check and install packages
for (package_name in packages_to_check) {
if (!package_name %in% rownames(installed.packages())) {
install.packages(package_name)
- cat(paste("Package", package_name, "installed.\n"))
+ cat(paste("Package", package_name, "installed.\n"))
} else {
- cat(paste("Package", package_name, "is already installed.\n"))
+ cat(paste("Package", package_name, "is already installed.\n"))
}
library(package_name, character.only = TRUE)
}
-Package stars is already installed.
-
-
-Loading required package: abind
-
-
-Loading required package: sf
-
-
-Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
-
-
-Package httr is already installed.
+Package stars is already installed.
+Package httr is already installed.
Package jsonlite is already installed.
Package tmap is already installed.
-
-Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
-remotes::install_github('r-tmap/tmap')
-
-
+
#in case tmap does not install
-#remotes::install_github('r-tmap/tmap')#in case tmap does not install
+#remotes::install_github('r-tmap/tmap')Reading the Data
Check what days are available: LANCE NRT
Based on availability, edit the year_day variable YYYY-DD. Example: ‘2022-01’
-
+
#add the year and date you want to search for (YYYY-DD, 2022-01)
-year_day <- '2023-336'#add the year and date you want to search for (YYYY-DD, 2022-01)
+year_day <- '2023-336'Determine tiles of interest: MODIS NRT Tile Map
-
Based on availability, edit the tile_code variable:
-
+
#add tile code from the map above (written as h00v00)
-tile_code <- 'h05v05'#add tile code from the map above (written as h00v00)
+tile_code <- 'h05v05'This is the NRT Flood F3 (MCDWD_L3_F3) API URL:
-
+
API_URL <- paste0('https://nrt3.modaps.eosdis.nasa.gov/api/v2/content/details?products=MCDWD_L3_F3_NRT&archiveSets=61&temporalRanges=')API_URL <- paste0('https://nrt3.modaps.eosdis.nasa.gov/api/v2/content/details?products=MCDWD_L3_F3_NRT&archiveSets=61&temporalRanges=')We can combine the API URL above with the year_day provided and print the available datasets:
-
+
#pasting together URL and year_day
-url <- paste0(API_URL, year_day)
-print(url)#pasting together URL and year_day
+url <- paste0(API_URL, year_day)
+print(url)
-
[1] "https://nrt3.modaps.eosdis.nasa.gov/api/v2/content/details?products=MCDWD_L3_F3_NRT&archiveSets=61&temporalRanges=2023-336"[1] "https://nrt3.modaps.eosdis.nasa.gov/api/v2/content/details?products=MCDWD_L3_F3_NRT&archiveSets=61&temporalRanges=2023-336"Load Data
Access the NASA Earthdata with the GET function:
-
-
if(!file.exists("modis_nrt_flood.nc")) {
- # Make the GET request
- response <- httr::GET(url)
- # Check response status from the GET function and check the contents from the parsed data.
- print(response)
- if (http_status(response)$category == "Success") {
- # Parse the response JSON
- data <- content(response, as = "text", encoding = "UTF-8")
- data_parsed <- jsonlite::fromJSON(data)
- #filter for the tile code
- content_items <- data_parsed$content[grepl(tile_code, data_parsed$content$name, ignore.case = TRUE), ]
- #check the content items
- print(content_items)
- #Select the URL from the 'downloadsLink' column in the content_items list:
- download_link <- content_items$downloadsLink
- print(download_link)
-
-
- } else {
- print("Request failed with status code", http_status(response)$status_code)
- }
-
-} else{
- download_link <- "modis_nrt_flood.nc"
-}
-
+Response [https://nrt3.modaps.eosdis.nasa.gov/api/v2/content/details?products=MCDWD_L3_F3_NRT&archiveSets=61&temporalRanges=2023-336]
- Date: 2023-12-14 16:37
- Status: 200
- Content-Type: application/json;charset=UTF-8
- Size: 112 kB
-
- archiveSets cksum collections dataDay
-57 61 1275757885 modis-nrt-c6.1 2023-336 = 2023-12-02
- downloadsLink
-57 https://nrt3.modaps.eosdis.nasa.gov/api/v2/content/archives/MCDWD_L3_F3_NRT.A2023336.h05v05.061.tif
- fileId md5sum mtime
-57 2051979420 64a8ba49c9893388c1da7217df15ec77 1701612811
- name products resourceType
-57 MCDWD_L3_F3_NRT.A2023336.h05v05.061.tif MCDWD_L3_F3_NRT File
- self size
-57 /api/v2/content/details/MCDWD_L3_F3_NRT.A2023336.h05v05.061.tif 1080978
-[1] "https://nrt3.modaps.eosdis.nasa.gov/api/v2/content/archives/MCDWD_L3_F3_NRT.A2023336.h05v05.061.tif"if(!file.exists("modis_nrt_flood.nc")) {
+ # Make the GET request
+ response <- httr::GET(url)
+ # Check response status from the GET function and check the contents from the parsed data.
+ print(response)
+ if (http_status(response)$category == "Success") {
+ # Parse the response JSON
+ data <- content(response, as = "text", encoding = "UTF-8")
+ data_parsed <- jsonlite::fromJSON(data)
+ #filter for the tile code
+ content_items <- data_parsed$content[grepl(tile_code, data_parsed$content$name, ignore.case = TRUE), ]
+ #check the content items
+ print(content_items)
+ #Select the URL from the 'downloadsLink' column in the content_items list:
+ download_link <- content_items$downloadsLink
+ print(download_link)
+
+
+ } else {
+ print("Request failed with status code", http_status(response)$status_code)
+ }
+
+} else{
+ download_link <- "modis_nrt_flood.nc"
+}Use the “read_stars()” function from the “stars” R Library to read the geoTiff raster. The raster is assigned to the “x” variable:
-
-
raster_df <- stars::read_stars(download_link)
-
+Warning: ignoring unrecognized unit: noneraster_df <- stars::read_stars(download_link)Set the Coordinate reference system (CRS) to “EPSG:4326”
-
+
-
+
@@ -301,12 +4683,12 @@ my_raster <- st_set_crs(raster_df, st_crs("EPSG:4326"))my_raster <- st_set_crs(raster_df, st_crs("EPSG:4326"))Warning in `st_crs<-.dimensions`(`*tmp*`, value = value): replacing CRS does
not reproject data: use st_transform, or st_warp to warp to a new CRSst_crs(my_raster)st_crs(my_raster)Coordinate Reference System:
User input: EPSG:4326
wkt:
-GEOGCRS["WGS 84",
- ENSEMBLE["World Geodetic System 1984 ensemble",
- MEMBER["World Geodetic System 1984 (Transit)"],
- MEMBER["World Geodetic System 1984 (G730)"],
- MEMBER["World Geodetic System 1984 (G873)"],
- MEMBER["World Geodetic System 1984 (G1150)"],
- MEMBER["World Geodetic System 1984 (G1674)"],
- MEMBER["World Geodetic System 1984 (G1762)"],
- MEMBER["World Geodetic System 1984 (G2139)"],
- ELLIPSOID["WGS 84",6378137,298.257223563,
- LENGTHUNIT["metre",1]],
+GEOGCRS["WGS 84",
+ ENSEMBLE["World Geodetic System 1984 ensemble",
+ MEMBER["World Geodetic System 1984 (Transit)"],
+ MEMBER["World Geodetic System 1984 (G730)"],
+ MEMBER["World Geodetic System 1984 (G873)"],
+ MEMBER["World Geodetic System 1984 (G1150)"],
+ MEMBER["World Geodetic System 1984 (G1674)"],
+ MEMBER["World Geodetic System 1984 (G1762)"],
+ MEMBER["World Geodetic System 1984 (G2139)"],
+ ELLIPSOID["WGS 84",6378137,298.257223563,
+ LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
- PRIMEM["Greenwich",0,
- ANGLEUNIT["degree",0.0174532925199433]],
+ PRIMEM["Greenwich",0,
+ ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
- AXIS["geodetic latitude (Lat)",north,
+ AXIS["geodetic latitude (Lat)",north,
ORDER[1],
- ANGLEUNIT["degree",0.0174532925199433]],
- AXIS["geodetic longitude (Lon)",east,
+ ANGLEUNIT["degree",0.0174532925199433]],
+ AXIS["geodetic longitude (Lon)",east,
ORDER[2],
- ANGLEUNIT["degree",0.0174532925199433]],
+ ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
- SCOPE["Horizontal component of 3D system."],
- AREA["World."],
+ SCOPE["Horizontal component of 3D system."],
+ AREA["World."],
BBOX[-90,-180,90,180]],
- ID["EPSG",4326]]Load Data
Visualizing NRT Flood Data
Plot the raster to quickly view it:
-
+
plot(my_raster, axes = TRUE)plot(my_raster, axes = TRUE)
-
downsample set to 10downsample set to 3
-
Create NRT
Create a classified legend; however, the NRT Flood data is stored in decimal numbers (aka floating-point). Create class breaks dividing the data by these breaks, and corresponding colors and labels:
-
+
class_breaks <- c(-Inf, 0.1, 1.1, 2.1, 3.1)
-colors <- c( "gray", "blue", "yellow","red")
-
-labels = c("0: No Water", "1: Surface Water", "2: Recurring flood", "3: Flood (unusual)")class_breaks <- c(-Inf, 0.1, 1.1, 2.1, 3.1)
+colors <- c( "gray", "blue", "yellow","red")
+
+labels = c("0: No Water", "1: Surface Water", "2: Recurring flood", "3: Flood (unusual)")Add a title for the plot that includes the year, day of year, and tile code:
-
+
title = paste("NRT Flood", year_day, tile_code)title = paste("NRT Flood", year_day, tile_code)Generate a plot from the tmap library using the tm_shape() function. With style as “cat,” meaning categorical. T
-
+
-
+
+
+
tmap_mode("view")tmap_mode("view")tmap mode set to interactive viewing## tmap mode set to plotting
-tm_shape(my_raster, style="cat" )+
- tm_raster(palette = c(colors),
- title = title,
- breaks = class_breaks,
- labels = labels )+
- tm_basemap(server = "Esri.WorldImagery") +
- tm_layout(legend.outside = TRUE)## tmap mode set to plotting
+tm_shape(my_raster, style="cat" )+
+ tm_raster(palette = c(colors),
+ title = title,
+ breaks = class_breaks,
+ labels = labels )+
+ tm_basemap(server = "Esri.WorldImagery") +
+ tm_layout(legend.outside = TRUE)stars object downsampled to 1000 by 1000 cells. See tm_shape manual (argument raster.downsample)
-
-
+
+
Create NRT
Value 2 (Recurring flood) is not populated in the beta release.↩︎
-
Value 2 (Recurring flood) is not populated in the beta release.↩︎