Build Sea Level Rise census Tracts and Dasymetric Tracts

Create rasters by census tract with sea level rise and dasymetric data in them.

Sea level rise and population distribution

For each census tract, create a series of sea level rise grids for 0-10 feet of sea level rise. On the same grid coordinates, cut the dasymetric grid for the lower 48.

The grids will be saved at 90 meter grid spacing, converted to WGS84 and lat-long coordinates. The water depths will be in feet. Population estimates for a grid cell less than 0.1 person will be set to NA.

Sea level rise data comes from https://coast.noaa.gov/slrdata/Depth_Rasters/index.html

To download the files, a convenient file of URL’s is supplied: wget –input-file=URLlist_SC.txt

Dasymetric grid comes from https://catalog.data.gov/dataset/enviroatlas-2020-dasymetric-population-for-the-conterminous-united-states-v1-12

Tract boundaries come from the 2020 census data via tidycensus.

Process outline

• Minimize coordinate projections
• Attach dasym
• Get tracts for a state
• project tracts to googlecrs
• read in a tile
• resample tile to 30 m (factor=9)
• project tile to googlecrs
• for each SL tile, find intersecting tracts
• get extent of SL tile
• project extent to Albers
• cut out dasym
• project cutout to googlecrs and resample to tile grid
• cut SL tile for each intersecting tract and save
• cut dasym tile by tract and save

Testing stub

# #| eval: false
#| warning: false
#| echo: false
#| message: false

##############################
#   Attach dasymetric file and get the geodetics for it
##############################

State <- "LA"
file <- "2020_Dasymetric_Population_CONUS_V1_1.tif"

x <- terra::rast(paste0(Dasym_path, "2020_Dasymetric_Population_CONUS_V1_1.tif")) 

##############################
# - read in a tile and resample
##############################

#   Read in tile and change grid spacing to 30 m from 10 m
SL_orig <- terra::aggregate(
              terra::rast(paste0(SL_path, "LA_CentralEast_slr_depth_3ft.tif")), 
                         fact=3, fun="median", na.rm=TRUE)
#   Get rid of flags
SL_orig <- terra::ifel(SL_orig > 9, NA, SL_orig)
#   convert depths to feet
SL_orig <- SL_orig*3.28084

##############################
#   Project dasym file onto SL coordinates
##############################

Dasym <- terra::project(x, SL_orig, method="sum")

##############################
#   Make a 2-layer file and save
##############################

Combo <- c(Dasym, SL_orig) %>% 
    tidyterra::rename(Pop='2020_Dasymetric_Population_CONUS_V1_1',
                      Depth='LA_CentralEast_slr_depth_3ft') 

# terra::writeRaster(Combo, paste0(Combo_path, file), overwrite=TRUE)

##############################
# - Get tracts for a state and convert to NAD83
##############################

Tracts <- readRDS(paste0(Census_path, "ACS_2023_", State, ".rds"))

Tracts <- Tracts %>% sf::st_transform(crs=terra::crs(Dasym)) %>% 
  filter(!sf::st_is_empty(.)) %>% #  remove empty collections
  sf::st_make_valid()

#   Make some plots

# tmap::tmap_options(basemaps="OpenStreetMap")
tmap_options(basemap.server = c("OpenStreetMap"))

tmap_mode <- tmap::tmap_mode("view") # set mode to interactive plots

tmap::tm_shape(terra::ifel(Dasym > 10, 10, Dasym)) +
  tmap::tm_raster(col_alpha=0.7,
                  # palette = "brewer.reds",
                  col.scale = tm_scale(values="brewer.reds"),
                  col.legend = tm_legend(title="Dasym") )

                  # title = "Dasym")

tmap::tm_shape(terra::as.polygons(terra::ext(Combo))) +
  tmap::tm_polygons(fill_alpha=0.0, col="blue", lwd=3) 

tmap::tm_shape(SL_orig) +
  tmap::tm_raster(col_alpha=0.7,
                  # palette = "brewer.blues",
                  col.scale = tm_scale(values="brewer.blues"),
                  col.legend = tm_legend(title="SL_orig") )

                  # title = "SL_orig")

##############################
# - for each SL tile, find intersecting tracts
##############################

#     Pull out only tracts that intersect the spatraster

terra::gdalCache(size = 2000)

foo <- prioritizr::intersecting_units(Tracts, SL_orig)
Tract_trim <- Tracts[foo,]

##############################
#  - For each tract, calculate statistics
##############################

val_by_tract <- exactextractr::exact_extract(
                     progress=FALSE,
                     x=Combo,
                     y=Tract_trim,
                     include_cols=c("GEOID","PopE")) %>% 
                 dplyr::bind_rows() %>%
                 tibble::as_tibble()

Values_by_node <- val_by_tract %>% 
  filter(Pop>0) %>% #   drop nodes with zero population
  replace_na(list(Depth=-1)) %>% 
  mutate(Depth_bins=cut_width(Depth, width=0.5, center=0.25, closed="left")) 
  
