What this tutorial does

This standalone practical explores 14 environmental and human-activity covariates for a spatial multi-criteria decision analysis (MCDA). We calculate summary statistics, draw maps and distribution plots, and compare covariates at shared locations to identify potentially redundant information.

We ask three questions:

What values occur across the study area? Summarize valid raster cells inside the country boundaries.
What values occur at recorded case locations? Summarize values at confirmed records, as a separate descriptive example.
Which covariates vary together across the study area? Compare values at the same spatial locations using correlation.

The correlation exercise uses study-area locations independently of the case records. It describes relationships between covariates, not relationships with an outcome. No existing MCDA outputs are needed.

Run the code in order, using terra and base R. The companion .R file contains the same code with explanatory comments. The HTML contains computed tables and figures. Summary statistics describe distributions; correlation helps us discuss redundancy when choosing criteria, alongside relevance and data quality.

1 Find the project and load the files

Start R in the project folder or one of its subfolders. The short loop below looks upward for 03Input, so the tutorial still works from its new subfolder. It does not change your working directory.

library(terra)
#> Warning: package 'terra' was built under R version 4.4.3
#> terra 1.9.34

project <- getwd()
while (!dir.exists(file.path(project, "03Input"))) {
  parent <- dirname(project)
  if (parent == project) stop("Start R inside the MCDA project folder.")
  project <- parent
}

input <- file.path(project, "03Input/01GisDatabase")
output <- file.path(project, "@covariate-summary/V2/outputs")
dir.create(file.path(output, "figures"), recursive = TRUE, showWarnings = FALSE)

files <- c(
  bio1_25m = "Bioclim/bio1_25m_clipped_Extent.tif",
  bio2_25m = "Bioclim/bio2_25m_clipped_Extent.tif",
  bio5_25m = "Bioclim/bio5_25m_clipped_Extent.tif",
  bio10_25m = "Bioclim/bio10_25m_clipped_Extent.tif",
  bio15_25m = "Bioclim/bio15_25m_clipped_Extent.tif",
  bio18_25m = "Bioclim/bio18_25m_clipped_Extent.tif",
  dem = "DEM/DEM_Extent.tif",
  ndvi = "NDVI/wg2414mx_clipped_Extent.tif",
  footprint = "Footprint/Human_Foot_Print_Extent.tif",
  worldpop = "Worldpop/WorldPop_PopDensity_Extent.tif",
  livestock_ctl = "GLW/GLW4-2020.D-DA.CTL_clipped_Extent.tif",
  livestock_shp = "GLW/GLW4-2020.D-DA.SHP_clipped_Extent.tif",
  corridors = "Biodiversity/Corridors/Corridors_Distance_Extent.tif",
  species_richness = "Biodiversity/Species_Richness/Combined_SR_2025_clipped_Extent.tif"
)

labels <- c(
  "BIO1 Annual mean temperature", "BIO2 Mean diurnal range",
  "BIO5 Maximum temperature of warmest month",
  "BIO10 Mean temperature of warmest quarter",
  "BIO15 Precipitation seasonality", "BIO18 Precipitation of warmest quarter",
  "Elevation DEM", "Maximum NDVI",
  "Human footprint", "Population density", "Cattle density", "Sheep density",
  "Distance to wildlife corridor", "Combined species richness"
)
names(labels) <- names(files)
paths <- file.path(input, "01Raster", files)
names(paths) <- names(files)
stopifnot(all(file.exists(paths)))

countries <- vect(file.path(input,
                 "02Vector/Administrative/Country_Boundaries_Extent_v2.shp"))
cases <- vect(file.path(input, "02Vector/RVF/RVF_EMPRESi_All_251030_Extent.shp"))
cases <- cases[which(cases$Diagnosi_1 == "Confirmed"), ]
if (any(!is.valid(countries))) countries <- makeValid(countries)
stopifnot(nrow(cases) > 0, all(is.valid(countries)))
nrow(cases)
#> [1] 261

Each layer is addressed by its name and filename. We retain the stored raster units: no temperature, NDVI or other product encoding is guessed or converted.

The study-area summaries below give every valid native cell equal weight. They are not area-weighted, and different native resolutions remain in use. The case summaries weight each record equally; repeated records at the same location remain repeated observations. A case distribution reflects recording and sampling as well as environmental conditions. It is not a case-control test.

2 Calculate familiar statistics for one variable

Start with human footprint as a worked example. extract() reads raster values at points. We project the points to the raster’s coordinate reference system (CRS) before extraction.

