Explore covariate distributions and spatial correlations

Standalone MCDA teaching practical · V2 · terra and base R

Runnable R script · Editable R Markdown

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:

1. What values occur across the study area? Summarize valid raster cells inside the country boundaries.
2. What values occur at recorded case locations? Summarize values at confirmed records, as a separate descriptive example.
3. 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.

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)
}

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.

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

9 Optional view of a strongly skewed variable

Population density may have a long right tail. Keeping the original plot is important, but a second view on log1p(x) can help reveal the low end without removing high values. log1p(x) means log(1 + x) and is valid here because these extracted density values are nonnegative.

x <- extracted_cases$worldpop
x <- x[is.finite(x)]
stopifnot(all(x >= 0))
png(file.path(output, "figures", "worldpop_case_log_view.png"),
    width = 1400, height = 650, res = 140)
par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
boxplot(x, main = "Original population density", ylab = "Stored value", col = "moccasin")
boxplot(log1p(x), main = "Same records after log1p", ylab = "log(1 + stored value)", col = "lightblue")
dev.off()
#> quartz_off_screen 
#>                 2

This supplementary plot changes the display variable, not the original summary tables or MCDA. Do not interpret its spread as being in the original density units.

10 What does correlation between maps mean?

If two covariates tend to have high values at the same locations and low values at the same locations, they have positive correlation. If one tends to be high where the other is low, they have negative correlation. Both can indicate overlapping information for MCDA. Correlation is unitless and lies between −1 and +1.

Spearman is a useful first view for skewed environmental layers. We also compute Pearson and inspect scatterplots. Neither measure establishes causation or importance to the decision. R computes Spearman by correlating ranks; tied values receive average ranks. See R's correlation documentation.

Two different spatial questions. Here we measure correlation between covariates at the same locations. Spatial autocorrelation asks whether nearby locations have similar values within one covariate. We do not calculate that second quantity here. Nearby observations are often similar, so our many points are not independent replicates. We use correlations descriptively and do not report ordinary correlation-test p-values.

11 Extract every covariate at shared locations

Never correlate two independent samples, or join native raster value vectors by row number: their cells may represent different places. Instead, create one set of locations and extract every covariate there.

We use centres of a 20 km grid in an equal-area CRS (EPSG:6933) inside the country boundaries. This spreads locations approximately evenly per unit area. It avoids giving more weight to a layer simply because its pixels are smaller. The grid spacing is a classroom choice, not an estimated independence distance.

These are point extractions, not 20 km averages. We read the native cell that contains each point with method = "simple"; the rasters are not resampled or aggregated. Native resolution, missingness and the chosen point spacing still influence the results. Boundary cells are included when their centre is inside a country. See terra extraction.

boundary_equal_area <- project(countries, "EPSG:6933")
grid_spacing_m <- 20000
sampling_grid <- rast(ext(boundary_equal_area), resolution = grid_spacing_m,
                      crs = crs(boundary_equal_area))
inside <- rasterize(boundary_equal_area, sampling_grid, field = 1,
                    background = NA, touches = FALSE)
locations <- as.points(inside, values = FALSE, na.rm = TRUE)

# One row per location; one column per covariate.
spatial_values <- data.frame(location_id = seq_len(nrow(locations)))
for (name in names(paths)) {
  r <- rast(paths[name])
  p <- project(locations, crs(r))
  spatial_values[[name]] <- extract(r, p, method = "simple", ID = FALSE)[, 1]
}

x_all <- spatial_values[, names(paths)]
x_all[] <- lapply(x_all, function(x) { x[!is.finite(x)] <- NA; x })
coverage <- data.frame(covariate = names(paths),
                       valid_locations = colSums(!is.na(x_all)),
                       missing_locations = colSums(is.na(x_all)))
coverage
#>                         covariate valid_locations missing_locations
#> bio1_25m                 bio1_25m           21438                 0
#> bio2_25m                 bio2_25m           21438                 0
#> bio5_25m                 bio5_25m           21438                 0
#> bio10_25m               bio10_25m           21438                 0
#> bio15_25m               bio15_25m           21438                 0
#> bio18_25m               bio18_25m           21438                 0
#> dem                           dem           21438                 0
#> ndvi                         ndvi           21438                 0
#> footprint               footprint           21381                57
#> worldpop                 worldpop           21410                28
#> livestock_ctl       livestock_ctl           21370                68
#> livestock_shp       livestock_shp           21370                68
#> corridors               corridors           21438                 0
#> species_richness species_richness           21438                 0

