Before you start

This introductory practical uses one covariate: precipitation seasonality (BIO15). You will load and inspect a raster, load administrative boundaries and a disease shapefile, prepare a study-area raster, calculate summaries, and make simple plots.

No sampling is needed here. The supplied BIO15 file has 621,504 cells. A numeric vector of that size takes only about 5 MB, so using all its values is practical on a typical teaching computer. V4 uses every valid study-area cell for the statistics, histogram, density estimate and boxplot. It also summarizes all usable confirmed records. The analysis remains descriptive and uses one raster.

Run one code block at a time. There are no custom functions or loops. You only need terra and base R for the analysis.

In RStudio, open the course project with the folder containing 03Input as your working directory. Alternatively, choose Session → Set Working Directory → Choose Directory and select that folder. Keep all shapefile components (.shp, .shx, .dbf, .prj, and any others) together; you load the .shp file.

# Run install.packages("terra") once if terra is not already installed.
library(terra)
Warning: package 'terra' was built under R version 4.4.3
terra 1.9.34
getwd()                          # Where is R looking for files?
[1] "/Users/mpetric/Documents/Documents/@RND/OHS/OHS-internshp/covariates/mcda-thresholds-and-parameters"
dir.exists("03Input")            # This should be TRUE.
[1] TRUE
stopifnot(dir.exists("03Input"))

output <- "@covariate-summary/V4/outputs_bio15"
dir.create(output, recursive = TRUE, showWarnings = FALSE)

The HTML contains worked results. The companion .R script contains the same analysis commands and explanatory comments. Outputs are saved in the V4 folder.

1 Load and inspect the raster

A raster is a grid of cells. Each cell in this file contains a precipitation-seasonality value. rast() opens the file as a SpatRaster; it does not immediately load all cell values into memory.

bio15 <- rast("03Input/01GisDatabase/01Raster/Bioclim/bio15_25m_clipped_Extent.tif")
names(bio15) <- "precipitation_seasonality"  # Give the layer a readable name.

bio15                         # Overview of the raster.
class       : SpatRaster
size        : 624, 996, 1  (nrow, ncol, nlyr)
resolution  : 0.04166667, 0.04166667  (x, y)
extent      : -25.5, 16, 2, 28  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source      : bio15_25m_clipped_Extent.tif
name        : precipitation_seasonality
min value   :                 16.372847
max value   :                186.730728
nlyr(bio15)                    # Number of layers: one.
[1] 1
res(bio15)                     # Cell width and height in CRS units.
[1] 0.04166667 0.04166667
ext(bio15)                     # Geographic extent.
SpatExtent : -25.5, 16, 2, 28 (xmin, xmax, ymin, ymax)
crs(bio15, describe = TRUE)    # Coordinate reference system (CRS).
    name authority code  area             extent
1 WGS 84      EPSG 4326 World -180, 180, -90, 90

Read the printed output before continuing. This file uses longitude and latitude, so its resolution is in degrees, not metres. Here 0.04166667 degrees equals 2.5 arcminutes; the filename does not mean 25-metre cells. The values describe precipitation seasonality, not rainfall amounts in millimetres.

What is BIO15? WorldClim defines BIO15 as precipitation seasonality, expressed as a coefficient of variation. Conceptually, it compares variation in monthly precipitation with mean monthly precipitation, usually expressed as a percentage. Low values indicate relatively even monthly precipitation; high values indicate stronger seasonality. Values can exceed 100.

The raster already contains BIO15. We do not calculate it again from monthly rainfall in this lesson. A mean or SD across its cells describes how seasonality varies between places. Calculating 100 * sd(cell_values) / mean(cell_values) would describe relative spatial spread; it would not recreate the monthly CV that defines BIO15. We therefore keep our summaries to familiar descriptive statistics below and retain the supplied raster values without rescaling them.

2 Load administrative boundaries and disease records

vect() opens a vector file as a SpatVector. Boundaries are polygons; the disease records are points. Vector layers also have an attribute table.

