Choosing recording locations from modeled sound data and urban land cover?

---

For noise, I downloaded the noise raster data from NPS (link). The data are based on this study. The data I pulled is the impact level – which is defined as the difference between existing sound levels and natural sound levels (which is why the values were so much lower than you expected).

For urban classification, I used the NLCD classified land cover from MRLC, downloaded using their viewer tool (link) medium- and high-intensity development for the NLCD classified land cover data.

---

Process raster data

1. Read raster files into R

library(tidyverse)

noise_big <-
  raster::raster(
    'CONUS_sumDay_L50dBA_imp.tiff')

nlcd_big <-
  raster::raster(
    'NLCD_CNBBEbkS1lK4LBmchsNP/NLCD_2016_Land_Cover_L48_20190424_CNBBEbkS1lK4LBmchsNP.tiff')

nlcd_legend <-
  read_csv('NLCD_CNBBEbkS1lK4LBmchsNP/NLCD_landcover_legend_2018_12_17_CNBBEbkS1lK4LBmchsNP.csv') %>%
  dplyr::filter(!is.na(Legend))

2. Define a regional subset by pulling the inner metro DC counties. Note: I figured counties would be the way to go about this. My first thought was to use all counties in the DMV, but most are pretty rural, so I though we might choose just the inner core counties of the DC metro region. Here’s how they look:

county_names_dmv <-
  c("District of Columbia",
    "Montgomery",
    "Prince George's",
    "Alexandria",
    "Arlington",
    "Fairfax County",
    "Fairfax",
    "Falls Church")

counties_dmv <-
  purrr::map_dfr(
    c('Virginia', 'Maryland', 'District of Columbia'),
    ~ tigris::counties(state = ., cb = TRUE)) %>% 
  dplyr::filter(NAME %in% county_names_dmv) %>% 
  sf::st_transform(crs = "+init=epsg:4326") %>%
  # Removing a Montgomery County in Virginia:
  dplyr::filter(GEOID != 51121) %>% 
  dplyr::mutate(NAME = factor(NAME))

counties_sp <-
  sf::as_Spatial(counties_dmv)

county_pal <- 
  leaflet::colorFactor(
    viridis::viridis(8), 
    counties_dmv$NAME)

leaflet::leaflet(counties_dmv) %>%
  leaflet::addTiles() %>%
  leaflet::addPolygons(
    color = "#777777", 
    weight = 1.25, 
    opacity = 1, 
    fillOpacity = 0.4,
    fillColor = ~county_pal(NAME))

3. Crop rasters to the extent of the counties file (to reduce memory load):

noise_cropped <-
  noise_big %>% 
  raster::crop(
    counties_dmv %>% 
      sf::st_transform(projection(noise_big)) %>% 
      sf::as_Spatial() %>% 
    raster::extent()) %>% 
  raster::projectRaster(crs = '+init=epsg:4326') %>% 
  raster::trim()

nlcd_cropped <-
  nlcd_big %>% 
  raster::crop(
    counties_dmv %>% 
      sf::st_transform(projection(noise_big)) %>% 
      sf::as_Spatial() %>% 
    raster::extent()) %>% 
  raster::projectRaster(crs = '+init=epsg:4326') %>% 
  raster::trim()

4. Mask the rasters to the county file:

noise_masked <-
  raster::mask(
    noise_cropped, 
    sf::as_Spatial(counties_dmv))

nlcd_masked <-
  raster::mask(
    nlcd_cropped, 
    sf::as_Spatial(counties_dmv))

5. Resample the nlcd raster (30m resolution) to match the origin and resolution of noise raster (270m):

nlcd_resampled <-
  raster::resample(
    x = nlcd_masked,
    y = noise_masked,
    method = 'ngb')

Here’s what the noise raster looks like:

noise_palette <- 
  leaflet::colorNumeric(
    c("#0000ff", "#ffff00", "#ff0000"), 
    raster::values(noise_masked),
  na.color = "transparent")

leaflet::leaflet(counties_dmv) %>%
  leaflet::addTiles() %>%
  leaflet::addRasterImage(
    noise_resampled,
    colors = noise_palette,
    opacity = 0.85) %>%
  leaflet::addLegend(
    pal = noise_palette,
    values = raster::values(noise_masked),
    title = "Anthropogenic noise") %>% 
  leaflet::addPolygons(
    color = "#777777", 
    weight = 1.25, 
    opacity = 1, 
    fillOpacity = 0)

---

Binary rasters: Noisy vs. loud

I thought a good way to establish cutoffs for the noise data might be the quartiles of the distribution.

tibble(noise = raster::values(noise_masked)) %>%
  na.omit() %>% 
  ggplot2:::ggplot(aes(x = noise)) +
  ggplot2::geom_histogram(fill = "#dcdcdc", color = '#000000', bins = 50) +
  ggplot2::theme_bw()

I’m questioning that now:

quiet_cutoff <-
  quantile(noise_masked, .25)

