2  Handling of spatial data in R

2.1 Introduction

Geostatistical analysis primarily deals with point-referenced data. However, geostatistical modelling often requires more than point data. Raster and areal (polygon) data are frequently used to build essential covariates, enhance spatial understanding and provide environmental context. For instance, population density, climate, land use and elevation are often key covariates in spatial analyses and usually come in raster or polygon formats. Efficiently handling these different data types is critical for successful model building. This chapter provides a comprehensive guide to managing spatial data in R using the sf and terra packages to support a complete geostatistical workflow.

2.2 Spatial data handling in geostatistical analysis

Throughout a geostatistical analysis, handling spatial data takes place at multiple key stages:

  1. Retrieving external spatial data sources. Before starting the modelling process, one often needs to acquire spatial datasets from external sources to improve model accuracy. These could include satellite-derived environmental variables, population data from WorldPop or climate data from WorldClim. Acquiring these datasets and ensuring they have compatible formats and resolutions is an essential early step in spatial analysis.

  2. Importing and standardising spatial data. Once external data are obtained, they must be imported into R. This involves handling various file formats such as shapefiles or GeoPackages for vector data and GeoTIFFs for raster data. Different datasets might use different coordinate reference systems (CRSs), requiring CRS standardisation to ensure that all datasets align correctly. Failure to standardise CRSs can result in spatial misalignment and incorrect results.

  3. Extracting covariate information for modelling. Covariate extraction is an essential step in geostatistical modelling. After importing spatial data, it is necessary to extract relevant covariate information both for the sampled locations, where we have observed data, and the prediction locations, where we wish to make predictions. This step involves linking raster or polygonal covariates, such as climate, population or land-cover data, to geostatistical data points.

  4. Prediction and creation of a spatial grid. For predictive geostatistical models, a regular grid of points is often created over the study region. Covariates must then be assigned to each grid point to make predictions across the entire region of interest. Creating predictive grids and linking them to the necessary covariates is key to generating continuous spatial predictions from geostatistical models.

  5. Visualising spatial data. Visualisation is crucial for exploring spatial data, interpreting model results and communicating findings. Whether working with point-referenced data, polygons or rasters, clear and effective visualisation helps reveal patterns that inform the modelling process. It can also highlight covariate trends, spatial clusters and uncertainty in predictions.

2.3 Accessing covariates for disease mapping

Covariates can play a crucial role in understanding and modelling spatial variation in disease risk. In geostatistical analysis, incorporating environmental, demographic and climatic variables might improve the predictive power of models or reduce the uncertainty in predictions. These covariates can influence the spread of infectious diseases, the distribution of disease vectors and the socio-economic conditions that affect health outcomes.

A wide range of open-access spatial datasets provide covariates for disease mapping. These include population density, climate, land cover and human infrastructure data, which often come in raster or polygon formats. Administrative boundaries are also crucial for public health analysis because they help organise and aggregate data at different levels, such as countries, regions or districts. For example, aggregating health outcomes or covariate data by administrative units allows researchers to identify geographic disparities and allocate resources accordingly. Open-source platforms such as geoBoundaries provide administrative boundary data for many countries at various levels of granularity. These data are available in formats compatible with geospatial analysis tools in R and can be integrated into geostatistical workflows.

Table 2.1 below summarises key sources of covariates useful in disease mapping. Each dataset offers specific types of data, from satellite-derived environmental variables to gridded population estimates and administrative boundaries, which can be accessed through various R packages or APIs.

Table 2.1: Some key sources of covariates useful in disease mapping and R packages to retrieve them.
Source Data Type Description R Package
WorldClim Climate (temperature, rainfall) Global climate data geodata
MODIS Remote Sensing Satellite imagery (e.g., vegetation indices) MODIStsp
OpenStreetMap (OSM) Human settlements, roads Global geographic features osmdata
Google Earth Engine Satellite imagery Large-scale environmental data analysis rgee
WorldPop Population data Gridded population counts and densities RiskMap

2.3.1 Example: Downloading administrative boundaries

Administrative boundaries provide an essential spatial structure for many types of geostatistical analyses, particularly in disease mapping where data is often aggregated by administrative units such as regions, districts or provinces. The geoBoundaries datasets, which are accessible via the geobounds package, provide openly available administrative boundary data for nearly every country, allowing researchers to integrate these boundaries into their analyses seamlessly.

Boundary datasets are distributed under different licences, which specify how they may be reused, and attribution identifies the organisation or individuals who produced them. These conditions are especially important when publishing maps. Depending on the licence, the source and licence may need to be acknowledged in the figure caption, acknowledgements or data-availability statement. The gb_get_metadata() function can be used to inspect this information before a map is shared or published. By default, geobounds downloads boundaries from the open gbOpen release.

library(geobounds)

# Download the country and first-level administrative boundaries
liberia_admin0 <- gb_get_adm0("Liberia")
liberia_admin1 <- gb_get_adm1("Liberia")
# Plot the administrative boundaries
plot(liberia_admin0$geometry)
plot(liberia_admin1$geometry)
(a) Admin level 0
(b) Admin level 1
Figure 2.1: Administrative boundaries for Liberia retrieved using the geobounds package.