boundaries <- vect("03Input/01GisDatabase/02Vector/Administrative/Country_Boundaries_Extent_v2.shp")
disease <- vect("03Input/01GisDatabase/02Vector/RVF/RVF_EMPRESi_All_251030_Extent.shp")

geomtype(boundaries)
[1] "polygons"
geomtype(disease)
[1] "points"
head(as.data.frame(boundaries)[, c("NAME", "ISO3")])
                      NAME ISO3
1                  Algeria  DZA
2                    Benin  BEN
3             Burkina Faso  BFA
4                 Cameroon  CMR
5               Cape Verde  CPV
6 Central African Republic  CAF
names(disease)                    # Names of the attribute columns.
 [1] "Event.ID"   "Disease"    "Serotype"   "Region"     "Subregion" 
 [6] "Country"    "Admin.leve" "Locality"   "Latitude"   "Longitude" 
[11] "Diagnosis." "Diagnosi_1" "Animal.typ" "Species"    "Observatio"
[16] "Report.dat" "Humans.aff" "Human.deat" "field_19"   "observat_1"
[21] "reportDate"
table(disease$Diagnosi_1)          # Counts by diagnostic status.

Confirmed    Denied 
      261         4 
# Keep only records whose status is Confirmed.
disease <- disease[which(disease$Diagnosi_1 == "Confirmed"), ]
nrow(disease)
[1] 261

$ selects an attribute column. The square brackets select rows. In the supplied file, this keeps 261 confirmed records and excludes four denied records.

3 Match the coordinate reference systems

Layers must use compatible coordinates to overlay or extract values correctly. We transform the two vector layers into the raster’s CRS. The raster stays on its original grid.

boundaries <- project(boundaries, crs(bio15))
disease <- project(disease, crs(bio15))

project() transforms coordinates. Simply assigning a different CRS label does not transform them. If these files already share a CRS, this step leaves their locations unchanged.

4 Crop and mask the raster

crop() cuts the raster to a rectangle covering the boundaries. mask() then sets cells outside the polygons to NA. NA means missing or excluded; it is different from a valid zero value.

bio15_crop <- crop(bio15, boundaries)
bio15_study <- mask(bio15_crop, boundaries, touches = FALSE)
bio15_study
class       : SpatRaster
size        : 624, 993, 1  (nrow, ncol, nlyr)
resolution  : 0.04166667, 0.04166667  (x, y)
extent      : -25.375, 16, 2, 28  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source(s)   : memory
varname     : bio15_25m_clipped_Extent
name        : precipitation_seasonality
min value   :                 16.372847
max value   :                186.730728

Here touches = FALSE retains cells whose centres fall inside a boundary. The original bio15 object is still available. Creating a new object with <- lets you inspect the result without replacing the original object.

5 Draw the map

Draw the raster first, then add the polygon outlines and disease points. png() starts saving a figure; dev.off() finishes the file. The same plotting commands can be used without those two lines to draw in RStudio’s Plots pane.

png(file.path(output, "bio15_map.png"), width = 1500, height = 1100, res = 150)
plot(bio15_study, main = "Precipitation seasonality (BIO15) and confirmed records",
     col = hcl.colors(100, "YlGnBu", rev = TRUE), maxcell = ncell(bio15_study),
     plg = list(title = "BIO15"))
lines(boundaries, col = "grey30", lwd = 0.7)
points(disease, pch = 1, cex = 0.45, col = "black")
dev.off()
quartz_off_screen 
                2

Circles mark confirmed records. Setting maxcell to the raster’s cell count allows terra to plot the full native grid without reducing the number of input cells for display. The exported image still has a finite pixel resolution.

Precipitation seasonality (BIO15) map
Precipitation seasonality (BIO15) map

6 Read all valid cells and calculate full summaries

values() reads the raster values into R. mat = FALSE gives us a simple vector because this raster has one layer. is.finite() removes NA and any non-finite values. Valid zeros are retained.

