Climate change already impacts the Sahel region with an increase in the occurrence of extreme events such as floods. In Niger, these floods events usually happens along the Niger river. In Niamey, the heavy rains of the end of August and early September resulted in extreme flooding along the river path.

homes in Niger inundated by heavy flooding.Source: Getty Images

homes in Niger inundated by heavy flooding.Source: Getty Images

In this post, we will develop a reproducible workflow for a rapid impact assessment of flooding using Sentintel S2 data.

R packages required for the analysis

The following packages will be used for this analysis.

The packages rhdx and snapbox are not yet on CRAN, you will need the remotes package to get them.

remotes::install_github("dickoa/rhdx")
remotes::install_github("anthonynorth/snapbox")

This analysis is based on a UN-SPIDER tutorial on flood mapping using Sentintel S2 data. In the original article, QGIS was used to map flood extent, here, we are going to do the same analysis using R.

Area of interest

We will restrict this analysis to Niamey and its surrounding where the flood impacted many people lives. The bounding box of the study area can be downloaded here.

bbox_url <- "https://tinyurl.com/yyobbc58"
bbox <- read_sf(bbox_url)
ggplot() +
  layer_mapbox(bbox,
               map_style = mapbox_satellite_streets())

Permanent Water in Niamey

In order to map flood extent, we need to compare the waterbody mapped to a baseline of permanent water. JRC Global Surface Water is an excellent source to access surface water. Permanent water surfaces are defined using the water occurence layer which is the frequency with which water was present in a given pixel from March 1984 to December 2019. Water surfaces are permanent in pixel when it reaches the frequency 100. It means that for this pixel contained water surfaces over the 36 years timespan. In this analysis, the definition is relaxed and permanent water surfaces is defined by water with 80% occurence over 36 years. The data is available online where it can be downloaded.

occurence <- "https://storage.googleapis.com/global-surface-water/downloads2019v2/occurrence/occurrence_0E_20Nv1_1_2019.tif"
occ <- read_stars(occurence) %>%
  st_crop(bbox) %>%
  st_as_stars()
occ

#> stars object with 2 dimensions and 1 attribute
#> attribute(s), summary of first 1e+05 cells:
#>  occurrence_0E_20Nv1_1_2019.tif 
#>  1      :95865                  
#>  91     :  324                  
#>  90     :  299                  
#>  89     :  268                  
#>  88     :  151                  
#>  (Other): 2964                  
#>  NA's   :  129                  
#> dimension(s):
#>    from    to offset    delta refsys point values x/y
#> x  7889  9434      0  0.00025 WGS 84 FALSE   NULL [x]
#> y 25624 26866     20 -0.00025 WGS 84 FALSE   NULL [y]

Let’s transform the occ raster layer to build build our permanent water layer.

mat <- occ[[1]]
str(mat)

#>  Factor[1:1546, 1:1243] w/ 256 levels "1","2","3","4",..: NA 1 1 1 1 1 1 1 1 1 ...
#>  - attr(*, "colors")= chr [1:256] "#FFFFFFFF" "#FEFCFCFF" "#FEF9FAFF" "#FEF7F7FF" ...

The data is of factor type, to easily subset the data, we need to convert the type to numeric. Particularly, selecting pixels with an occurence of 80 will be achieved more efficiently with numeric vectors.

mat <- as.numeric(as.character(mat))
mat <- matrix(mat, nrow = dim(occ)[1])
mat[mat < 80] <- NA
mat[!is.na(mat)] <- 1L
occ$water_perm <- mat
occ

#> stars object with 2 dimensions and 2 attributes
#> attribute(s), summary of first 1e+05 cells:
#>  occurrence_0E_20Nv1_1_2019.tif   water_perm    
#>  1      :95865                   Min.   :1      
#>  91     :  324                   1st Qu.:1      
#>  90     :  299                   Median :1      
#>  89     :  268                   Mean   :1      
#>  88     :  151                   3rd Qu.:1      
#>  (Other): 2964                   Max.   :1      
#>  NA's   :  129                   NA's   :98221  
#> dimension(s):
#>    from    to offset    delta refsys point values x/y
#> x  7889  9434      0  0.00025 WGS 84 FALSE   NULL [x]
#> y 25624 26866     20 -0.00025 WGS 84 FALSE   NULL [y]