2.3.2 Example: Downloading population data

Population distribution is often relevant to geostatistical models for public health research. In this section, we demonstrate how to retrieve high-resolution population data for Liberia from WorldPop. WorldPop estimates the number of people in each cell of a fine grid by distributing census totals with statistical models that use covariates such as settlement maps, roads and night-time lights. Its products differ in ways that matter for analysis:

  • Unconstrained estimates place people in every land cell, so sparsely populated or uninhabited areas receive small counts. They are available for every year from 2000 to 2020, for the whole population and by age and sex, at about 100 m and 1 km.
  • Constrained estimates place people only in cells mapped as built-up settlements. They locate people more sharply, particularly at fine resolution, but can miss small rural settlements. They are available from 2015 onwards.
  • UN-adjusted estimates rescale national totals to match United Nations population estimates.
  • Bottom-up estimates, available for a few countries, are modelled from household surveys rather than from census totals, and come with measures of uncertainty.

Unconstrained estimates suit population weights on grids of a few kilometres and the analysis of surveys carried out before 2015. Constrained estimates suit fine-scale analyses of recent years, and UN-adjusted totals suit reporting absolute numbers of people. Whatever the product, check whether it contains counts per cell or population densities, because this determines how its values are aggregated.

The RiskMap function get_worldpop() finds a product in the WorldPop catalogue, downloads it and returns it as a raster of people per cell. Here we download the unconstrained estimates for Liberia in 2014 at about 100 m, saving the GeoTIFF file in the data folder.

# Load the RiskMap package
library(RiskMap)

# Download WorldPop population counts for Liberia in 2014
lbr_pop_100 <- get_worldpop("LBR", year = 2014, dest_dir = "data")

The ages and sex arguments select age groups, whose counts are added together. For example, get_worldpop("LBR", year = 2011, ages = c(20, Inf), resolution = "1km") returns the number of adults aged 20 and over in each cell of about 1 km, which we use in Chapter 4. The constrained and un_adjusted arguments select the other products described above.

2.3.3 Example: Downloading data using Google Earth Engine

Google Earth Engine (GEE) is a powerful platform that hosts petabytes of environmental and geospatial data, including climate, topography, and land cover datasets. The rgee R package provides a convenient interface to GEE from within R.

🔗 Getting started with rgee: To use rgee, you need a Google account and must complete a one-time setup process. Follow the official guide here: https://r-spatial.github.io/rgee/. For further tutorials and use cases, see: https://csaybar.github.io/rgee-examples/.

In this example, we retrieve elevation data for Liberia from the SRTM 90 m dataset. This is the native source product. We request an output scale of 1,000 m below to reduce the size of the example, so the downloaded raster used in the book has an approximately 1 km grid rather than the native 90 m grid. The steps include initialising Earth Engine, accessing the desired data product with the ee$Image("ASSET_ID") function (the ASSET_ID can be found on the dataset page under Earth Engine Snippet), converting the boundary of Liberia into a format readable by GEE, clipping the global elevation raster to that boundary and downloading the result as a raster object in R.

Note

This example assumes you have already completed rgee setup and authenticated using ee_Initialize().

# Load required packages
library(rgee)
library(terra)

# Initialize Earth Engine (first time will prompt browser login)
ee_Initialize(
  project = "YOUR-GOOGLE-CLOUD-PROJECT-ID",
  drive = TRUE,
  quiet = TRUE
)

# Convert Liberia boundary to Earth Engine Geometry
liberia_ee <- sf_as_ee(liberia_admin0)

# Access the elevation band of the global SRTM dataset
elev <- ee$Image("CGIAR/SRTM90_V4")$select("elevation")

# Inspect general information about the data product
ee_print(elev)

# Clip the elevation data to Liberia
elev_liberia <- elev$clip(liberia_ee)

# Download clipped elevation data as a raster using Google Drive
elev_rast <- ee_as_rast(
  image = elev_liberia,
  region = liberia_ee$geometry(),
  via = "drive",
  scale = 1000, # Request an output grid of approximately 1 km
  quiet = TRUE
)

# Mask cells whose centres fall outside Liberia
elev_rast <- mask(elev_rast, liberia_admin0, touches = FALSE)

# Plot the result
plot(elev_rast, main = "Elevation in Liberia (m)")
plot(liberia_admin0$geometry, add = TRUE)
Figure 2.2: Elevation in metres on the approximately 1 km grid requested from the SRTM 90 m source product.

2.4 Importing and standardising spatial data

In geostatistical analysis, importing and standardising spatial data is critical to ensuring that data from different sources align and can be used effectively. Vector data, such as points, lines and polygons, and raster data can come in various formats and may use different coordinate reference systems (CRSs). Accurate spatial analysis requires the data to be imported correctly and to use consistent projections and formats. This section explains how to import vector and raster data into R, introduces CRSs and shows how to ensure that different datasets align for subsequent geostatistical analysis.

2.4.1 Importing vector data

Vector data are typically stored in formats such as shapefiles (.shp) or geopackages (.gpkg). The sf (simple features) package in R is the most common tool for handling vector data.

