#' crop_worldpop_tif
#'
#' \code{.zip} files downloaded from worldpop are national-scale. This function
#' crops them to the designated cities, and converts them to much smaller geotif
#' files.
#'
#' @param city Name of city (one of "accra", "kathmandu", or "bristol")
#' @param expand Relative factor by which to expand cropping limits beyond
#' rectangular bounding box of city
#' @return Nothing Files are written to the corresponding directories
#'
#' @note This should be run from somewhere within the \code{whodata} repository
#' structure; anywhere else will likely fail.
#'
#' @export
crop_worldpop_tif <- function (city = "accra", expand = 0.1)
{
bb <- osmdata::getbb (city)
bb_exp <- t (apply (bb, 1, function (i)
i + c (-expand, expand) * diff (i)))
wd <- file.path (here::here (), city, "popdens")
files <- list.files (wd)
files <- file.path (wd, files [grep ("\\.zip", files)])
if (length (files) == 0)
{
message ("No files to convert for ", city)
return (NULL)
}
pb <- txtProgressBar (style = 3)
count <- 1
n <- length (files)
for (f in files)
{
zfiles <- unzip (f, list = TRUE)$Name
f1 <- zfiles [grep ("tif$", zfiles, ignore.case = TRUE)]
unzip (f, files = f1)
r <- raster::raster (f1) %>% raster::crop (bb_exp)
f0 <- paste0 (tools::file_path_sans_ext (f1), ".zip")
file.remove (file.path (wd, f0))
file.remove (f1) # the full geotif dumped in current wd
# faster convert file names to lower case, so
f1 <- paste0 (tools::file_path_sans_ext (f1), ".tif")
f1 <- file.path (wd, f1)
raster::writeRaster (r, f1)
setTxtProgressBar (pb, count / n)
count <- count + 1
}
close (pb)
}
#' crop_global_tif
#'
#' There are no worldpop data for Europe, but the EC Joint Research Centre Data
#' Catalogue offers 250m resolution global population density \code{tif} files.
#' These are huge, and presumed to be downloaded into the root directory of this
#' project. This function crops them to designated cities (actuall only Bristol
#' at present) and deposits the results in the appropriate \code{[city]/popdens}
#' directory.
#'
#' The main JRCDC page is
#' \url{http://data.jrc.ec.europa.eu/dataset/jrc-ghsl-ghs_pop_gpw4_globe_r2015a},
#' leading to the download of GHS_POP_GPW42015_GLOBE_R2015A_54009_250_v1_0.zip,
#' which is 1GB, and is not archived in this repo.
#'
#' @param city Name of city (one of "accra", "kathmandu", or "bristol")
#' @param expand Relative factor by which to expand cropping limits beyond
#' rectangular bounding box of city
#' @return Nothing Files are written to the corresponding directories
#'
#' @note This should be run from somewhere within the \code{whodata} repository
#' structure; anywhere else will likely fail. It also takes a r-e-a-l-l-y long
#' time because of the raster re-projection.
#'
#' @export
crop_global_tif <- function (city = "bristol", expand = 0.1)
{
if (city != "bristol")
stop ("currently hard-coded for bristol only")
bb <- osmdata::getbb (city)
bb_exp <- t (apply (bb, 1, function (i)
i + c (-expand, expand) * diff (i)))
globes <- list.files ()
globes <- globes [grep ("\\.zip", globes)]
zfiles <- unzip (globes, list = TRUE)$Name
f1 <- zfiles [grep ("tif$", zfiles, ignore.case = TRUE)]
unzip (globes, files = f1)
# This file is in World Molleweide (EPSG:54009) projection, for which the
# rough bounds of bristol can be obtained from
# https://epsg.io/map#srs=54009&x=-188375.09900285&y=6060835.47073653&z=7
# These are used to first crop the map in this projection, then re-project
# the cropped version, and then crop again to the nominatim EPSG84 bounds.
bb_moll <- c (-230000.0, -160000.0, 6000000.0, 6060000.0)
bb_moll <- c (-205000.0, -180000.0, 6015000.0, 6040000.0)
message ("cropping global raster to Bristol ... ", appendLF = FALSE)
r <- raster::raster (f1) %>% raster::crop (matrix (bb_moll, nrow = 2))
crs <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
message ("done\nRe-projecting to WGS84 (this will take a long time!) ... ",
appendLF = FALSE)
r <- raster::projectRaster (r, crs = crs) %>% raster::crop (bb_exp)
message ("done")
fnew <- file.path (here::here (), city, "popdens")
if (!file.exists (fnew))
dir.create (fnew)
fnew <- file.path (fnew, basename (f1))
raster::writeRaster (r, fnew)
invisible (file.remove (f1))
invisible (file.remove (dirname (f1)))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.