In choisy/srtmVN: Elevations Data for Vietnam from SRTM90

knitr::knit_hooks$set(
  margin = function(before, options, envir) {
    if (before) par(mgp = c(1.5, .5, 0), bty = "n", plt = c(.105, .97, .13, .97))
    else NULL
  },
  prompt = function(before, options, envir) {
    options(prompt = if (options$engine %in% c("sh", "bash")) "$ " else "> ")
  })

knitr::opts_chunk$set(cache = TRUE, autodep = TRUE, message = FALSE, warning = FALSE,
                      margin = TRUE, dev.args = list(pointsize = 11), fig.height = 3.5,
                      fig.width = 4.24725, fig.retina = 2, fig.align = "center",
                      eval = FALSE)

options(width = 137)

Introduction

This script shows how to generate average elevations by province. We consider 2 measures of average elevation. The first one, elevation1 is the simple mean elevation of all the pixels inside a given province. The second one, elevation2 is an average of the elevation, weighted by local population density, in these same pixels. Here we show how to generate this data frame.

Packages

library(sp)
library(raster)
library(sptools)

Data sets

Human population density raster:

popdensity <- worldpopVN::getpop()

Not projected, with 0.000833 resolution:

popdensity

Elevation raster

srtm <- srtmVN::getsrtm()

Not projected, with 0.000833 resolution:

srtm

The provinces polygons:

provinces <- sf::as_Spatial(gadmVN::gadm())

Not projected:

proj4string(provinces)

Making sure that the 2 rasters are on the same grid

We need to resample srtm on popdensity (takes about 2'):

srtm2 <- resample(srtm, popdensity)

we can now compare srtm2 with popdensity:

srtm2

popdensity

Splitting the rasters by provinces

Next, we need to split these 2 rasters by the provinces of the polygons (takes about 20' each):

system.time(srtm_split <- split_on_poly(srtm2, provinces))
saveRDS(srtm_split, "srtm_split.rds")

and:

popdensity_split <- split_on_poly(popdensity, provinces)

Calculating population weights

Let's consider the following scaling function:

scaling <- function(x) {
  v <- values(x)
  values(x) <- v / sum(v, na.rm = TRUE)
  x
}

And let's run it on the population raster of each province:

popdensity_split2 <- purrr::map(popdensity_split, scaling)

Putting together and saving to disk

srtm_prov <- data.frame(
  province = provinces$province,
  elevation1 = unlist(purrr::map(srtm_split, ~ mean(values(.), na.rm = TRUE))),
  elevation2 = unlist(purrr::map2(popdensity_split2, srtm_split, ~ sum(values(.x) * values(.y), na.rm = TRUE)))
)

usethis::use_data(srtm_prov, overwrite = TRUE)

Comparing the 2 measures of elevation

plot(elevation2 ~ elevation1, srtm_prov, xlab = "weighted by population density", ylab = "non weighted")
abline(0, 1)

choisy/srtmVN documentation built on Aug. 21, 2019, 2:14 p.m.