# Load the sf package
library(sf)

# Import a shapefile (e.g., administrative boundaries)
admin <- st_read("path_to_your_shapefile/Admin_Boundaries.shp")

# Inspect the data
print(admin)
plot(admin)  # Basic plot of the shapefile

The st_read() function reads various spatial data formats and automatically recognises file types.

2.4.2 Importing raster data

Raster data consists of a grid of cells, where each cell holds a value representing a spatial attribute such as elevation or temperature. The terra package in R is designed to work with raster data and has superseded the older raster package due to better performance and greater functionality. Let’s see now how to import a GeoTIFF file using terra. The function get_worldpop() in Section 2.3.2 returned the population raster directly and also saved it as a GeoTIFF file, which rast() can read again at any time.

# Import a raster file
lbr_pop_100 <- rast("data/lbr_ppp_2014.tif")
# Load the terra package
library(terra)

# Inspect the raster data
print(lbr_pop_100)
class       : SpatRaster
size        : 5040, 4945, 1  (nrow, ncol, nlyr)
resolution  : 0.0008333333, 0.0008333333  (x, y)
extent      : -11.48625, -7.365417, 4.352084, 8.552084  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source      : lbr_pop_100.tif
name        : population
min value   :   0.006427
max value   :  92.716873
# Basic plot of the raster (log-scale)
plot(log(lbr_pop_100), main = "Population per grid cell (log scale)")
plot(liberia_admin0$geometry, add = TRUE)
Figure 2.3: Estimated population per approximately 100 m grid cell in Liberia in 2014, shown on the log scale.

The output of print() provides a detailed summary of the raster’s properties. The rast() function reads raster files in formats such as GeoTIFF and ASCII grid. We then use plot() to visualise the raster quickly.

2.4.3 Understanding coordinate reference systems (CRS)

When working with spatial data, especially from different sources, one of the most critical tasks is ensuring that all datasets share the same coordinate reference system (CRS). A CRS defines how two-dimensional map data correspond to locations on the three-dimensional Earth. If different layers, such as raster and vector data, have different CRSs, they may not align correctly when plotted or analysed together, leading to inaccurate analyses or visualisations.

CRSs fall into two broad groups.

  • Geographic CRSs use angular coordinates, usually longitude and latitude, to represent locations on the Earth’s surface. A common example is WGS 84 (EPSG:4326), which is widely used by GPS and global datasets.

  • Projected CRSs represent locations on a flat surface in linear units such as metres. Every projection introduces distortion, so a CRS should be chosen according to the study area and the property that matters most, such as distance, area or direction. Examples include Universal Transverse Mercator (UTM) and Albers Equal Area projections.

Longitude and latitude are angular coordinates, but this does not mean that every spatial function returns distances in degrees. CRS-aware functions can account for the Earth’s curvature. For example, sf::st_distance() calculates great-circle distances for longitude-latitude geometries by default and reports them in linear units. Ordinary numerical calculations do not have this spatial context. If longitude and latitude are extracted as columns or as a coordinate matrix and Euclidean distance is then applied, the result treats degrees as if they were Cartesian units.

This distinction is important for geostatistical modelling. RiskMap currently extracts the coordinates from an sf object and uses Euclidean calculations internally, so the modelling coordinates should use a suitable projected CRS. A local study contained within one UTM zone can often use that zone. Studies spanning several zones or large regions may require a regional projection chosen for the distances or areas being analysed. For global analyses, methods that work directly with spherical or ellipsoidal geometry may be more appropriate.

2.4.4 EPSG codes

An EPSG code is a unique identifier for a coordinate reference system or another object in the EPSG Geodetic Parameter Dataset. The dataset is maintained by the Geodesy Subcommittee of the Geomatics Committee of the International Association of Oil & Gas Producers. The familiar EPSG name reflects its origin in the former European Petroleum Survey Group. These identifiers make it easier for software and users to refer to the same definitions consistently.

Some commonly used EPSG codes are:

  • EPSG:4326: This code represents WGS84, the most commonly used geographic CRS, which uses latitude and longitude to describe locations on the Earth’s surface. It is the default CRS for global datasets and GPS systems.

  • EPSG:326XX identifies WGS 84 / UTM zones in the northern hemisphere. Liberia lies within UTM zone 29N, represented by EPSG:32629. UTM is designed for local and regional work within a zone, but it is not automatically appropriate for every study area.

  • EPSG:3857 represents the Web Mercator projection, which is widely used for web mapping services, including Google Maps, OpenStreetMap and Bing Maps. This projected CRS uses metres as the unit of distance and is optimised for displaying maps on a two-dimensional plane, although it distorts area and distance, especially at high latitudes. It is well suited to interactive online mapping but not to precise distance-based geostatistical analyses.

For a more in-depth discussion on how to choose an appropriate CRS, we recommend referring to the excellent treatment in Geocomputation with R (Lovelace et al. 2020), especially Section 7.6: Which CRS?.

2.4.5 Convert a data frame to an sf object