# Use the same complete rows for every pair of covariates.
keep <- complete.cases(x_all)
x_shared <- x_all[keep, , drop = FALSE]
stopifnot(nrow(x_shared) > 2, all(vapply(x_shared, sd, numeric(1)) > 0))

sampling_summary <- data.frame(grid_spacing_m = grid_spacing_m,
  candidate_locations = nrow(x_all), complete_locations = sum(keep),
  excluded_locations = sum(!keep), retained_percent = 100 * mean(keep))
sampling_summary
#>   grid_spacing_m candidate_locations complete_locations excluded_locations
#> 1          20000               21438              21336                102
#>   retained_percent
#> 1         99.52421

coordinates <- as.data.frame(crds(project(locations, "EPSG:4326")))
names(coordinates) <- c("longitude", "latitude")
write.csv(cbind(spatial_values, coordinates, complete = keep),
          file.path(output, "correlation_locations.csv"), row.names = FALSE)
write.csv(coverage, file.path(output, "correlation_coverage.csv"), row.names = FALSE)
write.csv(sampling_summary, file.path(output, "correlation_sampling_summary.csv"),
          row.names = FALSE)

Complete rows make every coefficient refer to the same geographic sample. However, excluding a location with any missing covariate may remove particular parts of the region. Read the coverage table and inspect the map below. If many locations are excluded, investigate coverage before using the correlations. When you change the set of variables, this shared footprint may change too.

png(file.path(output, "figures", "correlation_locations.png"),
    width = 1500, height = 1100, res = 150)
plot(boundary_equal_area, col = "grey97", border = "grey50",
     mar = c(6, 5, 4, 2),
     main = "Shared locations for covariate correlation")
points(locations[keep, ], pch = 16, cex = 0.16, col = "steelblue")
if (any(!keep)) points(locations[!keep, ], pch = 16, cex = 0.22, col = "darkorange3")
legend(x = xmin(boundary_equal_area) + 0.03 * diff(c(xmin(boundary_equal_area), xmax(boundary_equal_area))),
       y = ymin(boundary_equal_area) + 0.15 * diff(c(ymin(boundary_equal_area), ymax(boundary_equal_area))),
       legend = c("All 14 values available", "One or more values missing"),
       col = c("steelblue", "darkorange3"), pch = 16, bty = "n", cex = 0.85)
mtext(paste("Equal-area coordinates in metres; grid spacing", grid_spacing_m / 1000, "km"),
      side = 1, line = 4, cex = 0.8)
dev.off()
#> quartz_off_screen 
#>                 2

12 Calculate and plot the correlation matrices

cor() compares all pairs of columns. No standardization is needed for these coefficients: changing a variable's units by a positive linear conversion does not change its correlation. Nonlinear transformations can change Pearson; strictly increasing transformations preserve Spearman ranks.

pearson <- cor(x_shared, method = "pearson")
spearman <- cor(x_shared, method = "spearman")
write.csv(pearson, file.path(output, "correlation_pearson.csv"))
write.csv(spearman, file.path(output, "correlation_spearman.csv"))
round(spearman, 2)
#>                  bio1_25m bio2_25m bio5_25m bio10_25m bio15_25m bio18_25m   dem
#> bio1_25m             1.00     0.25     0.54      0.55      0.70      0.06 -0.50
#> bio2_25m             0.25     1.00     0.76      0.69      0.21     -0.62  0.22
#> bio5_25m             0.54     0.76     1.00      0.98      0.23     -0.67 -0.04
#> bio10_25m            0.55     0.69     0.98      1.00      0.20     -0.69 -0.04
#> bio15_25m            0.70     0.21     0.23      0.20      1.00      0.17 -0.24
#> bio18_25m            0.06    -0.62    -0.67     -0.69      0.17      1.00 -0.25
#> dem                 -0.50     0.22    -0.04     -0.04     -0.24     -0.25  1.00
#> ndvi                -0.01    -0.62    -0.66     -0.68      0.05      0.85 -0.20
#> footprint            0.10    -0.48    -0.52     -0.55      0.19      0.69 -0.29
#> worldpop             0.13    -0.50    -0.53     -0.54      0.20      0.74 -0.25
#> livestock_ctl        0.28    -0.33    -0.38     -0.42      0.44      0.65 -0.21
#> livestock_shp        0.35    -0.34    -0.33     -0.36      0.45      0.62 -0.32
#> corridors            0.05     0.48     0.64      0.66     -0.17     -0.67  0.15
#> species_richness    -0.03    -0.68    -0.76     -0.77      0.12      0.90 -0.24
#>                   ndvi footprint worldpop livestock_ctl livestock_shp corridors
#> bio1_25m         -0.01      0.10     0.13          0.28          0.35      0.05
#> bio2_25m         -0.62     -0.48    -0.50         -0.33         -0.34      0.48
#> bio5_25m         -0.66     -0.52    -0.53         -0.38         -0.33      0.64
#> bio10_25m        -0.68     -0.55    -0.54         -0.42         -0.36      0.66
#> bio15_25m         0.05      0.19     0.20          0.44          0.45     -0.17
#> bio18_25m         0.85      0.69     0.74          0.65          0.62     -0.67
#> dem              -0.20     -0.29    -0.25         -0.21         -0.32      0.15
#> ndvi              1.00      0.72     0.78          0.68          0.64     -0.69
#> footprint         0.72      1.00     0.91          0.76          0.78     -0.67
#> worldpop          0.78      0.91     1.00          0.78          0.78     -0.69
#> livestock_ctl     0.68      0.76     0.78          1.00          0.89     -0.66
#> livestock_shp     0.64      0.78     0.78          0.89          1.00     -0.62
#> corridors        -0.69     -0.67    -0.69         -0.66         -0.62      1.00
#> species_richness  0.88      0.77     0.81          0.69          0.67     -0.76
#>                  species_richness
#> bio1_25m                    -0.03
#> bio2_25m                    -0.68
#> bio5_25m                    -0.76
#> bio10_25m                   -0.77
#> bio15_25m                    0.12
#> bio18_25m                    0.90
#> dem                         -0.24
#> ndvi                         0.88
#> footprint                    0.77
#> worldpop                     0.81
#> livestock_ctl                0.69
#> livestock_shp                0.67
#> corridors                   -0.76
#> species_richness             1.00