loud_cutoff <-
  quantile(noise_masked, .75)

quiet_cutoff

     25% 
10.73598 

loud_cutoff

     75% 
13.83388 

tibble(noise = raster::values(noise_masked)) %>%
  na.omit() %>% 
  ggplot2:::ggplot(aes(x = noise)) +
  ggplot2::geom_histogram(fill = "#dcdcdc", color = '#000000', bins = 50) +
  ggplot2::geom_vline(xintercept = quiet_cutoff, color = '#0000ff') +
  ggplot2::geom_vline(xintercept = loud_cutoff, color = '#ff0000') +
  ggplot2::theme_bw()

This yields the following areas for quiet sites:

noise_quiet <-
  raster::reclassify(
    noise_masked,
    rcl = 
      matrix(
        c(0, quiet_cutoff, 1, quiet_cutoff, Inf, NA),
        byrow = TRUE,
        nrow = 2),
    include.lowest = TRUE)

quiet_palette <- 
  leaflet::colorNumeric(
    c("#0000ff"), 
    raster::values(noise_quiet),
  na.color = "transparent")

leaflet::leaflet(counties_dmv) %>%
  leaflet::addTiles() %>%
  leaflet::addRasterImage(
    noise_quiet,
    colors = quiet_palette,
    opacity = 0.5) %>%
  leaflet::addPolygons(
    color = "#777777", 
    weight = 1.25, 
    opacity = 1, 
    fillOpacity = 0)

And these areas for loud sites:

noise_loud <-
  raster::reclassify(
    noise_masked,
    rcl = 
      matrix(
        c(0, loud_cutoff, NA, loud_cutoff, Inf, 1),
        byrow = TRUE,
        nrow = 2))

loud_palette <- 
  leaflet::colorNumeric(
    c("#ff0000"), 
    raster::values(noise_loud),
  na.color = "transparent")

leaflet::leaflet(counties_dmv) %>%
  leaflet::addTiles() %>%
  leaflet::addRasterImage(
    noise_loud,
    colors = loud_palette,
    opacity = 0.5) %>%
  leaflet::addPolygons(
    color = "#777777", 
    weight = 1.25, 
    opacity = 1, 
    fillOpacity = 0)

Maybe it’s okay? Hard to say … the spatial distributions look pretty doable, and reasonable for sampling, but is such a small difference biologically meaningful? Maybe there’s a better cutoff than quartiles? For example, here’s what the upper and lower 10% look like:

q10 <-
  quantile(noise_masked, .1)

q90 <-
  quantile(noise_masked, .9)

tibble(noise = raster::values(noise_masked)) %>%
  na.omit() %>% 
  ggplot2:::ggplot(aes(x = noise)) +
  ggplot2::geom_histogram(fill = "#dcdcdc", color = '#000000', bins = 50) +
  ggplot2::geom_vline(xintercept = q10, color = '#0000ff') +
  ggplot2::geom_vline(xintercept = q90, color = '#ff0000') +
  ggplot2::theme_bw()

And here’s the quiet areas with a 10% cutoff:

noise_quiet10 <-
  raster::reclassify(
    noise_masked,
    rcl = 
      matrix(
        c(0, q10, 1, q10, Inf, NA),
        byrow = TRUE,
        nrow = 2),
    include.lowest = TRUE)

quiet_palette10 <- 
  leaflet::colorNumeric(
    c("#0000ff"), 
    raster::values(noise_quiet10),
  na.color = "transparent")

leaflet::leaflet(counties_dmv) %>%
  leaflet::addTiles() %>%
  leaflet::addRasterImage(
    noise_quiet10,
    colors = quiet_palette10,
    opacity = 0.5) %>%
  leaflet::addPolygons(
    color = "#777777", 
    weight = 1.25, 
    opacity = 1, 
    fillOpacity = 0)

And loud areas with a 10% cutoff:

noise_loud90 <-
  raster::reclassify(
    noise_masked,
    rcl = 
      matrix(
        c(0, q90, NA, q90, Inf, 1),
        byrow = TRUE,
        nrow = 2))

loud_palette90 <- 
  leaflet::colorNumeric(
    c("#ff0000"), 
    raster::values(noise_loud90),
  na.color = "transparent")

leaflet::leaflet(counties_dmv) %>%
  leaflet::addTiles() %>%
  leaflet::addRasterImage(
    noise_loud90,
    colors = loud_palette90,
    opacity = 0.5) %>%
  leaflet::addPolygons(
    color = "#777777", 
    weight = 1.25, 
    opacity = 1, 
    fillOpacity = 0)

Long story short: Not many options for quiet sites within a reasonable driving distance at the 10% cutoff. Loud areas look reasonable. My house is classified within the upper 10% and I can verify that it’s louder than heck at my house, though mostly it’s the damned ice cream trucks and staging for weekend street races – not sure if they included those in their model. At the 25% cutoff, lots of options – but, whether it’s biologically meaningful remains a question.