In geospatial analysis, data are often provided in tabular formats such as CSV files that contain spatial coordinates, for example latitude and longitude. To use these data effectively in R, the data frame must be converted into an sf object, which is the standard format for working with spatial vector data in R. We demonstrate this conversion using the liberia data frame from the RiskMap package.

# Load the Liberia data set from RiskMap
data("liberia")

# Convert the data frame to an sf object
liberia_sf <- st_as_sf(liberia,
                       coords = c("long", "lat"),
                       crs = 4326)

# Inspect the new sf object
liberia_sf
plot(liberia_sf)

The st_as_sf() function converts the data frame into an sf object. The coords argument identifies the columns containing the coordinates. The crs argument assigns WGS 84 (EPSG:4326), which is the CRS already used by these longitude and latitude values. Assigning a CRS describes existing coordinates but does not change them. By contrast, st_transform() calculates new coordinates in another CRS and should only be used after the source CRS has been assigned correctly.

The resulting sf object can be used for spatial joins, distance calculations and mapping. The coordinate columns have been replaced by a geometry column. Set remove = FALSE in st_as_sf() if the original columns should also be retained.

2.4.6 Working with CRSs in R

When working with data from multiple sources, such as environmental layers, population data, or administrative boundaries, ensuring that all datasets share the same CRS is essential for accurate spatial analysis. This section covers the core tasks involved in managing CRSs in R: checking the CRS of spatial data to ensure datasets are compatible and reprojecting spatial data into a common CRS when necessary.

2.4.6.1 Checking the CRS of spatial data

Before performing any spatial operation, it is crucial to check the CRS of your spatial datasets. Knowing whether your data uses geographic coordinates (e.g., WGS84) or a projected coordinate system (e.g., UTM) helps ensure that they are aligned and ready for analysis. Both sf and terra provide functions to retrieve and inspect the CRS, ensuring datasets are spatially aligned before analysis. The st_crs() function retrieves the CRS information for vector data.

# Check the CRS of the vector data
st_crs(liberia_sf)
Coordinate Reference System:
  User input: EPSG:4326 
  wkt:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

The same can be achieved for raster data with the crs() function from the terra package.

# Check the CRS of the raster data
crs(lbr_pop_100, proj = TRUE, describe = TRUE)

Two CRS descriptions can use different text while referring to equivalent coordinate systems. Compare complete CRS objects rather than comparing an EPSG number or a stored text field.

st_crs(liberia_sf) == st_crs(lbr_pop_100)
[1] TRUE

A matching CRS is necessary for overlaying these objects, although it does not by itself guarantee that their extents or resolutions are suitable for the intended operation.

2.4.6.2 Reprojecting spatial data to a common CRS

If datasets have different CRSs, one or more of them may need to be transformed before they can be overlaid or analysed together. The st_transform() function transforms vector coordinates to a specified CRS. RiskMap’s propose_utm() function identifies the UTM zone containing the centre of an object. This is a useful candidate for a compact study area, but the analyst must still check that the study lies suitably within that zone.

# Obtain the UTM CRS proposed for Liberia
liberia_utm_crs <- st_crs(propose_utm(liberia_sf))
liberia_utm_crs
Coordinate Reference System:
  User input: EPSG:32629 
  wkt:
PROJCRS["WGS 84 / UTM zone 29N",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 29N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-9,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Navigation and medium accuracy spatial referencing."],
        AREA["Between 12°W and 6°W, northern hemisphere between equator and 84°N, onshore and offshore. Algeria. Côte D'Ivoire (Ivory Coast). Faroe Islands. Guinea. Ireland. Jan Mayen. Liberia, Mali. Mauritania. Morocco. Portugal. Sierra Leone. Spain. United Kingdom (UK). Western Sahara."],
        BBOX[0,-12.01,84.01,-6]],
    ID["EPSG",32629]]
# Reproject the vector data
liberia_sf_utm <- st_transform(liberia_sf,
                               crs = liberia_utm_crs)

Reprojecting raster data is more complex than reprojecting vector data because it creates a new grid. The new grid may have a different origin, cell size, dimensions and alignment, so values must be estimated from the cells of the original raster. The method must match what the values represent. Bilinear or cubic interpolation is suitable for continuous fields such as elevation or temperature. Nearest-neighbour resampling is appropriate for categorical rasters because it does not create new category codes. Counts and totals are extensive quantities whose values depend on the area represented and are added when areas are combined. They therefore require a total-preserving method rather than interpolation. By contrast, continuous fields and densities are intensive quantities whose values do not simply increase when the represented area becomes larger.

The project() function from the terra package can reproject a raster. The WorldPop raster contains counts, so the "sum" method uses weighted contributions from the input cells to preserve the total population. When several rasters must align, it is preferable to supply a template raster that defines the desired CRS, extent, resolution and origin. Here we supply only a CRS because this example does not need to align the output with another raster.

lbr_pop_100_utm <- project(lbr_pop_100,
                           crs(liberia_sf_utm),
                           method = "sum")

The numerical resolution before and after reprojection cannot be compared directly when the units change from degrees to metres. More importantly, resolution describes the output grid rather than the effective detail of the source data. Reprojecting a coarse raster onto a finer grid does not create new spatial information. Projection distortion and resampling may also introduce edge effects or some loss of information.