As a final step to have a permanent water layer, the raster layer sais need to turned into a spatial polygon for further analysis. In order to do this transformation, the st_as_sf function can be used to turn it the occ raster to a sf object.

permanent_water <- st_as_sf(occ["water_perm"], merge = TRUE) %>%
  filter(water_perm == 1) %>%
  st_union()

ggplot() +
  layer_mapbox(bbox,
               map_style = mapbox_satellite_streets())  +
  geom_sf(data = permanent_water, fill = "blue", size = 0.2)

Flood extent using Sentinel S2

A quick survey of the literature suggests that using the Normalized Difference Water Index (NDWI) index allow to easily spot waterbody content. In this analysis, ready to use bottom of atmosphere reflectances products (BOA) will also be used to make it simpler.

It’s defined by the following formulas:

$$NDWI = \dfrac{X_{green} - X_{nir}}{X_{green} + X_{nir}}$$

Where $X_{green}$ is the green band (B3 at at 842 nm wavelength) and $X_{nir}$ is the near-infrared (NIR) band (B8, at 842 nm wavelength).

We will use the amazing sen2r R package to download, pre-process and automatically compute the NDWI. In order to download Sentinel S2 data, you will need to register to the Copernicus Open Access Hub and provide your credentials to the sen2r::write_scihub_login function. This analysis will focus on the situation as of September 15th.

write_scihub_login(username = "xxxxx",
                   password = "xxxxx")


out_dir <- "./data/raster/s2out"
safe_dir <- "./data/raster/s2safe"

if (!dir.exists(out_dir))
  dir.create(out_dir, recursive = TRUE)

if (!dir.exists(safe_dir))
  dir.create(safe_dir, recursive = TRUE)

###
out_path <- sen2r(
  gui = FALSE,
  step_atmcorr = "l2a",
  extent = bbox,
  extent_name = "Niamey",
  timewindow = c(as.Date("2020-09-10"), as.Date("2020-09-15")),
  list_prods = "BOA",
  list_indices = "NDWI2",
  list_rgb = "RGB432B",
  path_l1c = safe_dir,
  path_l2a = safe_dir,
  path_out = out_dir,
  parallel = 8,
  overwrite = TRUE
)

There’s two versions of the NDWI index and the one we need is named NDWI2 in the sen2r package.

ndwi <- read_stars(file.path("data/raster/s2out",
                             "NDWI2",
                             "S2A2A_20200915_022_Niamey_NDWI2_10.tif"),
                   proxy = FALSE)
ndwi

#> stars object with 2 dimensions and 1 attribute
#> attribute(s), summary of first 1e+05 cells:
#>  S2A2A_20200915_022_Niamey_NDWI2_10.tif 
#>  Min.   :-7223                          
#>  1st Qu.:-4323                          
#>  Median :-3757                          
#>  Mean   :-3674                          
#>  3rd Qu.:-3067                          
#>  Max.   : 1271                          
#> dimension(s):
#>   from   to  offset delta                refsys point values x/y
#> x    1 4193  388660    10 WGS 84 / UTM zone 31N FALSE   NULL [x]
#> y    1 3450 1503090   -10 WGS 84 / UTM zone 31N FALSE   NULL [y]

Looking at the ndwi data range, we can rescale the data by dividing by 10000 to simplify the rest of the analysis and map it to visualize this layer.

ndwi[[1]] <- ndwi[[1]]/10000

ggplot() +
  geom_stars(data = ndwi) +
  scale_fill_viridis_c(guide = guide_legend(title = "NDWI",
                                            direction = "horizontal")) +
  theme(legend.position = "bottom")

Binarize the NDWI

The next step, is to turn this raster layer into a binary data. The theoretical threshold to select waterbody is set at 0 but we can be more precise and context-specific. We will thus sample NDWI values in pixels we know contains water surface using the permanent_water layer and average these values to get our data.

permanent_water <- st_transform(permanent_water, st_crs(ndwi))
ndwi_water <- st_crop(ndwi,
                      permanent_water)
me <- median(ndwi_water[[1]], na.rm = TRUE)
xbar <- mean(ndwi_water[[1]], na.rm = TRUE)
c(median = me, average = xbar)