Stats_by_tract <-  Values_by_node %>% 
  group_by(GEOID, Depth_bins) %>% 
    summarize(PopE=first(PopE),
              Pop=sum(Pop*coverage_fraction),
              Num=sum(coverage_fraction)
              )
    
##############################
#  Stats for entire tile, and plot
##############################

Tract_pops <- Stats_by_tract %>% 
  group_by(GEOID) %>% 
    summarise(Pop_tot=sum(Pop),
             Area_tot=sum(Num)
             ) %>% 
  ungroup()

Stats_by_tract <- Stats_by_tract %>% 
  left_join(., Tract_pops, by="GEOID") %>% 
  mutate(Pop_pct=Pop/Pop_tot,
         Area_pct=Num/Area_tot)  %>% 
  ungroup()

Stats_by_tract2 <- Stats_by_tract %>% 
  filter(Pop>=1) %>% 
  filter(Depth_bins!="[-1,-0.5)") %>% 
  group_by(Depth_bins) %>% 
    summarise(Affected=sum(Pop)
              ) %>% 
  ungroup()

Stats_by_tract2 %>% 
  ggplot(aes(x=Depth_bins, y=Affected)) +
  geom_bar(stat="identity") +
  labs(title="Number of people affected by flooding",
       subtitle="Sea Level Rise = 3 feet, LA Central East Tile",
       x="Inundation Depth Range (feet)",
       y="Population affected")

Let’s do Louisiana

First let’s look at the tile extents

It appears that they must overlap a great deal, since the actual data appears to follow parish/county boundaries. But let’s make sure that is the case.

The tile extents are defined by the county lines, so much of the area around the edges of the tiles is set to NA, so we must work with only the tracts entirely contained within the tile extents, to avoid getting bogus results.

For the map, I set the NA values equal to 6, so that they will appear on the map. The values in the sea level grid are meters of water depth, with values greater than 9 representing current water surfaces.

Thankfully, the map shows that the tiles were designed very well, such that when tile boundaries cut tract boundaries, the grid nodes within the partial tract are set to null.

# #| eval: false 
#| warning: false
#| echo: false
#| message: false

#   Modified to use resampled tiles
State <- "LA"
infiles <- list.files(path=Combo_path, pattern=paste0("^",State,"_.*\\.tif$")) 
infiles <- infiles[stringr::str_detect(infiles, "Depth")]
Basenames <- infiles %>% 
  stringr::str_remove(., "_Depth_.*tif") %>%  unique()
Depth_names <- as.character(0:20/2) %>% stringr::str_remove(., "\\.0$") %>% 
  stringr::str_replace(., "\\.", "_")
SL <-Depth_names[21] 

i <- 1
for (Tilename in Basenames){
  memories <- gc()
      #   Read in tile 
    SL_temp <-  terra::rast(paste0(Combo_path, Tilename, "_Depth_",SL ,".tif")) 
    
    if (i==1) {
      Extents <- terra::as.polygons(terra::ext(SL_temp), crs=terra::crs(SL_temp))
    } else {
      Extents <- rbind(Extents, 
                       terra::as.polygons(terra::ext(SL_temp),
                       crs=terra::crs(SL_temp)))
    }
    i <-  i + 1
}

memories <- gc() # free memory
SL_temp <- terra::aggregate(
              terra::rast(paste0(Combo_path, Basenames[1], "_Depth_",SL ,".tif")), 
                         fact=3, fun="median", na.rm=TRUE)
SL_water <- terra::ifel(SL_temp > 9, NA, SL_temp)
SL_temp <- terra::ifel(is.nan(SL_temp), 6, SL_temp)

hist(SL_water[!is.na(SL_water)], xlab="Water Depth (meters)",
     main="Water depth distribution for one tile of 10 foot SL rise")

hist(SL_temp[!is.na(SL_temp)], xlab="Water Depth (meters)",
     main="Water depth distribution for one tile of 10 foot SL rise")

#   Get tracts *contained* within tile extent

Tracts <- readRDS(paste0(Census_path, "ACS_2023_", State, ".rds"))
Extent <- sf::st_as_sf(terra::as.polygons(terra::ext(SL_temp), crs=terra::crs(SL_temp)))
Extent <- sf::st_transform(Extent, crs=sf::st_crs(Tracts))
foobar <- sf::st_within(Tracts, Extent, sparse=FALSE)
Tract_trim <- Tracts[foobar,]

tmap::tmap_options(basemap.server = c("OpenStreetMap"))

tmap::tmap_mode("view") # set mode to interactive plots

tmap::tm_shape(Extents) +
  tmap::tm_polygons(fill_alpha=0.0, col="blue", lwd=3)+
tmap::tm_shape(SL_temp) +
  tmap::tm_raster(col_alpha=0.7,
                  col.scale = tm_scale(values="brewer.blues"),
                  col.legend = tm_legend(title="Inches Flooding") ) +
                  # palette = "brewer.blues",
                  # title = "Inches flooding") +
tmap::tm_shape(Tracts) +
  tmap::tm_polygons(fill_alpha=0.0, col="pink", lwd=1) +