As a general rule, avoid reprojecting a raster when transforming the vector layer to the raster’s CRS achieves the same valid analysis. This avoids unnecessary resampling and alteration of the raster values. Raster reprojection is still appropriate when an analysis requires a particular projected grid or when several rasters must share the same geometry. In these cases, choose the output grid deliberately and use a resampling method suited to the variable. For a geostatistical workflow, covariates can be extracted in the raster’s native CRS before the point data are transformed to a suitable projected CRS for model fitting.

2.5 Extracting covariate data

In geostatistical models, the inclusion of relevant covariates (environmental, demographic, or climatic) can potentially enhance predictive accuracy. Covariate data often come from raster or polygon sources, and extracting these values for point locations is essential to link spatial context to point-referenced data. These covariates could include variables like temperature, elevation, land cover, or population density, which influence the spatial distribution of diseases. In this section, we will cover how to extract covariates at point locations from both polygon and raster layers.

2.5.0.1 Extracting covariates from polygon layers

Polygon layers represent discrete spatial features such as administrative boundaries, land cover categories, or protected areas. These layers typically include associated attributes (e.g. region names, population totals, land type), which can serve as covariates in geostatistical models. A common task in spatial data analysis is to assign attributes from a polygon layer to a set of point-referenced data, for example, determining the administrative region or land use category that each survey location falls within. This is done using a spatial join, which overlays the point and polygon layers and attaches the attributes of the enclosing polygon to each point. Below is an example where we assign administrative region names (admin1) to a set of point locations in Liberia:

# Perform a spatial join to transfer polygon attributes to the points
points_with_admin <- st_join(liberia_sf, liberia_admin1["shapeName"])

# The points now include the attributes of the polygon layer
head(points_with_admin)

This spatial join enriches each point with the corresponding polygon attribute, in this case, the admin1 name (shapeName). These attributes can then be used as categorical covariates in geostatistical models to account for regional variation in disease risk or programmatic coverage.

2.5.1 Extracting covariates from raster layers

Raster layers can represent continuous fields, such as elevation or temperature, categorical variables, such as land-cover class, or quantities associated with cells, such as population count. Covariate values can be extracted for specific points using the extract() function from the terra package. By default, each point receives the value of the raster cell that contains it. This is a point extraction, so its scientific interpretation is the value represented by that cell rather than an average over the surrounding area.

# Extract raster values at the point locations
covariate_values <- extract(lbr_pop_100, liberia_sf)

# Combine the extracted values with the point data
liberia_sf$pop_cell <- covariate_values$lbr_ppp_2014

# View the updated dataset
head(liberia_sf)

In this example, extract() assigns the estimated population count in the approximately 100 m cell containing each survey location. This differs from population density, which expresses population per unit area. A cell count also depends on the raster resolution because changing the cell size changes the count associated with a point even if the underlying population distribution is unchanged.

2.5.1.1 Missing raster values

Raster extraction can return NA for several reasons. A sampling location may fall outside the raster extent, the source product may contain a missing cell, or preprocessing may have masked the location. Before replacing a missing value, check the raster extent and CRS, inspect the source product and plot the location with the raster. A missing value outside the product’s coverage or across a coastline should generally remain missing; depending on the analysis, the corresponding observation may need to be excluded or a different covariate source used.

For an isolated gap in a smoothly varying continuous field, such as elevation, a nearby valid cell may sometimes be a reasonable substitute. The search_radius argument of extract() searches for the nearest non-missing raster cell within a specified distance. For longitude-latitude rasters the radius is given in metres. The following controlled example makes the elevation cells under two survey locations missing and then searches within 2 km:

elevation_with_gap <- elev_rast
survey_cells <- cellFromXY(elevation_with_gap,
                           st_coordinates(liberia_sf[1:2, ]))
elevation_with_gap[survey_cells] <- NA

nearest_elevation <- extract(elevation_with_gap,
                             liberia_sf[1:2, ],
                             search_radius = 2000)
nearest_elevation

When a replacement is found, the result also reports its distance and cell number. Set the radius in advance from subject-matter knowledge rather than increasing it until a value is found. This approach is unsuitable for broad areas of missing data, locations outside source coverage or situations where a barrier makes the nearest cell scientifically unrelated. For categorical rasters, do not average category codes, retain the missing value unless a defensible rule supports using a nearby category or a neighbourhood mode.

Finding a nearby cell is an imputation that attempts to recover the value of the covariate at the point. A buffer summary answers a different question by defining the covariate as exposure over a neighbourhood. The summary must match both the scientific question and the type of raster quantity. A mean or area-weighted mean can represent average elevation, temperature or population density around a location. A sum is appropriate when the aim is to estimate an extensive quantity, such as the total population within the neighbourhood. Category codes should not be averaged. Proportions, a majority class or a complete frequency table may instead be meaningful.