r <- rast(paths["footprint"])
case_points <- project(cases, crs(r))
footprint_cases <- extract(r, case_points)[, 2]
x <- footprint_cases[is.finite(footprint_cases)]

length(x)                     # Valid observations
#> [1] 259
mean(x)                       # Arithmetic average
#> [1] 1670.846
median(x)                     # Middle observation
#> [1] 1197
var(x)                        # Variance, using denominator n - 1
#> [1] 1578749
sd(x)                         # Standard deviation
#> [1] 1256.483
IQR(x)                        # Width of the middle 50%
#> [1] 1121.5
100 * sd(x) / mean(x)         # Arithmetic CV in percent
#> [1] 75.20041
100 * mad(x) / median(x)      # Robust relative-spread ratio in percent
#> [1] 59.20491

Read these measures carefully.

Measure	What it tells us	Important qualification
Mean and median	Centre of the distribution	A mean well above the median often signals a right tail.
Variance	Squared dispersion around the mean	Units are squared. Large numerical variance does not imply greater importance across different units.
Standard deviation	Dispersion in the original units	SD is the square root of variance.
IQR	Q75 minus Q25	Less influenced by extremes than the full range.
CV	100 × SD / mean	Relative spread; meaningful interpretation requires an appropriate ratio scale and a mean away from zero.
Robust CV used here	100 × MAD / median	R’s `mad()` uses a default scaling factor of 1.4826. This ratio also depends on a meaningful zero and a nonzero median.

R’s var() and sd() use n − 1, even when we supply all raster cells. We use that convention consistently throughout this tutorial; these numbers alone do not establish sampling uncertainty or independent observations. See R’s SD documentation.

For the supplied footprint case values, CV is approximately 75.20% and the robust ratio is approximately 59.20%. Compare these with the printed R results.

3 Use one short summary function

We calculate the same measures for every variable. is.finite() excludes NA, NaN and infinite values. Zero is a valid value unless product metadata say it is a fill code. Undefined ratios become NA, not zero.

percent_ratio <- function(top, bottom) {
  if (!is.finite(bottom) || bottom <= 0) return(NA_real_)
  100 * top / bottom
}

summarize_values <- function(x) {
  missing <- sum(!is.finite(x))
  x <- x[is.finite(x)]
  n <- length(x)

  if (n == 0) {
    return(data.frame(n = 0, missing = missing, min = NA, Q25 = NA,
      median = NA, mean = NA, Q75 = NA, max = NA, variance = NA,
      sd = NA, IQR = NA, MAD = NA, CV_percent = NA, robust_CV_percent = NA))
  }

  data.frame(n = n, missing = missing, min = min(x),
    Q25 = unname(quantile(x, 0.25)), median = median(x), mean = mean(x),
    Q75 = unname(quantile(x, 0.75)), max = max(x),
    variance = var(x), sd = sd(x), IQR = IQR(x), MAD = mad(x),
    CV_percent = percent_ratio(sd(x), mean(x)),
    robust_CV_percent = percent_ratio(mad(x), median(x)))
}

# These test the handling of zero denominators and missing values.
summarize_values(c(0, 0, 0))
#>   n missing min Q25 median mean Q75 max variance sd IQR MAD CV_percent
#> 1 3       0   0   0      0    0   0   0        0  0   0   0         NA
#>   robust_CV_percent
#> 1                NA
summarize_values(c(1, 2, 3, NA))
#>   n missing min Q25 median mean Q75 max variance sd IQR    MAD CV_percent
#> 1 3       1   1 1.5      2    2 2.5   3        1  1   1 1.4826         50
#>   robust_CV_percent
#> 1             74.13

We calculate arithmetic CV columns as a teaching exercise, but do not rank all covariates by CV. In particular:

BIO1, BIO5 and BIO10 are temperatures with an arbitrary zero if expressed in Celsius. Their CVs should not be interpreted as comparable relative variability.
BIO2 is a temperature range, which differs from an absolute temperature: its zero represents no diurnal range.
Elevation depends on its vertical datum; NDVI can include negative values and may be stored with a scale factor.
Density and distance may support a relative-spread interpretation if their units and zero are meaningful.
Index variables, including footprint and precipitation seasonality, require attention to the meaning of the index. The CV of BIO15 values is not the original coefficient of variation used to define BIO15.

NIST’s CV guidance explains why a meaningful zero matters. Changing Celsius to Fahrenheit can change a CV without changing the physical temperature pattern.

4 Draw a map with terra

The following function creates one PNG. It plots raw covariate values, overlays boundaries, and adds recorded cases. Each map has its own value scale; these are covariate maps, not standardized risk maps.