#>     median    average 
#> 0.05240000 0.05036571

All sampled pixels within the permanent water layer can be visualized using an histogram in order its distribution and variability.

hist(ndwi_water[[1]], main = "", xlab = "NDWI")
abline(v = me, col = "red")
abline(v = xbar, col = "blue")

The distribution is unimodal and symmetric around the mode with low variability. Either the average or the median can be used as threshold. For the next step, the median will be used since it’s more robust than the average in general.

ndwi$water_extent <- 1L * (ndwi[[1]] >= me)
ndwi

#> stars object with 2 dimensions and 2 attributes
#> attribute(s), summary of first 1e+05 cells:
#>  S2A2A_20200915_022_Niamey_NDWI2_10.tif  water_extent   
#>  Min.   :-0.7223                         Min.   :0e+00  
#>  1st Qu.:-0.4323                         1st Qu.:0e+00  
#>  Median :-0.3757                         Median :0e+00  
#>  Mean   :-0.3674                         Mean   :7e-05  
#>  3rd Qu.:-0.3067                         3rd Qu.:0e+00  
#>  Max.   : 0.1271                         Max.   :1e+00  
#> dimension(s):
#>   from   to  offset delta                refsys point values x/y
#> x    1 4193  388660    10 WGS 84 / UTM zone 31N FALSE   NULL [x]
#> y    1 3450 1503090   -10 WGS 84 / UTM zone 31N FALSE   NULL [y]

For this data to be used it’s need to be polygonized, we will follow some steps (comments in the code) to make sure the output is good enough for our analysis:

  1. Polygonize the raster
  2. Select the polygons containing water surfaces
  3. Remove polygons smaller than 500 square meters (avoid artificial water bodies, wet rooftops, etc)
  4. Dissolve and join all polygons into a single polygon
water_extent <- ndwi["water_extent"] %>%
  ## Polygonize the raster
  st_as_sf(merge = TRUE) %>%
  ## select only water extent
  filter(water_extent == 1) %>%
  ## Remove polygons smaller than 500 sqm
  mutate(area = st_area(geometry)) %>% ## compute areas
  filter(area > set_units(500, m^2)) %>%
  ## Dissolve and join all polygons into a single polygon
  st_buffer(dist = set_units(0.5, m)) %>%
  st_union() %>%
  st_simplify(dTolerance = set_units(10, m)) %>%
  st_set_geometry(value = ., x = data.frame(id = 1))

water_extent

#> Simple feature collection with 1 feature and 1 field
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: 388659.5 ymin: 1468590 xmax: 430590.5 ymax: 1502860
#> projected CRS:  WGS 84 / UTM zone 31N
#>   id                       geometry
#> 1  1 MULTIPOLYGON (((430470.5 14...

The water_extent derived from NDWI can now be analyzed and visualized.

ggplot() +
  layer_mapbox(bbox,
               map_style = mapbox_satellite_streets())  +
  geom_sf(data = water_extent, fill = "blue", size = 0.2)

The final step required to obtain the flood extent layer is to remove the permanent water extent (permanent_water) from the waterbody derived from NDWI (water_extent).

flood_extent <- st_difference(water_extent,
                              permanent_water)


ggplot() +
  layer_mapbox(bbox,
               map_style = mapbox_satellite_streets())  +
  geom_sf(data = flood_extent, fill = "blue", size = 0.2)

Some statistics can be computed out of the flood_extent and one important metric is the total area of flooded zones within our area of interest.

flood_extent %>%
  st_area() %>%
  set_units(km^2) -> flood_area
flood_area

#> 44.84563 [km^2]

44.85 square kilometers of a total area of 1437 square kilometers is flooded as of September 15th.

A second figure of interest is the number of people potentially affected, to estimate this figure we need overlay our flood polygon (flood_extent) and some population count raster layer. In this analysis, we used the Facebook population density maps data for Niger. Facebook population density data are available on HDX and the rhdx can be used to access it.

dir <- dirname(tempdir())
pull_dataset("highresolutionpopulationdensitymaps-ner") %>%
  get_resource(1) %>%
  download_resource() %>%
  unzip(exdir = dir)

ner_pop <- read_stars(file.path(dir, "population_ner_2018-10-01.tif")) %>%
  st_crop(bbox)

names(ner_pop) <- "ner_pop"

ggplot() +
  geom_stars(data = ner_pop) +
  scale_fill_viridis_c(guide = guide_legend(title = "Population count 2018",
                                            direction = "horizontal")) +
  theme(legend.position = "bottom")

We now have the population grid, for each pixel of 30x30m it estimates the population count at that given pixel. To have the population affected, one can sum the values of all pixels falling into our flood extent layer.

ner_pop %>%
  st_as_stars() %>%
  st_crop(st_transform(flood_extent, 4326)) %>%
  pull(ner_pop) %>%
  sum(na.rm = TRUE) -> pop_affected
pop_affected

#> [1] 4512.976

We estimate around 4513 affected persons in Niamey and its surroundings. We can also estimate the number of affected buildings. A good source for buildings data is OpenStreetMap (OSM) and Niger has one vibrant OSM community, so the data is of good quality particularly in cities like Niamey. Once again, the data is available on HDX.

ner_bu <- pull_dataset("hotosm_niger_buildings") %>%
  get_resource(1) %>%
  read_resource()

ner_bu <- st_transform(ner_bu, st_crs(flood_extent))

ner_bu <- st_crop(ner_bu, flood_extent)
ner_bu_centroid <- st_point_on_surface(ner_bu)

impacted_bu <- st_intersection(ner_bu_centroid, flood_extent)
nrow(impacted_bu)

#> [1] 341

We estimated around 341 potentially affected building in the study area.

Limitations

The goal of this analysis was to show how to conduct a reproducible rapid flood mapping and damage assessment using R and Sentintel S2 data. We made several choices during the process and they all have an impact on the final results. For example, other spectral indices could have been used, some authors suggested that Modified NDWI (MDNWI) is more accurate in urban settings and others advocated for the use the Normalized Difference Flood Index for this type of analysis (NDFI). The package sen2r is quite flexible and allow to easily add all the above indices. For the impact and damage results, the population grid used in this analysis is based on 2018 population statistics. The analysis would give different results if other population estimate layers such as Worldpop population count grid were used. Finally, SAR imagery can also be used to push the analysis further and estimate for example indicators such as flood water depth for better planing of a response plan.

Session info for this analysis.

Session info
devtools::session_info()