The diagonal is 1 because each variable is compared with itself. The two halves repeat the same pairs. In the heatmaps, blue means positive and red means negative; stronger colour indicates a stronger relationship. Read the numbers as well as the colours. The same −1 to +1 colour scale is used for both figures.

draw_correlation <- function(m, title, filename) {
  png(filename, width = 1700, height = 1550, res = 150)
  on.exit(dev.off())
  par(mar = c(9, 9, 4, 2))
  n <- ncol(m)
  colours <- colorRampPalette(c("#b2182b", "white", "#2166ac"))(101)
  image(seq_len(n), seq_len(n), t(m[n:1, ]), col = colours,
        zlim = c(-1, 1), axes = FALSE, xlab = "", ylab = "", main = title)
  axis(1, at = seq_len(n), labels = colnames(m), las = 2, cex.axis = 0.8)
  axis(2, at = seq_len(n), labels = rev(rownames(m)), las = 2, cex.axis = 0.8)
  for (row in seq_len(n)) {
    for (column in seq_len(n)) {
      value <- m[row, column]
      text(column, n + 1 - row, sprintf("%.2f", value), cex = 0.72,
           col = if (abs(value) > 0.6) "white" else "black")
    }
  }
  mtext("Red: negative | White: near zero | Blue: positive", side = 3, line = 0.5, cex = 0.85)
}

draw_correlation(spearman, "Spearman: similarity in covariate ranks",
                 file.path(output, "figures", "correlation_spearman.png"))
draw_correlation(pearson, "Pearson: linear covariate relationships",
                 file.path(output, "figures", "correlation_pearson.png"))

13 Identify pairs to review and inspect scatterplots

For this exercise, flag a pair when absolute Spearman or Pearson correlation is at least 0.70. This is a discussion threshold, not a universal rule or a statistical significance test. Strong negative relationships deserve review too. The full table retains all pairs so students can try another threshold.

indices <- which(upper.tri(spearman), arr.ind = TRUE)
pairs_table <- data.frame(
  variable_1 = rownames(spearman)[indices[, 1]],
  variable_2 = colnames(spearman)[indices[, 2]],
  spearman = spearman[indices], pearson = pearson[indices],
  n_locations = nrow(x_shared))
pairs_table$review <- abs(pairs_table$spearman) >= 0.70 |
                      abs(pairs_table$pearson) >= 0.70