The following example creates a 5 km buffer around each survey location and calculates the buffer totals for the first two locations. Because the WorldPop layer contains population counts rather than densities, we sum the cell values. Cells that cross the buffer boundary should contribute in proportion to the area that lies inside it. The argument weights = TRUE estimates these overlap fractions approximately and is generally faster, whereas exact = TRUE calculates them more accurately but can be substantially slower. We use exact = TRUE here to demonstrate the more precise calculation. The distinction matters most when cells are large relative to the buffer or the boundary intersects many cells. For a count raster, either method assumes that the quantity is distributed uniformly within each intersected cell when allocating a fraction of its value. The same call can be applied to all the buffers when these totals are required for an analysis.

# Create buffers around each point (e.g., 5 km radius)
buffered_points <- st_buffer(liberia_sf_utm, dist = 5000)

# Plot the buffers for visualisation
plot(st_geometry(buffered_points), border = "black")
plot(st_geometry(liberia_sf_utm), add = TRUE, pch = 19, col = "red",
     cex = .2)
Figure 2.4: Survey locations (red dots) surrounded by a 5 km buffer (dark circles). 
# Transform the buffers back to the raster CRS
buffered_points_wgs84 <- st_transform(buffered_points,
                                      crs = crs(lbr_pop_100))

# Sum population counts, accounting for partial cell overlap
population_5km <- extract(lbr_pop_100,
                          buffered_points_wgs84[1:2, ],
                          fun = sum,
                          exact = TRUE,
                          na.rm = TRUE)
population_5km

The choice between point and buffer extraction is therefore part of defining the covariate, not merely a computational preference. A buffer may be useful when disease risk is influenced by the surrounding environment rather than conditions in one raster cell. Its radius should be chosen from subject-matter knowledge and, where appropriate, examined in a sensitivity analysis.

Missing raster cells and study boundaries require particular care. With na.rm = TRUE, a mean describes only the non-missing part of a buffer, while a sum can underestimate a total if populated cells are missing. Buffers near a coastline or the edge of source coverage may contain much less observed land area than inland buffers. Before using the result, inspect raster coverage within each buffer and decide whether the intended exposure concerns the full circular neighbourhood, land only or another scientifically defined area.

2.6 Creating a predictive grid

A predictive grid is a regularly arranged set of prediction locations covering the study region. It may be represented by cell polygons or by points at the cell centres, depending on what the predictions represent and how they will be used. Raster pixels are cells with associated values, whereas a grid-centre point has a location but no area of its own.

Here, st_make_grid(what = "centers") creates points separated by 5 km rather than square grid cells. Each point represents a location at which the model will make a point prediction. Reducing the spacing produces more prediction locations and a finer-looking map, but it also increases the computational cost. It does not create information at scales finer than those supported by the data and model.

# Transform the Liberia boundary to the modelling CRS
liberia_admin0_utm <- liberia_admin0 |>
  st_transform(crs = liberia_utm_crs)

# Generate grid-centre points at 5 km spacing
pred_locations <- st_make_grid(liberia_admin0_utm,
                               cellsize = 5000,
                               what = "centers")

# Exclude locations that fall outside the study area
pred_locations <- st_intersection(pred_locations, liberia_admin0_utm)

# Visualise the result
plot(liberia_admin0_utm$geometry, col = "white")
plot(pred_locations, cex = .01, col = "red", pch = 19, add = TRUE)
Figure 2.5: Prediction locations at the centres of a regular grid with 5 km spacing. Points are generated in UTM coordinates and filtered to remain within Liberia.

The RiskMap create_grid() function automates these steps and will be used for the remainder of the book. It returns regularly spaced sf point locations within the supplied region. The spacing argument is interpreted in kilometres by default, while distance_units = "m" can be used to specify it in metres. If the region uses longitude and latitude, create_grid() transforms it to a suitable UTM coordinate system and reports the transformation. Transforming the region beforehand, as above, gives the user explicit control over the projected coordinate system.

The equivalent 5 km prediction grid can therefore be created as follows:

pred_locations_riskmap <- create_grid(liberia_admin0_utm, spacing = 5)

2.6.1 Building the grid from population counts

When predictions will be combined using population, for example to average prevalence over the people living in an area or to count the people affected, the grid can be built from the population raster itself. We aggregate the counts to the resolution of the grid and use the centres of the aggregated cells as prediction locations. The population of each cell is then attached to the grid by construction, with no extraction step.

The counts must be summed, not averaged or interpolated, because they are an extensive quantity (see Section 2.4.6). We first project the raster to the modelling CRS at 100 m with the total-preserving "sum" method, choosing a resolution that divides the 5 km spacing exactly. We then add up blocks of 50 by 50 cells with aggregate().

# Project the population counts to UTM at 100 m, preserving totals
lbr_pop_utm <- project(lbr_pop_100, crs(liberia_admin0_utm),
                       res = 100, method = "sum")

# Sum blocks of 50 by 50 cells into 5 km cells
lbr_pop_5km <- aggregate(lbr_pop_utm, fact = 50, fun = "sum",
                         na.rm = TRUE)

# Use the 5 km cells that overlap Liberia as prediction locations
pop_grid <- create_grid(liberia_admin0_utm, raster = lbr_pop_5km,
                        keep = "overlap")