#> ─ Session info ────────────────────────────────────────────────
#>  setting  value                                      
#>  version  R version 4.0.3 Patched (2020-10-10 r79320)
#>  os       Arch Linux                                 
#>  system   x86_64, linux-gnu                          
#>  ui       X11                                        
#>  language (EN)                                       
#>  collate  en_US.UTF-8                                
#>  ctype    en_US.UTF-8                                
#>  tz       UTC                                        
#>  date     2020-10-12                                 
#> 
#> ─ Packages ────────────────────────────────────────────────────
#>  package     * version     date       lib
#>  abind       * 1.4-5       2016-07-21 [1]
#>  assertthat    0.2.1       2019-03-21 [1]
#>  backports     1.1.10      2020-09-15 [1]
#>  base64enc     0.1-3       2015-07-28 [1]
#>  blob          1.2.1       2020-01-20 [1]
#>  brio          1.1.0       2020-08-31 [1]
#>  broom         0.7.1       2020-10-02 [1]
#>  callr         3.5.0       2020-10-08 [1]
#>  cellranger    1.1.0       2016-07-27 [1]
#>  class         7.3-17      2020-04-26 [1]
#>  classInt      0.4-3       2020-04-07 [1]
#>  cli           2.0.2.9000  2020-09-01 [1]
#>  codetools     0.2-16      2018-12-24 [1]
#>  colorspace    1.4-1       2019-03-18 [1]
#>  crayon        1.3.4.9000  2020-09-01 [1]
#>  crosstalk     1.1.0.1     2020-03-13 [1]
#>  crul          1.0.0       2020-07-30 [1]
#>  curl          4.3         2019-12-02 [1]
#>  data.table    1.13.0      2020-07-24 [1]
#>  DBI           1.1.0       2019-12-15 [1]
#>  dbplyr        1.4.4       2020-05-27 [1]
#>  desc          1.2.0.9000  2020-09-01 [1]
#>  devtools      2.3.2       2020-09-18 [1]
#>  digest        0.6.25      2020-02-23 [1]
#>  doParallel    1.0.15      2019-08-02 [1]
#>  downlit       0.2.0       2020-09-25 [1]
#>  dplyr       * 1.0.2       2020-08-18 [1]
#>  e1071         1.7-3       2019-11-26 [1]
#>  ellipsis      0.3.1.9000  2020-09-01 [1]
#>  evaluate      0.14.1      2020-09-01 [1]
#>  fansi         0.4.1       2020-01-08 [1]
#>  forcats     * 0.5.0       2020-03-01 [1]
#>  foreach       1.5.0       2020-03-30 [1]
#>  foreign       0.8-80      2020-05-24 [1]
#>  fs            1.5.0.9000  2020-09-01 [1]
#>  generics      0.0.2.9000  2020-09-03 [1]
#>  geojson       0.3.4       2020-06-23 [1]
#>  geojsonio     0.9.2.93    2020-09-02 [1]
#>  ggplot2     * 3.3.2       2020-06-19 [1]
#>  ggthemes    * 4.2.0       2019-05-13 [1]
#>  glue          1.4.2       2020-08-27 [1]
#>  gtable        0.3.0.9000  2020-09-01 [1]
#>  haven         2.3.1       2020-06-01 [1]
#>  hms           0.5.3       2020-01-08 [1]
#>  hoardr        0.5.2.93    2020-08-01 [1]
#>  htmltools     0.5.0       2020-06-16 [1]
#>  htmlwidgets   1.5.2       2020-10-03 [1]
#>  httpcode      0.3.0       2020-04-10 [1]
#>  httr          1.4.2.9000  2020-09-01 [1]
#>  hugodown      0.0.0.9000  2020-10-04 [1]
#>  iterators     1.0.12      2019-07-26 [1]
#>  jqr           1.1.0.92    2020-09-02 [1]
#>  jsonlite      1.7.1       2020-09-07 [1]
#>  KernSmooth    2.23-17     2020-04-26 [1]
#>  knitr         1.30        2020-09-22 [1]
#>  lattice       0.20-41     2020-04-02 [1]
#>  lazyeval      0.2.2       2019-03-15 [1]
#>  leafem        0.1.3       2020-07-26 [1]
#>  leaflet       2.0.3       2019-11-16 [1]
#>  lifecycle     0.2.0.9000  2020-09-01 [1]
#>  lubridate     1.7.9       2020-06-08 [1]
#>  lwgeom        0.2-5       2020-09-01 [1]
#>  magrittr      1.5.0.9000  2020-09-24 [1]
#>  maptools      1.0-2       2020-08-24 [1]
#>  mapview     * 2.9.3       2020-09-08 [1]
#>  memoise       1.1.0       2017-04-21 [1]
#>  modelr        0.1.8       2020-05-19 [1]
#>  munsell       0.5.0       2018-06-12 [1]
#>  pillar        1.4.6       2020-07-10 [1]
#>  pkgbuild      1.1.0.9000  2020-09-01 [1]
#>  pkgconfig     2.0.3       2019-09-22 [1]
#>  pkgload       1.1.0.9000  2020-09-01 [1]
#>  png           0.1-7       2013-12-03 [1]
#>  prettyunits   1.1.1       2020-01-24 [1]
#>  processx      3.4.4       2020-09-03 [1]
#>  ps            1.4.0       2020-10-07 [1]
#>  purrr       * 0.3.4       2020-04-17 [1]
#>  R6            2.4.1.9000  2020-09-01 [1]
#>  rappdirs      0.3.1       2016-03-28 [1]
#>  raster        3.3-13      2020-07-17 [1]
#>  Rcpp          1.0.5       2020-07-06 [1]
#>  RcppTOML      0.1.6       2019-06-25 [1]
#>  readr       * 1.4.0       2020-10-05 [1]
#>  readxl        1.3.1       2019-03-13 [1]
#>  remotes       2.2.0.9000  2020-09-01 [1]
#>  reprex        0.3.0       2019-05-16 [1]
#>  rgdal         1.5-17      2020-10-08 [1]
#>  rgeos         0.5-5       2020-09-07 [1]
#>  rhdx        * 0.1.0.9000  2020-09-03 [1]
#>  rlang         0.4.8       2020-10-08 [1]
#>  rmarkdown     2.4         2020-09-30 [1]
#>  rprojroot     1.3-2       2018-01-03 [1]
#>  rstudioapi    0.11        2020-02-07 [1]
#>  rvest         0.3.6       2020-07-25 [1]
#>  satellite     1.0.2       2019-12-09 [1]
#>  scales        1.1.1.9000  2020-09-01 [1]
#>  sen2r       * 1.3.8       2020-08-26 [1]
#>  sessioninfo   1.1.1.9000  2020-09-01 [1]
#>  sf          * 0.9-7       2020-10-10 [1]
#>  snapbox     * 0.2.1       2020-09-09 [1]
#>  sp            1.4-4       2020-10-07 [1]
#>  stars       * 0.4-4       2020-10-10 [1]
#>  stringi       1.5.3       2020-09-09 [1]
#>  stringr     * 1.4.0       2019-02-10 [1]
#>  stylebox      0.1.0       2020-09-09 [1]
#>  testthat      2.99.0.9000 2020-09-01 [1]
#>  tibble      * 3.0.3       2020-07-10 [1]
#>  tidyr       * 1.1.2       2020-08-27 [1]
#>  tidyselect    1.1.0       2020-05-11 [1]
#>  tidyverse   * 1.3.0       2019-11-21 [1]
#>  units       * 0.6-7       2020-06-13 [1]
#>  usethis       1.6.3       2020-09-17 [1]
#>  V8            3.2.0       2020-06-19 [1]
#>  vctrs         0.3.4.9000  2020-09-01 [1]
#>  webshot       0.5.2       2019-11-22 [1]
#>  withr         2.3.0       2020-09-22 [1]
#>  xfun          0.18        2020-09-29 [1]
#>  XML           3.99-0.5    2020-07-23 [1]
#>  xml2          1.3.2.9001  2020-09-01 [1]
#>  yaml          2.2.1       2020-02-01 [1]
#>  source                                
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.3)                        
#>  CRAN (R 4.0.3)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.3)                        
#>  CRAN (R 4.0.2)                        
#>  local                                 
#>  CRAN (R 4.0.3)                        
#>  CRAN (R 4.0.2)                        
#>  local                                 
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  local                                 
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  local                                 
#>  local                                 
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.3)                        
#>  local                                 
#>  local                                 
#>  CRAN (R 4.0.2)                        
#>  local                                 
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  local                                 
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  local                                 
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.3)                        
#>  CRAN (R 4.0.2)                        
#>  local                                 
#>  Github (r-lib/hugodown@18911fc)       
#>  CRAN (R 4.0.2)                        
#>  local                                 
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.3)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.3)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  local                                 
#>  CRAN (R 4.0.2)                        
#>  Github (r-spatial/lwgeom@99d07a7)     
#>  Github (tidyverse/magrittr@0221e18)   
#>  CRAN (R 4.0.2)                        
#>  Github (r-spatial/mapview@3b0810b)    
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  local                                 
#>  CRAN (R 4.0.2)                        
#>  local                                 
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.3)                        
#>  CRAN (R 4.0.2)                        
#>  local                                 
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.3)                        
#>  CRAN (R 4.0.2)                        
#>  local                                 
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.3)                        
#>  CRAN (R 4.0.2)                        
#>  local                                 
#>  CRAN (R 4.0.3)                        
#>  CRAN (R 4.0.3)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  local                                 
#>  CRAN (R 4.0.3)                        
#>  local                                 
#>  Github (r-spatial/sf@5caaca7)         
#>  Github (anthonynorth/snapbox@ba3141b) 
#>  CRAN (R 4.0.3)                        
#>  Github (r-spatial/stars@7f8ef33)      
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  Github (anthonynorth/stylebox@500218e)
#>  local                                 
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  local                                 
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  CRAN (R 4.0.2)                        
#>  local                                 
#>  CRAN (R 4.0.2)                        
#> 
#> [1] /usr/lib/R/library