Binary rasters: Urban vs. rural

Using the cut-off of medium and high-intensity development, I reclassified urban land cover as “1” and aggregated cells (aggregation factor = 7, function = max). Here’s the geographic distribution of urban land cover within the study area:

urban_rc <- nlcd_resampled

urban_rc[!values(urban_rc) %in% c(23, 24)] <- NA

urban_rc[values(urban_rc) %in% c(23, 24)] <- 1

urban_palette <-
  leaflet::colorNumeric(
    c("#ff0000"),
    raster::values(urban_rc),
  na.color = "transparent")

leaflet::leaflet(counties_dmv) %>%
  leaflet::addTiles() %>%
  leaflet::addRasterImage(
    urban_rc,
    colors = urban_palette,
    opacity = 0.75) %>%
  leaflet::addPolygons(
    color = "#777777",
    weight = 1.25,
    opacity = 1,
    fillOpacity = 0)

NA

One thing that I would say with the result (if you zoom in, especially to suburban areas) is that the nlcd does pretty poorly at predicting urban development, but it’s a starting point.

… and here’s the geographic distribution of non-urban land cover within the study area (unsurprisingly, almost everything):

nonurban_rc <- nlcd_resampled

nonurban_rc[values(nonurban_rc) %in% c(0, 11, 23, 24)] <- NA

nonurban_rc[!is.na(values(nonurban_rc))] <- 1

nonurban_palette <-
  leaflet::colorNumeric(
    c("#0000ff"),
    raster::values(nonurban_rc),
  na.color = "transparent")

leaflet::leaflet(counties_dmv) %>%
  leaflet::addTiles() %>%
  leaflet::addRasterImage(
    nonurban_rc,
    colors = nonurban_palette,
    opacity = 0.75) %>%
  leaflet::addPolygons(
    color = "#777777",
    weight = 1.25,
    opacity = 1,
    fillOpacity = 0)

---

Combining noisy, quiet, urban and not-urban:

Just a bit of raster math for this one. I’m back to using the quartile cutoff for noise.

We’ll start with the easy one. Here’s noisy urban sites:

leaflet::leaflet(counties_dmv) %>%
  leaflet::addTiles() %>%
  leaflet::addRasterImage(
    urban_rc*noise_loud,
    colors = urban_palette,
    opacity = 0.75) %>%
  leaflet::addPolygons(
    color = "#777777",
    weight = 1.25,
    opacity = 1,
    fillOpacity = 0)

Here’s quiet sites that are classified as urban:

urban_quiet <-
  urban_rc*noise_quiet

leaflet::leaflet(counties_dmv) %>%
  leaflet::addTiles() %>%
  leaflet::addRasterImage(
    urban_quiet,
    colors = urban_palette,
    opacity = 0.75) %>%
  leaflet::addPolygons(
    color = "#777777",
    weight = 1.25,
    opacity = 1,
    fillOpacity = 0)

The problem with the urban classification really shines here. If you explore around, you can see those few red pixels in the mix. What you can’t see them Shawn? I’ve made them into points to better show off the results:

urban_quiet_points <-
  raster::rasterToPoints(urban_quiet, spatial = TRUE)

leaflet::leaflet(counties_dmv) %>%
  leaflet::addTiles() %>%
  leaflet::addRasterImage(
    urban_quiet,
    colors = urban_palette,
    opacity = 0) %>%
  leaflet::addMarkers(urban_quiet_points@coords[,1], urban_quiet_points@coords[,2]) %>% 
  leaflet::addPolygons(
    color = "#777777",
    weight = 1.25,
    opacity = 1,
    fillOpacity = 0)

I’d be hard pressed to call any of these places urban (and you’d have to drive all over hell to get to them). Of course, there’s also the issue to the noise map. Because urban land cover / roads have such a big impact on modeled noise, it isn’t likely going to predict much in the way of quiet urban sites. This may be related to the truth on the ground and limitations of the model.

You can get to quite a few “non-urban” areas that are loud:

leaflet::leaflet(counties_dmv) %>%
  leaflet::addTiles() %>%
  leaflet::addRasterImage(
    nonurban_rc*noise_loud,
    colors = urban_palette,
    opacity = 0.75) %>%
  leaflet::addPolygons(
    color = "#777777",
    weight = 1.25,
    opacity = 1,
    fillOpacity = 0)

And, of course, quite a bit of “non-urban” areas that are quiet:

leaflet::leaflet(counties_dmv) %>%
  leaflet::addTiles() %>%
  leaflet::addRasterImage(
    nonurban_rc*noise_quiet,
    colors = urban_palette,
    opacity = 0.75) %>%
  leaflet::addPolygons(
    color = "#777777",
    weight = 1.25,
    opacity = 1,
    fillOpacity = 0)

Despite the limitations, if you decide to use the urban v. rural and loud v. quiet layers to select your sampling locations, it’ll be pretty easy to lay some points on these surfaces.