x <- values(bio15_study, mat = FALSE)
x <- x[is.finite(x)]
length(x)                         # All valid study-area cells; no sample.
[1] 419907
summary(x)                        # Minimum, quartiles, median, mean, maximum.
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  16.37   66.13  101.84   99.60  134.47  186.73 
sd(x)                             # Standard deviation.
[1] 38.73166
var(x)                            # Variance.
[1] 1500.141
IQR(x)                            # Width of the middle 50%.
[1] 68.34433
mad(x)                            # Robust measure of spread.
[1] 50.46549

There are 419,907 valid study-area cells in the supplied data. Cells outside the country polygons or with missing raster values are excluded. “Full” means all valid cells within this defined study area, not values outside it.

global() provides the same basic summaries directly from the raster. We add quartiles and robust measures from the full vector to make one exportable table.

raster_summary <- global(bio15_study,
                         c("notNA", "min", "max", "mean", "sd"), na.rm = TRUE)
raster_summary$Q25 <- unname(quantile(x, 0.25))
raster_summary$median <- median(x)
raster_summary$Q75 <- unname(quantile(x, 0.75))
raster_summary$variance <- var(x)
raster_summary$IQR <- IQR(x)
raster_summary$MAD <- mad(x)
round(raster_summary, 2)
                           notNA   min    max mean    sd   Q25 median    Q75
precipitation_seasonality 419907 16.37 186.73 99.6 38.73 66.13 101.84 134.47
                          variance   IQR   MAD
precipitation_seasonality  1500.14 68.34 50.47
Output Meaning
notNA Number of valid cells
min / max Smallest / largest value
mean / median Average / middle value
Q25 / Q75 Lower and upper quartiles
sd Standard deviation: spread in the original units
variance SD squared: spread in squared units
IQR Q75 minus Q25: width of the middle half of values
MAD Median absolute deviation from the median, multiplied by R’s default scaling factor of 1.4826

These summaries give every valid cell equal weight. Longitude–latitude cells do not all have equal ground area, so this is a cell-based summary, not an area-weighted one. R’s sd() and var() use denominator n − 1. That is a formula convention; it does not mean we selected a sample of cells here. Adjacent cells can be similar, so the cell count is not a count of independent observations.

For a much larger raster, reading every value might use too much memory. Then terra’s built-in global() summaries can still be useful, and sampling may help with plots. This small BIO15 file does not need that compromise.

7 Plot the full distribution

Use the same vector x in all three plots. Every valid study-area cell contributes. There is no call to spatSample() and no selection of rows for plotting.

png(file.path(output, "bio15_distribution.png"),
    width = 1800, height = 650, res = 150)
par(mfrow = c(1, 3))               # One row with three panels.
hist(x, breaks = 40, col = "lightblue", main = "Histogram: all valid cells",
     xlab = "BIO15")
plot(density(x), main = "Density: all valid cells", xlab = "BIO15",
     xlim = range(x), xaxs = "i", col = "steelblue", lwd = 2)
boxplot(x, col = "lightblue", main = "Boxplot: all valid cells",
        ylab = "BIO15")
dev.off()
quartz_off_screen 
                2

Look at the centre, spread and tails. A long tail can pull the mean away from the median. Ask whether most cells have low, intermediate or high seasonality. These plots describe the value distribution, not geographic clustering.

Distribution of all valid raster cells
Distribution of all valid raster cells

8 Extract BIO15 values at disease records

extract() answers: what value does the raster contain at each point? With method = "simple", it reads the cell containing the point. We use the study-area raster, so excluded or missing cells return NA. Near a coast, a point can fall in a cell whose centre was outside the polygons and was therefore masked. With the supplied files, 259 of the 261 records have a usable value after masking.

point_values <- extract(bio15_study, disease, method = "simple")
head(point_values)
  ID precipitation_seasonality