# Compare the total population before and after aggregation
c(original = global(lbr_pop_100, "sum", na.rm = TRUE)[[1]],
  grid = sum(pop_grid$population, na.rm = TRUE))
original     grid 
 3977215  3976799

The grid holds almost the whole population; the small difference comes from WorldPop cells just outside the boundary used here. The grid stores the population of each cell in the column population, so it can be used directly as weights, as we do in Chapter 4. We keep every cell that overlaps the country with keep = "overlap". Keeping only cells whose centre lies inside the boundary would leave out the people living in the edge cells, which along the coast of Liberia is a large part of the population.

2.7 Visualising spatial data

Visualisation is a key part of spatial data analysis because it allows spatial patterns and relationships to be explored and communicated effectively. So far, we have used basic plotting functions for quick inspection. We now introduce the ggplot2 package, which provides a flexible framework for combining geographic data in simple or complex figures. Functions such as geom_sf() and geom_raster() can display point locations, polygons and continuous raster data.

2.7.1 Visualising point data

Point data often represents the locations of observations (e.g., disease cases, sampling sites). ggplot2 allows you to plot these points and optionally colour them by a covariate (e.g., disease prevalence or population). Here is an example that uses the Liberia data.

# Load necessary libraries
library(ggplot2)

# Create a new variable with prevalence in the dataset
liberia_sf$prevalence <- liberia$npos / liberia$ntest

# Plot only the locations
ggplot(data = liberia_sf) +
  geom_sf(col = "black") +
  theme_minimal() +
  labs(title = "Survey locations")

# Colour the points according to prevalence
ggplot(data = liberia_sf) +
  geom_sf(aes(color = prevalence)) +
  scale_color_viridis_c(labels = scales::label_percent()) +
  theme_minimal() +
  labs(title = "Onchocerciasis in Liberia",
       color = "Prevalence (%)")
(a) Survey locations.
(b) Survey locations coloured by prevalence.
Figure 2.6: Using geom_sf() to plot point data.

Here, geom_sf() plots the spatial points. The expression aes(color = prevalence) specifies that the points should be coloured according to the prevalence covariate, providing a visual representation of spatial variation in disease risk. The scale_color_viridis_c() function applies the Viridis colour scale, which is suitable for continuous data and accessible to readers with colour-vision deficiencies. The labels = scales::label_percent() argument displays the colour-scale labels as percentages rather than raw decimal values. The theme_minimal() function removes unnecessary background elements, while labs() defines the plot title and colour-legend label.

2.7.2 Visualising polygon data

Polygon data typically represents administrative boundaries, land use, or other regional divisions. We can still use geom_sf() to create maps of polygons, optionally filling them by a covariate.

# Plot Liberia admin 1 level boundaries
ggplot(data = liberia_admin1) +
  geom_sf() +
  theme_minimal() +
  labs()

# We compute the area of each polygon
liberia_admin1$area <- as.numeric(st_area(liberia_admin1) / 1000 ^ 2)

# Colour the polygons according to this new variable
ggplot(data = liberia_admin1) +
  geom_sf(aes(fill = area), color = "black") +
  scale_fill_distiller(direction = -1) +
  theme_minimal() +
  labs(fill = "Area km^2")
(a) Admin 1 units for Liberia.
(b) Admin 1 units coloured according to their area.
Figure 2.7: Using geom_sf() to plot polygon data.

In this code, aes(fill = area) fills each polygon with a colour corresponding to its area. The color = "black" argument outlines the polygons in black. Setting fill = NA would make the polygons transparent while retaining their borders. The scale_fill_distiller(direction = -1) function applies a colour gradient from ColorBrewer, with direction = -1 reversing the gradient so that larger areas have darker colours.

2.7.3 Visualising raster data

Raster data can be visualised with ggplot2 by using the terra package’s as.data.frame() function to convert the raster into a data frame of coordinates and values.

# Convert the raster to a data frame for ggplot2
raster_df <- as.data.frame(log(lbr_pop_100), xy = TRUE)

# Rename the raster value column used for the fill scale
names(raster_df)[names(raster_df) == "population"] <- "log_population"

# Plot raster using ggplot2
ggplot(data = raster_df) +
  geom_raster(aes(x = x, y = y, fill = log_population)) +
  scale_fill_viridis_c() +
  coord_sf(crs = 4326) +
  theme_minimal() +
  labs(title = "Population per grid cell (log scale)",
       fill = "Log population",
       x = "",
       y = "")
Figure 2.8: Estimated population per approximately 100 m grid cell in Liberia in 2014, shown on the log scale using geom_raster().

In this example as.data.frame() converts the raster into a data frame with x and y coordinates and their corresponding raster values. The geom_raster() function then colours the cells according to their log population count.

2.7.4 Combining multiple spatial data types

It is often useful to combine different spatial data types, such as points, polygons and rasters, in a single visualisation. The ggplot2 package can overlay these layers to provide a more comprehensive view of the spatial data.