draw_map <- function(r, boundary, points, title, filename) {
  png(filename, width = 1600, height = 1100, res = 150)
  on.exit(dev.off())
  par(mar = c(4, 4, 4, 7))

  plot(r, main = title, col = hcl.colors(100, "YlGnBu", rev = TRUE),
       maxcell = 250000,
       plg = list(title = "Stored value"))
  lines(boundary, col = "grey30", lwd = 0.7)
  points(points, pch = 1, cex = 0.45, col = "black")
  mtext("Circles: confirmed records | Native values; display resolution reduced",
        side = 1, line = 2.7, cex = 0.8)
}

For display, terra limits the number of plotted cells with maxcell. The summary calculations still use all valid native cells inside the mask. A display approximation must not be confused with the values used for statistics. See terra plot.

5 Draw a density plot and boxplots

A density curve is a smoothed picture of the value distribution. A boxplot shows the median, the middle half of observations and the tails. Together they show more than a mean and SD alone.

Each figure compares a random sample of study-area cells with all valid case records. Sampling is only for the regional graphics, to keep plotting fast. We use up to 20,000 valid cells per variable; the summary table is not sampled. Case points are shown on the boxplot, with small horizontal jitter so overlapping records remain visible.

draw_distributions <- function(region, cases_x, title, filename) {
  region <- region[is.finite(region)]
  cases_x <- cases_x[is.finite(cases_x)]
  png(filename, width = 1600, height = 1050, res = 150)
  on.exit(dev.off())
  par(mfrow = c(1, 2), mar = c(5, 5, 4, 1), oma = c(2, 0, 2, 0))

  # A curve needs at least two observations and some variation.
  curves <- list()
  if (length(region) > 1 && sd(region) > 0) curves$region <- density(region)
  if (length(cases_x) > 1 && sd(cases_x) > 0) curves$cases <- density(cases_x)

  if (length(curves) > 0) {
    x_limits <- range(region, cases_x)
    y_max <- max(sapply(curves, function(d) max(d$y)))
    plot(NA, xlim = x_limits, ylim = c(0, 1.2 * y_max),
         xlab = "Stored raster value", ylab = "Probability density",
         main = "Distribution shape")
    if (!is.null(curves$region)) lines(curves$region, col = "steelblue", lwd = 2)
    if (!is.null(curves$cases)) lines(curves$cases, col = "darkorange3", lwd = 2)
    legend("topright", legend = c("Study-area sample", "Case records"),
           col = c("steelblue", "darkorange3"), lwd = 2, bty = "n", cex = 0.8)
  } else {
    plot.new()
    text(0.5, 0.5, "No density curve: constant or insufficient values")
  }

  boxplot(list("Study area" = region, "Case records" = cases_x),
          col = c("lightblue", "moccasin"), ylab = "Stored raster value",
          main = "Centre and spread", outline = TRUE, pch = 20, cex = 0.35)
  if (length(cases_x)) {
    points(jitter(rep(2, length(cases_x)), amount = 0.12), cases_x,
           pch = 16, cex = 0.4, col = adjustcolor("black", alpha.f = 0.35))
  }
  mtext(title, outer = TRUE, font = 2, cex = 1.1)
  mtext(paste("Study-area sample n =", length(region),
              "| valid case records n =", length(cases_x)),
        outer = TRUE, side = 1, line = 0.5, cex = 0.9)
}

Density uses R’s default Gaussian kernel and automatic bandwidth. Bandwidth changes the smoothness; small bumps are not automatically meaningful clusters. Gaussian smoothing can extend beyond a bounded variable’s support. We display only the observed value range without applying boundary correction, so the visible part of a curve need not integrate to exactly one. Each full curve is normalized separately: its height is not the number of cases or a disease probability.

In a base R boxplot, the whiskers reach the most extreme observations within 1.5 IQR of the hinges. Points beyond the whiskers are potential outliers, not necessarily errors. Base R hinges can differ slightly from the quantile() values in small samples. See R boxplots.

6 Process all 14 variables

For each file, we keep two distinct sets of values:

regional_values: cells inside the country mask, at native resolution;
case_values: extracted from the original raster at confirmed record locations.

Case extraction uses the original raster at each record location. Polygon masking uses cell centres (touches = FALSE). The regional missing count includes cells outside the polygons in the cropped rectangle, so it is not an estimate of missing coverage within the study area.

regional_summary <- data.frame()
case_summary <- data.frame()
extracted_cases <- as.data.frame(cases)
plot_manifest <- data.frame()
set.seed(2026)  # Repeatable random samples and jitter.

