First, let’s load necessary R packages

library(rgdal)
library(sf)
library(leaflet)
library(dplyr)
library(raster)
library(leafem)

This example shows basic flood analysis for a zip code in Florida
Data sets are downloaded from the following sources:
1. Florida Land Use: https://geodata.dep.state.fl.us/datasets/
2. Florida Zip Codes: https://download.fgdl.org/pub/state/zcta2000.zip
3. Florida Broward county DEM from EOX DAPA implementation (hosted on Euro Data Cube): https://dapa-eval-gmu.f874d3d9-442c-4b31-9b0a-de84208c3463.hub.eox.at/oapi/collections/DEM/dapa/area?bbox=-80.4765,25.9567,-80.0153,26.3555&fields=DEM
Then, let’s load the land use and zip code shapefiles

fl_lu <- st_read("D:/OGC16/Florida/Statewide_Land_Use_Land_Cover-shp/Statewide_Land_Use_Land_Cover.shp", 
        layer = "Statewide_Land_Use_Land_Cover")

fl_zip <- st_read("D:/OGC16/Florida/zcta2000/zcta2000.shp",
                  layer = "zcta2000")

fl_dem <- raster("https://dapa-eval-gmu.f874d3d9-442c-4b31-9b0a-de84208c3463.hub.eox.at/oapi/collections/DEM/dapa/area?bbox=-80.4765,25.9567,-80.0153,26.3555&fields=DEM")

Subset the shapefile with polygons to get the boundary for zip code 33314 and transform it to 4326

fl_33314 <- fl_zip[fl_zip$ZCTA5 %in% c("33314"), ]
fl_33314_4326 <- st_transform(fl_33314, st_crs(fl_lu))

Filter out polygons for specific land use (single family homes)

single_family_homes <- c("1110: Low Density, Fixed Single Family Units",
                         "1210: Medium Density, Fixed Single Family Units",
                         "1310: High Density, Fixed Single Family Units")
fl_lu_sfh <- fl_lu[fl_lu$LANDUSE_DE %in% single_family_homes, ]

Crop state DEM raster and single family homes land use by extent of zip code

fl_dem_33314 <- crop(fl_dem, extent(fl_33314_4326))
fl_lu_sfh_33314 <- st_intersection(fl_33314_4326, fl_lu_sfh)
fl_3314_dem_sfh <- mask(fl_dem_33314, fl_lu_sfh_33314)

Let’s take a look at the distribution of pixel values

hist(fl_3314_dem_sfh,
     main = "Distribution of elevation values in ZIP code 33314\n Single family homes land use",
     xlab = "Height (m)", ylab = "Number of Pixels",
     col = "springgreen",
     xlim=c(-2, 5),
     breaks=100)

Reclassify single family homes raster to identify areas that would be affected if water rises 1,2,3,4 feet

# 1 foot = 0.3048 m
reclass_vector <- c(-0.3048*2, -0.3048, -2,
                -0.3048, 0, -1,
                0, 0.3048, 1,
                0.3048, 0.3048*2, 2,
                0.3048*2, Inf, 3
)

Reshape the object into a matrix with columns and rows

reclass_matrix <- matrix(reclass_vector,
                    ncol = 3,
                    byrow = TRUE)

Reclassify the raster using the reclass object - reclass_matrix

fl_3314_dem_sfh_classified <- reclassify(fl_3314_dem_sfh, reclass_matrix)

Let’s plot the reclassified raster using leaflet for R

pal <- colorNumeric(c("red", "green", "blue"), 
                    values(fl_3314_dem_sfh_classified),
                    na.color = "transparent")
leaflet() %>%
        # Basemaps
        addTiles(group = "OSM (default)") %>%
        addProviderTiles(providers$Stamen.Toner, group = "Toner") %>%
        addProviderTiles(providers$Esri.NatGeoWorldMap, group = "NatGeo World Map") %>%
        addProviderTiles(providers$Esri.WorldImagery, group = "Satellite") %>%
        #Layers
        addRasterImage(fl_3314_dem_sfh_classified, opacity = 1, 
                       layerId = "Elevation", 
                       group = "Elevation")%>%
        #Pop up showing layer value
        addImageQuery(fl_3314_dem_sfh_classified, project = TRUE,
              layerId = "Elevation",position="bottomleft",suffix = " ft")%>%
        #Layer switcher
        addLayersControl(
                baseGroups = c("OSM (default)", "Satellite", "Toner", "NatGeo World Map"),
                overlayGroups = c("Elevation"),
                options = layersControlOptions(collapsed = FALSE)) %>%
        #Legend
        addLegend("bottomright", pal = pal, values = values(fl_3314_dem_sfh_classified),
                  title = "Elevation (in feet)") %>%
        addMouseCoordinates()