WorldPop data on Guatemala

Data from the WorldPop project for Guatemala will be used to obtain an estimation of the population in each municipality for 2015.

library(rgdal)
## Loading required package: sp
## rgdal: version: 1.2-16, (SVN revision 701)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.0, released 2017/04/28
##  Path to GDAL shared files: C:/Users/Guillermo/Documents/R/win-library/3.4/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/Guillermo/Documents/R/win-library/3.4/rgdal/proj
##  Linking to sp version: 1.2-5
library(raster)
library(sp)
library(colorRamps)
library(ggmap)
## Loading required package: ggplot2
library(ggplot2)
WorldPopGT <- raster("C:/DATOS/WorldPop/GTM_ppp_v2b_2015/GTM_ppp_v2b_2015.tif")

The municipalities shape file obtained from SEGEPLAN website is most probably outdated. I’ve used a shapefile in geojson format obtained from IGN (national geographic institute of Guatemala) which seems to be more up to date. This file was retrieve in mid 2017 and has been uploaded to the internet archive so it is accesible. The link is mentioned below:

# Old data from SEGEPLAN data. This only contains 334 municipalities. 
gtmMunis = readOGR("../../Guillermo/GIS/municipios_gtm/municipios_GTM.shp", encoding = "utf-8")
## OGR data source with driver: ESRI Shapefile 
## Source: "../../Guillermo/GIS/municipios_gtm/municipios_GTM.shp", layer: "municipios_GTM"
## with 337 features
## It has 12 fields
# This data comes from IGN GIS system. It contains 340 Municipalities. It seems to be less outdated.
# This file can be found in the archive: https://archive.org/download/IGNCartografiaBasicaDivisionPoliticaAdministrativaMunicipios/IGN-cartografia_basica-Division%20politica%20Administrativa%20%28Municipios%29.geojson
gtmMunisIGN = readOGR("../TempData/IGN-cartografia_basica-Division politica Administrativa (Municipios).geojson", encoding="utf-8")
## OGR data source with driver: GeoJSON 
## Source: "../TempData/IGN-cartografia_basica-Division politica Administrativa (Municipios).geojson", layer: "IGN-cartografia_basica-Division politica Administrativa (Municipios)"
## with 344 features
## It has 6 fields

A raw visual overview of the country. Spent a while tweaking the color scale in order to be able to visually recognize very low counts and high counts.

plot(log10(WorldPopGT), colNA="black", col=colorRampPalette(c("#121212", "#123620", "#108650", "#80D6C0", "#DFDF3B"),0.5)(110))

Estimating, for instance, the population in Guatemala department. Result is 3.4 million people.

print(sum(extract(WorldPopGT, gtmMunis[gtmMunis$Cod_Dep==1,], fun = sum, na.rm = TRUE ) ))
## Warning in .local(x, y, ...): Transforming SpatialPolygons to the CRS of
## the Raster
## [1] 3439633
print(sum(extract(WorldPopGT, gtmMunisIGN[gtmMunisIGN$COD_DEPT__ %in% c("01"),], fun = sum, na.rm = TRUE ) ))
## [1] 3436954

What about the population of the entire country? 15.9 million. National statistics institute of Guatemala projections estimate this number to be 16176034 for 2015. I would believe more in this data since it is based on other covariates, while INE projections are based on a raw population growth rate, which is probably estimated from official data on live births and deaths.

tempdata = as.array(WorldPopGT)
sum(tempdata, na.rm = T)
## [1] 15929699

Mixing googlemaps and the worldpop raster

Now I leave you with a nice snippet to explore this data with google maps images. Just change the location and zoom values:

gmap = get_map(location=c(-90.50,14.60), zoom=12, scale=1, maptype = "roadmap", source="google", color="bw")
gmapbbox = attr(gmap, "bb")
subpop <- crop(WorldPopGT, extent(gmapbbox$ll.lon, gmapbbox$ur.lon, gmapbbox$ll.lat, gmapbbox$ur.lat))
subpopDF = data.frame(rasterToPoints(subpop))
colnames(subpopDF) = c("lon", "lat", "z")
bm = ggmap(gmap)
bm = bm + geom_raster(data = subpopDF, aes(y=lat, x=lon, fill=z), alpha=0.6) + 
    scale_fill_gradientn(colours=rev(c("#FF1202FF", "#FDAD22FF", "#FDFD22FF", "#00000077")), values = c(0,0.1,0.5,1) ) + coord_cartesian()
plot(bm)

Now let’s estimate the population by municipality, which is what brought me here in the first place:

Putting the data in the municipalities database and exporting to a CSV: