How to reproduce sea level rise code?

[Robinlovelace]

This reprex was created from https://github.com/geocompr/geostats_18/blob/master/code/geocompr/slr.R - using the reprex::reprex() function.

if(!require(elevatr)) {
  install.packages("elevatr")
}
#> Loading required package: elevatr
library(elevatr)
library(raster)
#> Loading required package: sp
library(sf)
#> Linking to GEOS 3.5.1, GDAL 2.2.2, proj.4 4.9.2

# load data
if(!dir.exists("data")) {

  download.file("https://github.com/geocompr/geostats_18/releases/download/0.1/data.zip", "data.zip")
  unzip("data.zip")
  file.rename("pres/geocompr/data/", "data")
}

# vector data
plot(spData::world)
#> Warning: plotting the first 9 out of 10 attributes; use max.plot = 10 to
#> plot all

poland = dplyr::filter(spData::world, grepl(pattern = "Pola", name_long))
plot(st_geometry(poland))

# testing elevatr (from readme / vignette)
# data(lake)
# elevation_high = get_elev_raster(lake, z = 14, src = "aws")
# 120 MB"
# elevation = get_elev_raster(lake, z = 9, src = "aws")
# elevation_low = get_elev_raster(lake, z = 4, src = "aws")
# pryr::object_size(elevation)
# pryr::object_size(elevation_low)
# plot(elevation_low)
# mapview::mapview(elevation_low)
# cols = rainbow(5)
# plot(elevation, add = T, col = cols)

# getting a map of the szczecin coastline
# option 1 - via osmdata:
# library(osmdata)
# scz = getbb(place_name = "szczecin", format_out = "sf_polygon")
# sf::st_write(scz, "data/scz-osm.geojson")
scz = sf::read_sf("data/scz-osm.geojson")
(m = mapview::mapview(scz)) # not quite what I was after, but a good start

# option 2 - via (new+improved) mapedit
# scz_mapedit = mapedit::drawFeatures(map = m)
# sf::write_sf(scz_mapedit, "data/scz-mapedit.geojson")
scz_mapedit = sf::read_sf("data/scz-mapedit.geojson")
plot(scz_mapedit)

scz_sp = as(scz_mapedit, "Spatial")
# e = elevatr::get_elev_raster(locations = scz_sp, z = 5, src = "aws")
# raster::writeRaster(x = e, filename = "data/scz-elev-z5.tif")
e = raster("data/scz-elev-z5.tif")

# get points less than 10m above slr
e_mask = e_low = e < 10
plot(e_low)
plot(spData::world, add = T, col = NA)
#> Warning in plot.sf(spData::world, add = T, col = NA): ignoring all but the
#> first attribute

e_mask[e > 10] = NA
e_low = mask(e, e_mask)
mapview::mapview(e_low)
# out-takes
# e = raster::getData(name = "SRTM", lon = c(-180, 180), lat = c(-60, 60))
writeRaster(e_low, "scz-elev-low.tif")
#> Error in .getGDALtransient(x, filename = filename, options = options, : filename exists; use overwrite=TRUE

system.time({
  r_orig = raster::raster("data/20180816184319_429977190-14-degrees-polish-coast.tif")
  plot(r_orig)
  })

#>    user  system elapsed 
#>   3.048   0.464   3.527
summary(r_orig)
#> Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (0.22% of all cells)
#>         X20180816184319_429977190.14.degrees.polish.coast
#> Min.                                                    0
#> 1st Qu.                                                 0
#> Median                                                  8
#> 3rd Qu.                                                35
#> Max.                                                  176
#> NA's                                                    0
r_agg = aggregate(r_orig, 100)
res(r_agg)
#> [1] 0.02777778 0.02777778
res(e)
#> [1] 0.0220 0.0134
plot(r_agg)
plot(e_low, add = T)

e_resampled = raster::resample(e, r_agg)
plot(values(r_agg), values(e_resampled))

cor(values(r_agg), values(e_resampled))^2
#> [1] 0.90463

# r = extend(r_agg, e)
# r_low = mask(r, e_mask)

# proportion of sczezcin that is at risk from slr
e_low_crop = crop(e_low, scz)
plot(e_low_crop)

e_low_sf = spex::polygonize(e_low_crop)
e_low_intersect = st_intersection(e_low_sf, scz)
#> although coordinates are longitude/latitude, st_intersection assumes that they are planar
plot(scz)
plot(e_low_sf, add = T)
plot(e_low_intersect, add = TRUE, col = "blue")

sum(st_area(e_low_intersect)) / st_area(scz)
#> 0.4800977 1

# publish
scz_buff = st_buffer(scz, dist = 0)
#> Warning in st_buffer.sfc(st_geometry(x), dist, nQuadSegs): st_buffer does
#> not correctly buffer longitude/latitude data
#> dist is assumed to be in decimal degrees (arc_degrees).
library(tmap)
scz$col = "red"
m = tm_shape(e_low) +
  tm_raster(alpha = 0.5) +
  tm_shape(scz_buff) + # not working...
  tm_fill(col = "red") +
  tm_shape(e_low_intersect) +
  tm_polygons() +
  tm_shape(poland) +
  tm_borders()
m
#> Variable "scz.elev.z5" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

Created on 2018-08-20 by the reprex package (v0.2.0).

[Robinlovelace]

geostat18-geocomputation.pdf

Those are the slides associated with this issue.

[Robinlovelace]

geocompr-preview.pdf