# Combine points, polygons, and raster data in one plot
combined_map <- ggplot() +
  geom_raster(data = raster_df,
              aes(x = x, y = y, fill = log_population)) +
  geom_sf(data = liberia_admin1, fill = NA, color = "black") +
  geom_sf(data = liberia_sf, shape = 21, col = "black",
          fill = "white") +
  scale_fill_viridis_c() +
  theme_minimal() +
  labs(title = "Combined Spatial Data: Points, Polygons, and Raster",
       fill = "Log population",
       x = "", y = "")

combined_map
Figure 2.9: Overlay of population per grid cell, admin-1 boundaries and survey locations in a single ggplot2 map.

Adding a scale bar and north arrow can improve the readability of spatial visualisations. The ggspatial package provides tools for integrating these elements into maps. The following example adds a scale bar and north arrow to a map containing raster, polygon and point data.

# Load necessary libraries
library(ggspatial)

# Add scale bar and north arrow
combined_map +
  annotation_scale(location = "bl", width_hint = 0.25) +
  annotation_north_arrow(location = "tr", which_north = "true",
                         height = unit(.5, "cm"),
                         width = unit(.5, "cm"))
Figure 2.10: Enhanced map with scale bar and north arrow added using the ggspatial package.

2.8 FAQs

  • How should I choose a CRS?

Start from the spatial operation and the property that matters, such as distance, area or direction. Preserve the source CRS when it is suitable. For local distance-based geostatistical modelling, use a projected CRS that represents the study region appropriately. A UTM zone is often useful for a compact area contained within that zone, but it is not a universal choice. Always inspect the transformed layers together before analysis.

  • Why did raster extraction return NA?

The location may be outside the raster extent, the source cell may be missing or preprocessing may have masked it. First check the CRS and extent, then plot the location with the raster and inspect the source data. Do not replace the value automatically. A nearby valid cell may be defensible for an isolated gap in a smooth continuous field, but not for locations outside the source coverage or across a meaningful barrier.

  • Should I extract a raster value at a point or summarise a buffer?

Use a point value when the covariate at the exact sampling location matches the scientific question. Use a buffer when the surrounding environment defines the relevant exposure. Choose the buffer radius using subject-matter knowledge. The summary must also match the raster values, for example a mean for a continuous field or density and a sum for population counts.

  • How should I choose a raster or prediction-grid resolution?

Resolution should reflect the source data, the spatial scale of the process, the intended output and the available computing resources. A finer grid can make a map look more detailed, but it does not create information beyond that supported by the data and model. When combining rasters, also check their CRS, extent, origin and alignment rather than comparing resolution alone.

  • Where can I find spatial datasets and learn more about working with them in R?

The sources listed in Table 2.1 are useful starting points, but always inspect the product definition, units, reference period, resolution and licence before use. Geocomputation with R provides a comprehensive introduction to spatial analysis in R, while rspatialdata catalogues packages and services for obtaining spatial data.

2.9 Review questions

  • What are the main differences between raster and vector data, and in which scenarios is each type commonly used in geostatistical analysis?

  • Why is it important to ensure that spatial datasets share the same coordinate reference system (CRS)? What problems can arise if CRSs are mismatched?

  • Why does RiskMap require suitable projected coordinates for distance-based modelling even though CRS-aware functions such as sf::st_distance() can calculate distances from longitude and latitude?

  • Describe the purpose of creating a predictive grid. How does the resolution of this grid affect the model?

  • What is the difference between extracting raster values at a point location versus using a buffer around a point?

2.10 Exercises

  1. Download a shapefile of administrative boundaries for a country of your choice using the geobounds package. Once the shapefile has been imported into R, visualise it in two different ways, first using base R plotting functions and then using ggplot2. Compare the two visualisations and discuss the advantages and disadvantages of each approach.

  2. Obtain raster data for the same country, such as elevation or population density. You may use either get_worldpop() from RiskMap or access Google Earth Engine through the rgee package. Plot the raster with terra::plot() and again with ggplot2. Comment on how the raster is displayed under each method and the flexibility provided by each plotting tool.

  3. Check the CRS and geographic extent of your vector and raster datasets, then choose a projected CRS suitable for the study area and analysis. A UTM zone may be appropriate for a compact study contained within that zone. Transform the datasets, plot them together and confirm that they align. Explain why assigning the source CRS and transforming coordinates are different operations.

  4. Generate 50 random points within the boundary of your chosen country. At each of these locations, extract the values of the raster dataset (for example elevation or population). Add these values to the point dataset as new attributes. Reflect on the usefulness of point-level covariates for modelling and prediction.

  5. Create a 2 km buffer around each of the random points generated in the previous exercise. Decide whether the raster represents an intensive quantity, an extensive quantity or a category, and choose an appropriate buffer summary. Account for cells that only partially overlap each buffer where necessary. Compare the buffer summaries with the point extractions, inspect missing coverage and edge effects, and explain how the two covariates answer different scientific questions.

  6. Construct prediction locations over the country using grid-centre points spaced approximately 5 km apart. Overlay the points on the raster background to check their coverage and extract raster covariates at each location. Explain the difference between these points, polygonal grid cells and raster pixels, and describe how the spacing affects computation and the displayed prediction surface.