tmap::tm_shape(Tract_trim) +
  tmap::tm_polygons(fill_alpha=0.0, col="red", lwd=3) 

And here we do the real work

First lets resample the dasym file once and save it so we don’t pay that expense more than once.

#  Resample the dasymmetric file and save it

x <- terra::rast(paste0(Dasym_path, "2020_Dasymetric_Population_CONUS_V1_1.tif")) 

print("Read in dasym and resample")
#   Change grid spacing to 90 m from 30 m
x <- terra::aggregate(x, fact=3, fun="sum", na.rm=TRUE)


#   Save resampled file
terra::writeRaster(x, paste0(Dasym_path, 
                             "2020_Dasymetric_Population_CONUS_V1_1_resamp90m.tif"),
                   overwrite=TRUE)

Function to get max extent from a set of tiles

max_extent <- function(files, path) {
  # files <- list.files(path=SL_path, pattern=paste0("^",State,"_EKA_.*ft\\.tif$"))
  
  Tile_xmin <- 0
  Tile_xmax <- 0
  Tile_ymin <- 0
  Tile_ymax <- 0
  for (file in files) {
    Xtnt <- terra::ext(terra::rast(paste0(path, file)))
    # print(paste0(path, file))
    Tile_xmin <- min(Tile_xmin, Xtnt$xmin)
    Tile_xmax <- min(Tile_xmax, Xtnt$xmax)
    Tile_ymin <- max(Tile_ymin, Xtnt$ymin)
    Tile_ymax <- max(Tile_ymax, Xtnt$ymax)
    # print(paste(Tile_xmin, Tile_xmax, Tile_ymin, Tile_ymax))
  }
  
  return(terra::ext(c(Tile_xmin, Tile_xmax, Tile_ymin, Tile_ymax)))
}

For each state, save resampled tiles and a tract by tract summary

Especially for California, there is a lot of cruft to deal with. Some tiles have no water, the extents may not be consistent, so those cases must be dealt with.

[1] "=============================="
[1] "Tilename = DC"
[1] "=============================="

Fix CA files

Some of the California files have different extents, so we needed to test ways to deal with that. Also, two files (as shown below) were in projected coordinates instead of geographic, so we had to project them into geographic coordinates in order for the process to proceed.

State="CA"
files <- list.files(path=SL_path, pattern=paste0("^",State,"_MTR23_.*ft\\.tif$"))

Tile_xmax <- NULL
Tile_ymin <- NULL
Tile_ymax <- NULL
for (file in files) {
  print(file)
  print(terra::ext(terra::rast(paste0(SL_path, file))))
  Xtnt <- terra::ext(terra::rast(paste0(SL_path, file)))
  Tile_xmin <- min(Tile_xmin, Xtnt$xmin)
  Tile_xmax <- min(Tile_xmax, Xtnt$xmax)
  Tile_ymin <- max(Tile_ymin, Xtnt$ymin)
  Tile_ymax <- max(Tile_ymax, Xtnt$ymax)
}

Max_size <- terra::ext(c(Tile_xmin, Tile_xmax, Tile_ymin, Tile_ymax))


##############   fix rasters that were projected to XY

Template <- terra::rast(paste0(SL_path, files[6]))
Old_rast <- terra::rast(paste0(SL_path, files[7]))

New_rast <- terra::project(Old_rast, terra::crs(Template))
terra::writeRaster(New_rast, paste0(SL_path, files[7]),
                     overwrite=TRUE)                        

Old_rast <- terra::rast(paste0(SL_path, files[8]))

New_rast <- terra::project(Old_rast, terra::crs(Template))
terra::writeRaster(New_rast, paste0(SL_path, files[8]),
                     overwrite=TRUE)                        

Fix VA files

Some of the Virginia are projected to XY coordinates instead of geographic, so we had to project them into geographic coordinates in order for the process to proceed. The Mid and N files in particular.

State="VA"
files <- list.files(path=SL_path, pattern=paste0("^",State,"_N.*ft\\.tif$"))
Template <- terra::rast(paste0(SL_path, files[1]))

##############   fix rasters that were projected to XY

for (file in files) {
  Old_rast <- terra::rast(paste0(SL_path, file))
  footest <- stringr::str_extract(terra::crs(Old_rast), "^.*\\[")
  print(paste(file, footest))
  if (footest=="GEOGCRS[") {
    next
  } else {
    New_rast <- terra::project(Old_rast, terra::crs(Template))
    terra::writeRaster(New_rast, paste0(SL_path, file),
                         overwrite=TRUE)                        
  }
}

Check states for redo

Ugh. They all need to be redone to fix depth bins.

for (state in c("AL", "CT", "DC", "DE", "MA", "MS", "NH", "PA", "RI", "WA", "LA")) {
  print("======================")
  print(paste("State:", state))
  print("======================")
  files <- list.files(path=Combo_path, pattern=paste0("^",state,".*\\.rds$"))
  for (file in files){
    print(paste("------->", file))
    foobar <- readRDS(paste0(Combo_path, file))
    print(unique(foobar$Depth_bins))
  }
}