for (name in names(paths)) {
  r <- rast(paths[name])
  boundary <- project(countries, crs(r))
  point_layer <- project(cases, crs(r))

  regional_raster <- mask(crop(r, boundary), boundary, touches = FALSE)
  regional_values <- values(regional_raster, mat = FALSE)
  case_values <- extract(r, point_layer)[, 2]
  extracted_cases[[name]] <- case_values

  region_row <- summarize_values(regional_values)
  region_row$covariate <- name
  regional_summary <- rbind(regional_summary, region_row)
  case_row <- summarize_values(case_values)
  case_row$covariate <- name
  case_summary <- rbind(case_summary, case_row)

  # Graphics use a sample of valid regional values; statistics above use all.
  valid_region <- regional_values[is.finite(regional_values)]
  if (!length(valid_region)) stop("No valid study-area cells for ", name)
  chosen <- sample.int(length(valid_region), min(20000, length(valid_region)))
  regional_sample <- valid_region[chosen]

  map_file <- file.path(output, "figures", paste0("map_", name, ".png"))
  dist_file <- file.path(output, "figures", paste0("distribution_", name, ".png"))
  draw_map(regional_raster, boundary, point_layer, labels[name], map_file)
  draw_distributions(regional_sample, case_values, labels[name], dist_file)

  plot_manifest <- rbind(plot_manifest, data.frame(covariate = name,
    label = labels[name], regional_plot_n = length(regional_sample),
    case_plot_n = sum(is.finite(case_values)), map = map_file, distribution = dist_file))
  rm(regional_values, valid_region)
  message("Finished ", name)
}
#> Finished bio1_25m
#> Finished bio2_25m
#> Finished bio5_25m
#> Finished bio10_25m
#> Finished bio15_25m
#> Finished bio18_25m
#> Finished dem
#> Finished ndvi
#> Finished footprint
#> Finished worldpop
#> Finished livestock_ctl
#> Finished livestock_shp
#> Finished corridors
#> Finished species_richness

write.csv(regional_summary, file.path(output, "regional_summary.csv"), row.names = FALSE)
write.csv(case_summary, file.path(output, "case_summary.csv"), row.names = FALSE)
write.csv(extracted_cases, file.path(output, "case_values.csv"), row.names = FALSE)
write.csv(plot_manifest, file.path(output, "plot_manifest.csv"), row.names = FALSE)

The footprint raster is large. Exact quantiles and MAD require reading and processing many values, so the full run can take several minutes and several GB of memory. Do not silently substitute sampled statistics for these exact native-cell summaries. For a smaller classroom exercise, process one variable first by changing the loop to for (name in "footprint"); restore the full loop before the later 14-variable checks.

7 Read the results

At confirmed case records, the following table describes the extracted values. CV is an arithmetic output, subject to the interpretation cautions above.

case_display <- case_summary[, c("covariate", "n", "missing", "mean", "median",
                                 "sd", "variance", "CV_percent")]
case_display[, 4:8] <- round(case_display[, 4:8], 2)
case_display
#>           covariate   n missing     mean   median       sd     variance
#> 1          bio1_25m 261       0    28.37    28.95     1.60 2.570000e+00
#> 2          bio2_25m 261       0    13.91    14.26     1.06 1.120000e+00
#> 3          bio5_25m 261       0    40.12    41.70     2.94 8.660000e+00
#> 4         bio10_25m 261       0    32.23    33.31     2.09 4.370000e+00
#> 5         bio15_25m 261       0   142.18   146.45    19.35 3.742500e+02
#> 6         bio18_25m 261       0    67.43    39.00    67.20 4.516170e+03
#> 7               dem 261       0   172.54    96.00   166.80 2.782088e+04
#> 8              ndvi 261       0  2823.23  2603.00  1291.86 1.668905e+06
#> 9         footprint 259       2  1670.85  1197.00  1256.48 1.578749e+06
#> 10         worldpop 261       0   491.08    24.59  2569.77 6.603726e+06
#> 11    livestock_ctl 257       4    15.05    10.69    19.06 3.634100e+02
#> 12    livestock_shp 257       4    37.20    22.51    73.65 5.424280e+03
#> 13        corridors 261       0 94860.30 53326.08 97176.20 9.443213e+09
#> 14 species_richness 261       0   282.24   219.43   168.37 2.834702e+04
#>    CV_percent
#> 1        5.65
#> 2        7.62
#> 3        7.34
#> 4        6.49
#> 5       13.61
#> 6       99.66
#> 7       96.67
#> 8       45.76
#> 9       75.20
#> 10     523.29
#> 11     126.70
#> 12     197.97
#> 13     102.44
#> 14      59.65

