Thursday with Ashok
Download the Rmd file to follow along!
PRESS CTRL+SHIFT+O in RStudio to show the outline
Preparation
Are you in the correct directory before we begin?
The following line calls a function from the here package, which uses a lot of heuristics to make sure that project-specific files are available and that any path we need in the following code chunks are available to use.
if(!require(here)) install.packages("here")## Loading required package: here## here() starts at /home/bryce/ashokThursday20240620here::i_am("ashokThursday20240620.Rmd")## here() starts at /home/bryce/ashokThursday20240620If the output of the chunk is like the following, you’re ready to continue.
here() starts at /home/bryce/ashokThursday20240620
Packages required by spatialEpisim
# An important and useful package which we need right away.
if (!"magrittr" %in% installed.packages())
install.packages("magrittr")
stopifnot(require(magrittr))## Loading required package: magrittrquote({
# A package symbol/name, i.e. unqouted, with no trailing commas.
av
bslib
colorRamps
countrycode
cptcity
deSolve
dplyr # Loaded by the tidyverse meta-package
DT
fasterize
geodata
# ggplot2 # Loaded by the tidyverse meta-package
ggmap
ggspatial
htmltools
latex2exp
lattice
latticeExtra
leaflet
lubridate
magick
maps
markdown
ncdf4
plotly
purrr
raster
rasterVis
Rcpp
readr
readxl
# rgdal
rjson
rosm
rsconnect
rstudioapi
sp
sf
shiny
shinyalert
shinybusy
shinyjs
shinyhelper
shinyWidgets
stringr
shinyFeedback
shinyvalidate
terra
tidyverse
tinytex
viridis
writexl
}) %>%
lapply(as.character) %>%
unlist() %>%
tail(-1) ->
requiredPackages
## The Magrittr pipe is used here so that users without R 4.0 don't have issues.
requiredPackagesInstalled <- (requiredPackages %in% installed.packages())
stopifnot(all(requiredPackagesInstalled))If the previous chunk caused a stop, then install any missing required packages.
if(!all(requiredPackagesInstalled))
install.packages(requiredPackages[!requiredPackagesInstalled])If you are running a “too recent” version of R, then the package rgdal isn’t available. For example, it isn’t available for R version 4.3.
Load all required packages, stopping if there are any issues loading any.
suppressAllButErrors <-
function(...)
suppressWarnings(suppressPackageStartupMessages(suppressMessages(...)))
## The Magrittr pipe is used here so that users without R 4.0 don't have issues.
lapply(requiredPackages, require, character.only = TRUE) %>%
unlist() %>%
all() %>%
stopifnot() %>%
suppressAllButErrors()You should update the required packages; if any were installed some time ago then there may be a bug which we have not encountered, but would be produced by old package versions.
update.packages(oldPkgs = requiredPackages,
repos = "https://cloud.r-project.org") %>%
suppressAllButErrors()Let’s begin!
Some utilitiy functions
You should ignore these functions definitions; understanding them is not important right now. If you’re reading in RStudio, use the Outline to skip to the Analysis section or collapse the function definitions and proceed to the following prose. If you’re using reading in the knitted HTML document, click that link or this one.
# With `directOutput = FALSE`, that means you don't want the plot
# to go directly to your RStudio viewer/plot pane, or if you're adventurous and
# not using RStudio then you should be able to guess what you don't want to
# happen. If it is `FALSE`, it means (to this code) that you want to save the
# plot to a file.
createBasePlot <- function(selectedCountry, susceptible, directOutput) {
inputISO <- countrycode(selectedCountry,
origin = 'country.name',
destination = 'iso3c')
inputISOLower <- tolower(inputISO)
# print(Susceptible)
# Only a terra::rast() object can use the classify() function
# Susceptible <- terra::rast(Susceptible)
fname <- sprintf("%s_PopulationCount.png", inputISO)
PNGFileName <- here("plots", fname)
if(!directOutput) png(PNGFileName, width = 1024, height = 768)
x <- classify(susceptible, c(0, 10, 25, 50, 100, 250, 1000, 100000))
# plot(x, col=pal(8)[-1], xlab = "Longitude", ylab = "Latitude")
levs <- levels(x)[[1]]
#levs[7] <- "> 1000"
levels(x) <- levs
# ramp <- c('#D0D8FB', '#BAC5F7', '#8FA1F1', '#617AEC', '#0027E0',
# '#1965F0', '#0C81F8', '#18AFFF', '#31BEFF', '#43CAFF',
# '#60E1F0', '#69EBE1', '#7BEBC8', '#8AECAE', '#ACF5A8',
# '#CDFFA2', '#DFF58D', '#F0EC78', '#F7D767', '#FFBD56',
# '#FFA044', '#EE4F4D')
ramp <- c('#FFFFFF', '#D0D8FB', '#BAC5F7', '#8FA1F1', '#617AEC',
'#0027E0', '#1965F0', '#0C81F8', '#18AFFF', '#31BEFF',
'#43CAFF', '#60E1F0', '#69EBE1', '#7BEBC8', '#8AECAE',
'#ACF5A8', '#CDFFA2', '#DFF58D', '#F0EC78', '#F7D767',
'#FFBD56', '#FFA044', '#EE4F4D')
pal <- colorRampPalette(ramp)
aggrPlotTitle <-
sprintf("2020 UN-Adjusted Population Count \n for %s (1 sq. km resolution)",
if (selectedCountry %in% c("Czech Republic",
"Democratic Republic of Congo",
"Gambia",
"Netherlands")) {
paste("the", selectedCountry)
} else {
selectedCountry
})
terra::plot(x,
col = pal(8)[-1],
axes = TRUE,
buffer = TRUE,
box = TRUE,
cex.main = 1.5,
line.main = 1.25,
main = aggrPlotTitle,
xlab = expression(bold(Longitude)),
ylab = expression(bold(Latitude)),
line.lab = 2.25,
cex.lab = 1.5,
plg = list(title = expression(bold("Persons")),
title.cex = 1.25,
horiz = TRUE,
loc = "bottom",
yjust = 3.5,
x.intersp = 0.6,
inset = c(0, -0.2),
cex = 1.25),
pax = list(cex.axis = 1.7),
mar = c(8.5, 3.5, 4, 2.5))
terra::north(type = 2, xy = "bottomleft", cex = 1)
# if (selectedCountry == "Czech Republic") {
# sbar(100, type="bar", below="km", cex=0.9, xy="bottomright")
# } else if (selectedCountry == "Nigeria") {
# sbar(300, type="bar", below="km", cex=0.9, xy="bottomright")
# }
# x %>%
# plot(col = pal(8)[-1],
# axes = TRUE,
# main = aggrPlotTitle,
# plg = list(legend = c("0-10",
# "10-25",
# "25-50",
# "50-100",
# "100-250",
# "250-1000",
# ">1000"),
# horiz = TRUE,
# x = "bottom",
# title = "Persons per sq. km"))
# title(xlab = expression(bold(Longitude)),
# ylab = expression(bold(Latitude)),
# line = 2,
# cex.lab = 1.20)
#---------------------------------------#
# Source 2: From GADM: Level1Identifier #
#---------------------------------------#
geodata_path(here())
Level1Identifier <- geodata::gadm(country = inputISOLower,
path = geodata_path(),
version = 3.6)
plot(Level1Identifier, add = TRUE)
if(!directOutput) dev.off()
}
createSusceptibleLayer <- function(selectedCountry, rasterAgg = 0) {
#----------------------------------------------------------------#
# Source 1: WorldPop UN-Adjusted Population Count GeoTIFF raster #
#----------------------------------------------------------------#
inputISO <- countrycode(selectedCountry,
origin = 'country.name',
destination = 'iso3c')
inputISOLower <- tolower(inputISO)
url <- paste0("https://data.worldpop.org/GIS/Population/Global_2000_2020_1km_UNadj/2020/", inputISO, "/", inputISOLower, "_ppp_2020_1km_Aggregated_UNadj.tif")
tifFileName <- basename(url) # name of the .tif file
tifFolder <- "tif/" # .tif files should be stored in local tif/ folder
if (!file.exists(paste0(tifFolder, tifFileName)))
{
download.file(url, paste0(tifFolder, tifFileName), mode = "wb")
}
WorldPop <- rast(paste0(tifFolder, tifFileName))
#Gives the five number summary
print(summary(values(WorldPop)))
#Number of cells that have an NA value
print(sum(is.na(values(WorldPop))))
# print(as.raster(WorldPop))
WorldPop <- replace(WorldPop, is.na(WorldPop), 0) # Delete this line for clear plot. Check!!!
# WorldPop <- terra::rast(paste0(tifFolder, tifFileName))
# Use the above line if fully switching over to terra R package completely
# Note: rasterBasePlot.R was developed with the terra::rast()
# Error in (function (classes, fdef, mtable) :
# unable to find an inherited method for function ‘classify’ for signature ‘"RasterLayer"’
# print(WorldPop)
# print(nrow(WorldPop))
# print(ncol(WorldPop))
# print(ncell(WorldPop))
# print(res(WorldPop))
# print(ext(WorldPop))
if (rasterAgg == 0 || rasterAgg == 1) {
Aggregated <- WorldPop
} else {
Aggregated <- aggregate(WorldPop, fact = c(rasterAgg, rasterAgg), fun = sum, na.rm = TRUE)
}
# print(Susceptible)
returnList <- list("Susceptible" = WorldPop, "Aggregated" = Aggregated, "nRows" = nrow(WorldPop), "nCols" = ncol(WorldPop), "nCells" = ncell(WorldPop))
return(returnList)
}Analysis
The following object (a numeric) controls how much the raster is aggregated. Alternative, recommended values include five (5) or fifteen (15).
rasterAgg <- 10Choosing the country to analyze, we can get some codified information about administrative boundaries (borders for the nation and the provinces, states, or regions therein). Look at the output for the selected country; that’s the ISO three letter (3) code.
# Italy, Nigeria, Korea, or any other nation.
selectedCountry <- "Democractic Republic of Congo"
inputISO <- countrycode(selectedCountry,
origin = 'country.name',
destination = 'iso3c')
inputISOLower <- tolower(inputISO)
inputISO## [1] "COD"tifFilename <- sprintf("%s_ppp_2020_1km_Aggregated_UNadj.tif", inputISOLower)
tifPath <- here("tif", tifFilename)
URL <- paste(sep = "/",
"https://data.worldpop.org",
"GIS",
"Population",
"Global_2000_2020_1km_UNadj",
"2020",
inputISO,
tifFilename)
suppressAllButErrors(dir.create(here("tif")))
if (!file.exists(tifPath)) download.file(URL, sprintf("tif/%s", tifFilename), mode = "wb")
WorldPop <- raster(tifPath)
# raster(file.choose())What is a Raster object?
The first output of this block is the whole Raster object’s default
printing format. This is waht you’ll see when you look at a Raster
object, but there are other information available in the object which
are not printed by default. We can view that information by looking at
the object using other functions from the terra package, or by
using built-in R functions like summary.
WorldPop## class : RasterLayer
## dimensions : 2261, 2289, 5175429 (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : 12.20708, 31.28208, -13.45042, 5.39125 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : cod_ppp_2020_1km_Aggregated_UNadj.tif
## names : cod_ppp_2020_1km_Aggregated_UNadj
## values : 0, 83281.66 (min, max)res(WorldPop)## [1] 0.008333333 0.008333333extent(WorldPop)## class : Extent
## xmin : 12.20708
## xmax : 31.28208
## ymin : -13.45042
## ymax : 5.39125origin(WorldPop)## [1] -0.001249999 -0.000416459crs(WorldPop)## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs
## WKT2 2019 representation:
## GEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]],
## CS[ellipsoidal,2],
## AXIS["longitude",east,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## AXIS["latitude",north,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]]]nrow(WorldPop)## [1] 2261ncol(WorldPop)## [1] 2289ncell(WorldPop)## [1] 5175429dim(WorldPop)## [1] 2261 2289 1cellStats(WorldPop == 0, sum) # Number of grid cells that have a population count = 0## [1] 114610cellStats(!is.na(WorldPop), sum) # Number of grid cells that have a population count## [1] 2703157cellStats(is.na(WorldPop), sum) # Number of grid cells that are NA## [1] 2472272cellStats(WorldPop, sum) # Total estimated 2020 population count## [1] 89561404GADM: the first-level administrative boundary identifiers
The first level of administrative boundaries are whatever a country uses to sub-divide itself after the level of nation-state borders. These are things like provinces (as in Canada) or states (as in the United States).
A legacy data acquisition method with GADM
My team and I used
terra::readRDSin the past, but constructing the filename reliably and downloading the data, ensuring it’s locally available, etc. can be an issue, and the function we were using is now deprecated. Thegeodata::gadmfunction is much easier for getting and reading shape files from the GADM project. If the file has already been downloaded to the path you give it, then it is read from that path, and only downloaded from GADM if necessary. We like to use version 3.6 of the data.
geodata_path(here())
Level1Identifier <- geodata::gadm(country = inputISOLower,
path = geodata_path(),
version = 3.6)
## List of all States/Provinces/Regions
print(Level1Identifier$NAME_1)## [1] "Bas-Uélé" "Équateur" "Haut-Katanga" "Haut-Lomami"
## [5] "Haut-Uélé" "Ituri" "Kasaï" "Kasaï-Central"
## [9] "Kasaï-Oriental" "Kinshasa" "Kongo-Central" "Kwango"
## [13] "Kwilu" "Lomami" "Lualaba" "Maï-Ndombe"
## [17] "Maniema" "Mongala" "Nord-Kivu" "Nord-Ubangi"
## [21] "Sankuru" "Sud-Kivu" "Sud-Ubangi" "Tanganyika"
## [25] "Tshopo" "Tshuapa"Plotting a country’s boundaries
plot(Level1Identifier, main = "Level 1 Administrative Boundaries")
Plotting a Raster object
# plot(WorldPop)
#
# plot(WorldPop, col = terrain.colors(255))
#
# par(mfrow = c(1, 2))
# image(log(WorldPop), col = heat.colors(10), main = "heat: 2020 UN-Adjusted Population Count (log-scale) \n (each grid cell is 1 km x 1 km) \n")
# image(log(WorldPop), col = topo.colors(10), main = "topo: 2020 UN-Adjusted Population Count (log-scale) \n (each grid cell is 1 km x 1 km) \n")
#
# plot(log(WorldPop), xlab = "Longitude", ylab = "Latitude", main = "2020 UN-Adjusted Population Count (log-scale) \n (each grid cell is 1 km x 1 km) \n")
#
# createBasePlot(selectedCountry = selectedCountry,
# rasterAgg = 0,
# directOutput = TRUE)
suscLayer <- createSusceptibleLayer(selectedCountry, 0)## cod_ppp_2020_1km_Aggregated_UNadj
## Min. : 0.0
## 1st Qu.: 2.3
## Median : 6.7
## Mean : 33.1
## 3rd Qu.: 20.0
## Max. :83281.7
## NA's :2472272
## [1] 2472272createBasePlot(
selectedCountry = selectedCountry,
susceptible = suscLayer$Susceptible,
directOutput = TRUE # or FALSE
) 
Continue reading the PowerPoint®
When Ashok presents on Thursday, at this point you will return to watching the presentation. After Thursday, when referencing this RMarkdown file, you should switch to the place in the Presentation that continues with a heading like the following heading.
Aggregating a Raster object
WorldPop <- replace(WorldPop, is.na(WorldPop), 0)
WorldPop_aggr <- aggregate(WorldPop,
fact = c(rasterAgg, rasterAgg),
fun = sum,
na.rm = TRUE)
names(WorldPop_aggr) <- "Susceptible"
WorldPop_aggr## class : RasterLayer
## dimensions : 227, 229, 51983 (nrow, ncol, ncell)
## resolution : 0.08333333, 0.08333333 (x, y)
## extent : 12.20708, 31.29042, -13.52542, 5.39125 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : Susceptible
## values : 0, 1911429 (min, max)summary(getValues(WorldPop_aggr))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0 0.0 57.6 1722.9 913.2 1911428.7# Check that the population count before and after aggregation are the same
cellStats(WorldPop, sum) == cellStats(WorldPop_aggr, sum)## [1] TRUEContinue reading the PowerPoint®
When Ashok presents on Thursday, at this point you will return to watching the presentation. After Thursday, when referencing this RMarkdown file, you should switch to the place in the Presentation that continues with a heading like the following heading.
Exporting a Raster object
Ashok should explain why some of the following lines are commented out.
writeRaster(WorldPop_aggr,
# filename = sprintf("%s_unaggregated.nc", inputISO),
filename = sprintf("%s_aggr_0000.tif", inputISO),
# format = "CDF",
format = "GTiff",
varname = "Susceptible",
varunit = "Persons",
longname = "Susceptible",
overwrite = TRUE)Continue reading the PowerPoint®
When Ashok presents on Thursday, at this point you will return to watching the presentation. After Thursday, when referencing this RMarkdown file, you should return to the presentation file Ashok should have shared with you; this is the end of the RMarkdown file.