pairs_table <- pairs_table[order(-abs(pairs_table$spearman)), ]
write.csv(pairs_table, file.path(output, "correlation_pairs.csv"), row.names = FALSE)
review_pairs <- pairs_table[pairs_table$review, ]
review_pairs[, c("variable_1", "variable_2", "spearman", "pearson")]
#>       variable_1       variable_2   spearman     pearson
#> 6       bio5_25m        bio10_25m  0.9813550  0.97447113
#> 45     footprint         worldpop  0.9059150  0.33250798
#> 84     bio18_25m species_richness  0.9016610  0.82101309
#> 66 livestock_ctl    livestock_shp  0.8926644  0.47574631
#> 86          ndvi species_richness  0.8837323  0.95043624
#> 27     bio18_25m             ndvi  0.8473499  0.80870274
#> 88      worldpop species_richness  0.8075341  0.15199729
#> 65      worldpop    livestock_shp  0.7809999  0.10988028
#> 44          ndvi         worldpop  0.7778135  0.09077485
#> 55      worldpop    livestock_ctl  0.7771402  0.07768134
#> 64     footprint    livestock_shp  0.7762253  0.42597837
#> 87     footprint species_richness  0.7731940  0.72919172
#> 82     bio10_25m species_richness -0.7677001 -0.72972304
#> 3       bio2_25m         bio5_25m  0.7639146  0.82662212
#> 54     footprint    livestock_ctl  0.7588123  0.46113032
#> 81      bio5_25m species_richness -0.7572242 -0.72242954
#> 91     corridors species_richness -0.7569653 -0.62369834
#> 42     bio18_25m         worldpop  0.7372348  0.12572078
#> 36          ndvi        footprint  0.7233894  0.65360645
#> 5       bio2_25m        bio10_25m  0.6939029  0.71975667
#> 14     bio10_25m        bio18_25m -0.6865441 -0.73282504
#> 80      bio2_25m species_richness -0.6791437 -0.71183529
#> 25     bio10_25m             ndvi -0.6779258 -0.72419499
#> 13      bio5_25m        bio18_25m -0.6729763 -0.75022031
#> 24      bio5_25m             ndvi -0.6561419 -0.70836951
#> 12      bio2_25m        bio18_25m -0.6222245 -0.71845101

Inspect raw scatterplots and rank scatterplots for the three strongest Spearman pairs. A straight diagonal in the ranks indicates a strong monotonic relationship. If Pearson and Spearman differ, look for curvature, extremes or many tied values. A strong relationship may reflect a broad regional gradient rather than a direct connection between the two variables.

selected_pairs <- head(pairs_table, 3)
set.seed(2026)
display_rows <- sample.int(nrow(x_shared), min(3000, nrow(x_shared)))
png(file.path(output, "figures", "correlation_scatterplots.png"),
    width = 1700, height = 1800, res = 150)
par(mfrow = c(3, 2), mar = c(4, 4, 3, 1))
for (i in seq_len(nrow(selected_pairs))) {
  a <- selected_pairs$variable_1[i]
  b <- selected_pairs$variable_2[i]
  plot(x_shared[[a]][display_rows], x_shared[[b]][display_rows],
       pch = 16, cex = 0.35, col = adjustcolor("steelblue", alpha.f = 0.3),
       xlab = a, ylab = b,
       main = paste("Raw values | Pearson", round(selected_pairs$pearson[i], 2)))
  plot(rank(x_shared[[a]])[display_rows], rank(x_shared[[b]])[display_rows],
       pch = 16, cex = 0.35, col = adjustcolor("darkorange3", alpha.f = 0.3),
       xlab = paste("Rank of", a), ylab = paste("Rank of", b),
       main = paste("Ranks | Spearman", round(selected_pairs$spearman[i], 2)))
}
dev.off()
#> quartz_off_screen 
#>                 2

The plots display up to 3,000 shared rows for readability. Coefficients and ranks are calculated using all complete shared locations, not this plotting subset.

Worked reading of the supplied data. BIO5 and BIO10 have Spearman correlation about 0.98 and Pearson correlation about 0.97: their spatial patterns strongly overlap under both measures. Footprint and population density have Spearman about 0.91 but Pearson about 0.33: their ranks align much more closely than their raw values follow a straight line. The two scatterplot views help explain the difference. These examples prompt a redundancy review; neither automatically determines which criterion to retain.

14 Use the results to discuss the MCDA criteria

Correlation is evidence about redundancy, not a ranking of criterion importance. Giving separate weights to several very similar layers can give their shared pattern more influence in a weighted sum. Conversely, a weakly correlated variable can be irrelevant or poorly measured. No variable is automatically retained or deleted by this tutorial.

Use this sequence in the classroom:

1. Define the decision. Write one sentence explaining why each candidate criterion belongs in the MCDA.
2. Check the data. Examine units, dates, missingness, map patterns and outliers. A high-quality layer is preferable to a questionable proxy for the same concept.
3. Review flagged pairs. Read both correlation coefficients, their maps and their scatterplots. Decide whether they represent the same concept or distinct criteria that happen to vary together.
4. Make a reasoned choice. If two layers duplicate the same concept, consider retaining the clearer or better measured one. If both are needed, explain why and consider their combined influence when assigning weights.
5. Check the chosen set again. Filtering is not transitive: if A correlates with B and B with C, A need not correlate strongly with C. Review all retained pairs.