Across the study-area raster cells, the distribution can be very different. Here we mask to the country boundaries first, so the summary represents the defined study area rather than the whole rectangular raster file.

region_display <- regional_summary[, c("covariate", "n", "mean", "median",
                                       "sd", "variance", "CV_percent")]
region_display[, 3:7] <- round(region_display[, 3:7], 2)
region_display
#>           covariate         n      mean   median        sd     variance
#> 1          bio1_25m    419907     26.54    26.90      2.08 4.330000e+00
#> 2          bio2_25m    419907     13.49    14.27      2.12 4.500000e+00
#> 3          bio5_25m    419907     39.27    40.05      3.98 1.585000e+01
#> 4         bio10_25m    419907     31.52    32.00      3.21 1.030000e+01
#> 5         bio15_25m    419907     99.60   101.84     38.73 1.500140e+03
#> 6         bio18_25m    419907     89.14    34.00    114.76 1.316877e+04
#> 7               dem  10498022    387.61   338.00    244.40 5.973166e+04
#> 8              ndvi  10136228   3681.20  1855.00   2877.48 8.279866e+06
#> 9         footprint 100099315    712.11   501.00    756.04 5.715930e+05
#> 10         worldpop  10484389     51.57     1.97    428.71 1.837909e+05
#> 11    livestock_ctl    104731     11.01     0.49     22.22 4.935800e+02
#> 12    livestock_shp    104730     16.42     3.15     38.22 1.460530e+03
#> 13        corridors     72909 140509.45 56692.62 172398.14 2.972112e+10
#> 14 species_richness     73373    303.18   180.52    260.66 6.794533e+04
#>    CV_percent
#> 1        7.84
#> 2       15.72
#> 3       10.14
#> 4       10.18
#> 5       38.89
#> 6      128.74
#> 7       63.05
#> 8       78.17
#> 9      106.17
#> 10     831.30
#> 11     201.70
#> 12     232.71
#> 13     122.70
#> 14      85.98

A narrow case distribution within a broad regional distribution is a descriptive pattern worth discussing. It does not by itself demonstrate environmental preference or discrimination: there are no surveyed absence/control locations, and reporting effort is not modeled. Neither SD nor CV captures geographic clustering; two maps can have identical value summaries and different patterns.

For a robust view of the case distributions, inspect the quartiles and MAD too.

robust_display <- case_summary[, c("covariate", "Q25", "median", "Q75", "IQR",
                                   "MAD", "robust_CV_percent")]
robust_display[, 2:7] <- round(robust_display[, 2:7], 2)
robust_display
#>           covariate      Q25   median       Q75      IQR      MAD
#> 1          bio1_25m    27.49    28.95     29.47     1.98     1.33
#> 2          bio2_25m    13.65    14.26     14.34     0.70     0.40
#> 3          bio5_25m    37.93    41.70     42.19     4.26     2.05
#> 4         bio10_25m    30.28    33.31     33.76     3.48     1.45
#> 5         bio15_25m   138.75   146.45    151.79    13.04     8.20
#> 6         bio18_25m    24.00    39.00     64.00    40.00    22.24
#> 7               dem    15.00    96.00    368.00   353.00   133.43
#> 8              ndvi  2004.00  2603.00   3537.00  1533.00  1000.76
#> 9         footprint   948.00  1197.00   2069.50  1121.50   708.68
#> 10         worldpop     1.32    24.59    114.88   113.56    34.83
#> 11    livestock_ctl     6.45    10.69     19.96    13.51    10.72
#> 12    livestock_shp     6.18    22.51     41.49    35.31    26.02
#> 13        corridors 33853.10 53326.08 131741.23 97888.13 43381.84
#> 14 species_richness   131.87   219.43    438.18   306.32   183.94
#>    robust_CV_percent
#> 1               4.59
#> 2               2.80
#> 3               4.92
#> 4               4.35
#> 5               5.60
#> 6              57.02
#> 7             138.99
#> 8              38.45
#> 9              59.20
#> 10            141.65
#> 11            100.27
#> 12            115.58
#> 13             81.35
#> 14             83.83

8 Inspect the maps and distributions

For each variable, first find where low and high values occur on the map. Then look at the density and boxplots. Ask whether the case values span the whole regional distribution, whether a tail dominates the mean, and whether the plot is compressed by a few extremes. Each map retains its own stored-value scale. The exported PNG files can be inserted into QGIS Print Layout using Add Picture.