LANCE_MODIS_NRT_GlobalFlood_MCDWD.qmd

---
title: "MODIS NRT Global Flood Product"
author: "Juan F. Martinez"
date: "December 14, 2023"
---

## MODIS NRT Global Flood Product

### Introduction

The MODIS/Aqua+Terra Global Flood Product L3 Near Real Time (NRT) 250m Global Flood Product (MCDWD_L3_NRT) (beta) provides daily maps of flooding globally. The product is provided over 3 compositing periods (1-day, 2-day, and 3-day) to minimize the impact of clouds and more rigorously identify flood water (the best composite will depend on the cloudiness for a particular event).

<img src="https://www.earthdata.nasa.gov/sites/default/files/imported/Flood_mekong.png" alt="MODIS NRT FLOOD" width="600"/>

[NASA EARTHDATA](https://www.earthdata.nasa.gov/learn/find-data/near-real-time/modis-nrt-global-flood-product)

[CRM SEARCH](https://cmr.earthdata.nasa.gov/search/concepts/C2018599131-LANCEMODIS.html)

```{r message=FALSE, warning=FALSE}
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"))
  } else {
    cat(paste("Package", package_name, "is already installed.\n"))
  }
  library(package_name, character.only = TRUE)
}


```

```{r}
#in case tmap does not install
#remotes::install_github('r-tmap/tmap')

```

### Reading the Data

#### Check what days are available: [LANCE NRT](https://nrt3.modaps.eosdis.nasa.gov/archive/allData/61/MCDWD_L3_F3_NRT/)

Based on availability, edit the year_day variable YYYY-DD. Example: '2022-01'

```{r}

#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](https://www.earthdata.nasa.gov/s3fs-public/2023-01/MCDWD_GlobalTileMapHoriz_Website_865x2250.jpg?VersionId=lOQ_j47U8T.j7UUqibD7SM63FVHjM_V5)

<img src="https://www.earthdata.nasa.gov/s3fs-public/2023-01/MCDWD_GlobalTileMapHoriz_Website_865x2250.jpg?VersionId=lOQ_j47U8T.j7UUqibD7SM63FVHjM_V5" alt="A hand in water" width="100%"/>

Based on availability, edit the tile_code variable:

```{r}
#add tile code from the map above (written as h00v00)
tile_code <- 'h05v05'
```

This is the NRT Flood F3 (MCDWD_L3_F3) API URL:

```{r}
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:

```{r}
#pasting together URL and year_day
url <- paste0(API_URL, year_day)
print(url)
```

### Load Data

Access the NASA Earthdata with the GET function:

```{r}
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)
    
    stars::write_stars(raster_df, "modis_nrt_flood.nc")
   
  } 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:

```{r}

raster_df <- stars::read_stars(download_link)
```

Set the Coordinate reference system (CRS) to "EPSG:4326"

```{r}
my_raster <- st_set_crs(raster_df, st_crs("EPSG:4326"))
st_crs(my_raster)
```

### Visualizing NRT Flood Data

Plot the raster to quickly view it:

```{r}
plot(my_raster, axes = TRUE)
```

#### Create NRT Flood Classification

Refer to the [MODIS NRT Global Flood Product User Guide](https://www.earthdata.nasa.gov/s3fs-public/2023-01/MCDWD_UserGuide_RevC.pdf) for more information.

NRT Flood data has 5 classifications:

| Code | Definition          |
|:----:|:--------------------|
|  0   | No Water            |
|  1   | Surface Water       |
|  2   | Recurring flood[^1] |
|  3   | Flood (unusual)     |
| 255  | Insufficient data   |

[^1]: Value 2 (Recurring flood) is not populated in the beta release.

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:

```{r}
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:

```{r}
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

```{r}
tmap_mode("view")
## 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)
  


```