Correlations depend on the study boundary, native raster resolution and sampling locations. As an extension, change grid_spacing_m to 40000 and compare results in a different output folder. A relationship that changes substantially deserves closer inspection. Point spacing does not remove spatial autocorrelation.

These calculations use raw stored covariates. Later MCDA transformations, such as thresholds and saturation, can change their relationships; raw correlation is an initial screen, not a complete assessment of overlap in the final score layers.

worksheet <- data.frame(covariate = names(paths), description = unname(labels),
  reason_for_inclusion = "", data_quality_notes = "", correlated_partners = "",
  proposed_decision = "", justification = "")
for (i in seq_len(nrow(worksheet))) {
  name <- worksheet$covariate[i]
  partners <- c(review_pairs$variable_2[review_pairs$variable_1 == name],
                review_pairs$variable_1[review_pairs$variable_2 == name])
  worksheet$correlated_partners[i] <- paste(partners, collapse = "; ")
}
write.csv(worksheet, file.path(output, "covariate_selection_worksheet.csv"),
          row.names = FALSE)

15 Check the calculations

These checks confirm that the requested variables were processed, the sample is complete, the matrices are symmetric and the reported variance equals SD squared. They check arithmetic and output completeness, not the suitability of a final MCDA.

stopifnot(nrow(case_summary) == 14, nrow(regional_summary) == 14,
          all(case_summary$n + case_summary$missing == nrow(cases)),
          isTRUE(all.equal(case_summary$variance, case_summary$sd^2, tolerance = 1e-12)),
          all(file.exists(plot_manifest$map)), all(file.exists(plot_manifest$distribution)),
          !anyNA(x_shared), ncol(x_shared) == 14,
          nrow(pairs_table) == choose(14, 2),
          isTRUE(all.equal(spearman, t(spearman))),
          isTRUE(all.equal(pearson, t(pearson))),
          all(abs(spearman) <= 1 + 1e-12), all(abs(pearson) <= 1 + 1e-12),
          all(abs(diag(spearman) - 1) < 1e-12))
capture.output(sessionInfo(), file = file.path(output, "sessionInfo.txt"))
message("Checks passed: 14 covariates, 28 map/distribution figures and 91 correlation pairs.")
#> Checks passed: 14 covariates, 28 map/distribution figures and 91 correlation pairs.

Practice for the QGIS course

1. Put the footprint map and its distribution figure next to each other in a QGIS Print Layout. Explain the mean–median difference.
2. Compare regional and case-location distributions for two climate variables. Describe the patterns without interpreting them as causal.
3. Explain why comparing temperature and population density using raw variance or CV can be misleading.
4. Find a strongly positive and a strongly negative correlation. Locate the corresponding patterns on the maps.
5. Find a pair for which Pearson and Spearman differ. Inspect its scatterplot and suggest a reason.
6. Complete the selection worksheet. Propose a small set of criteria and justify every inclusion or exclusion using relevance, quality and redundancy.
7. Explain why shared spatial locations are necessary, and why a large number of sampled locations does not mean independent observations.

Files produced and running the tutorial

Open the HTML for the worked results. Run Covariate_Summary_Tutorial_V2.R from the project root or one of its subfolders to reproduce them. The .Rmd is the editable tutorial source; it can be rendered in RStudio with knitr, rmarkdown and Pandoc installed. The runnable .R requires only terra and base R.

Results are saved in @covariate-summary/V2/outputs, separately from other tutorials:

• regional_summary.csv: exact native-cell summaries inside country polygons.
• case_summary.csv and case_values.csv: descriptive case-location results.
• plot_manifest.csv and figures/: 14 maps, 14 density/boxplots, the log-view example, the sampling map, two correlation heatmaps and scatterplots.
• correlation_locations.csv: shared point coordinates, extracted values and inclusion flag.
• correlation_coverage.csv and correlation_sampling_summary.csv: missingness and retained sample size.
• correlation_pearson.csv, correlation_spearman.csv and correlation_pairs.csv: coefficients and pairs to review.
• covariate_selection_worksheet.csv: a discussion worksheet; decisions are left to students.
• sessionInfo.txt: the R environment used.

All tables and figures are computed from the supplied inputs. Only the descriptive statistics and correlation methods explained here are required to use the tutorial.