1  1                  152.2462
2  2                  152.6683
3  3                  155.2011
4  4                  152.7672
5  5                  151.0690
6  6                  153.4597
# The first column is the point ID; the second contains BIO15 values.
disease$bio15 <- point_values[, 2]
y <- disease$bio15
sum(is.na(y))                      # Records without a raster value.
[1] 2
y <- y[is.finite(y)]
length(y)                         # Records with a usable value.
[1] 259
summary(y)                        # Minimum, quartiles, median, mean, maximum.
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  38.49  138.47  146.40  142.07  151.79  164.61 
sd(y)
[1] 19.37869
var(y)
[1] 375.5336
IQR(y)
[1] 13.31913
mad(y)
[1] 8.265468
record_summary <- data.frame(n = length(y), missing = sum(!is.finite(disease$bio15)),
  min = min(y), Q25 = unname(quantile(y, 0.25)), median = median(y),
  mean = mean(y), Q75 = unname(quantile(y, 0.75)), max = max(y),
  sd = sd(y), variance = var(y), IQR = IQR(y), MAD = mad(y))
round(record_summary, 2)
    n missing   min    Q25 median   mean    Q75    max    sd variance   IQR
1 259       2 38.49 138.47  146.4 142.07 151.79 164.61 19.38   375.53 13.32
   MAD
1 8.27

These statistics describe record locations, not all cells. Each record counts once; repeated records at the same place still count separately. Differences between these values and regional summaries can reflect recording effort as well as environmental patterns. Extraction alone does not show that BIO15 causes disease or predicts where it will occur.

png(file.path(output, "bio15_at_records.png"), width = 1800, height = 650, res = 150)
par(mfrow = c(1, 3))
hist(y, breaks = 20, col = "moccasin", main = "Histogram: all usable records",
     xlab = "BIO15")
plot(density(y), main = "Density: all usable records", xlab = "BIO15",
     xlim = range(y), xaxs = "i", col = "darkorange3", lwd = 2)
boxplot(y, col = "moccasin", ylab = "BIO15",
        main = "Boxplot: all usable records")
dev.off()
quartz_off_screen 
                2
BIO15 at confirmed records
BIO15 at confirmed records

9 Save your results

write.csv(raster_summary, file.path(output, "raster_summary.csv"), row.names = FALSE)
write.csv(record_summary, file.path(output, "record_summary.csv"), row.names = FALSE)
write.csv(as.data.frame(disease), file.path(output, "disease_with_bio15.csv"),
          row.names = FALSE)

The three PNGs and three CSVs are now in @covariate-summary/V4/outputs_bio15. The disease CSV is an attribute table, not a spatial file. To save the raster and points for QGIS too, you can optionally run:

writeRaster(bio15_study, file.path(output, "bio15_study.tif"), overwrite = TRUE)
writeVector(disease, file.path(output, "disease_with_bio15.gpkg"), overwrite = TRUE)

The GeoTIFF stores your cropped and masked raster. overwrite = TRUE replaces an existing file with that name. These two optional exports are not run in the companion script.

Remember these commands

Task Command
Load a raster rast()
Load polygons or points vect()
Inspect resolution, extent and CRS res(), ext(), crs()
Transform vector coordinates project()
Cut to a rectangle crop()
Exclude cells outside polygons mask()
Read every raster value values()
Summarize a raster global()
Summarize a vector of values summary(), sd(), var(), IQR(), mad()
Read raster values at points extract()
Draw the layers plot(), lines(), points()

Try these questions before looking back at the code:

  1. What are the raster’s resolution units?
  2. Why do we use both crop() and mask()?
  3. What is the difference between zero and NA?
  4. How can you check that the plots use every valid study-area cell? Why is reading all values practical for this raster?
  5. How many confirmed records have a usable BIO15 value?
  6. What is the difference between precipitation seasonality at one cell and the SD of BIO15 across many cells?
  7. Load the three input files again, then make the map using only the command table as a reminder.