diff --git a/.Rbuildignore b/.Rbuildignore
index 6187d6e..b0ca1b0 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -10,4 +10,5 @@ extra/*
 ^docs$
 ^pkgdown$
 ^CRAN-SUBMISSION$
-revdep
\ No newline at end of file
+revdep
+^\.gitattributes$
diff --git a/.gitattributes b/.gitattributes
new file mode 100644
index 0000000..43ab72a
--- /dev/null
+++ b/.gitattributes
@@ -0,0 +1,17 @@
+# Normalise line endings: LF in the repository, native on disk.
+# Without this, files edited on Windows show up as entirely rewritten and
+# diffs become unreadable.
+* text=auto eol=lf
+
+# Files git must not treat as text
+*.png binary
+*.jpg binary
+*.jpeg binary
+*.pdf binary
+*.rda binary
+*.RData binary
+*.rds binary
+*.gpkg binary
+*.shp binary
+*.shx binary
+*.dbf binary
diff --git a/DESCRIPTION b/DESCRIPTION
index 59322b1..5729608 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -21,7 +21,6 @@ Description: Precision agriculture spatial data
 License: MIT + file LICENSE
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 7.3.3
 Imports: 
     data.table,
     e1071,
@@ -43,3 +42,4 @@ Suggests:
 URL: https://ppaccioretti.github.io/paar/, https://github.com/PPaccioretti/paar
 VignetteBuilder: knitr, rmarkdown
 BugReports: https://github.com/PPaccioretti/paar/issues
+Config/roxygen2/version: 8.0.0
diff --git a/NAMESPACE b/NAMESPACE
index ca9f941..3fc8680 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,11 +1,14 @@
 # Generated by roxygen2: do not edit by hand
 
 S3method(cbind,paar)
+S3method(print,multispati)
 S3method(print,paar)
 S3method(print,summary.paar)
+S3method(summary,multispati)
 S3method(summary,paar)
 export(compare_zone)
 export(depurate)
 export(fuzzy_k_means)
 export(kmspc)
 export(spatial_t_test)
+export(spca_biplot)
diff --git a/NEWS.md b/NEWS.md
index f7963df..fb9aa86 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,29 @@
 # paar (development version)
 
+* Fixed: in `spca_summary_results`, the `Prop` column of `kmspc()` was a
+  proportion between 0 and 1 while `Acum. Prop.` in the same table was a
+  percentage between 0 and 100. `spca_biplot()` formats `Prop` with a per-cent
+  sign, so the axis labels reported one hundredth of the explained variability,
+  for example `CS1 (0.5%)` where the axis actually explains 50.6%. Both columns
+  are now percentages. The denominator is the total inertia of the analysis
+  rather than the number of retained eigenvalues, which differ when the data
+  matrix is rank deficient. The selection of components from
+  `explainedVariance` is unaffected, as it already operated on the cumulative
+  percentage. The `Prop` and `Acum. Prop.` columns of the classical PCA summary
+  returned by `only_spca_results = FALSE` were also proportions and are now
+  percentages, so that both tables can be read side by side.
+* Fixed: `kmspc()` computed the spatial PCA results but discarded them unless
+  `all_results = TRUE`. With the default arguments `pca_results` came back
+  empty, even though `only_spca_results = TRUE` is meant to return precisely
+  those results. The spatial PCA summary and the loadings are now always
+  returned.
+* `kmspc()` also returns `scores`, the site coordinates on the retained
+  spatial components. These are the values the clustering step operates on and
+  were previously computed but not exposed.
+* New `spca_biplot()`, which draws a biplot of the spatial PCA from a
+  `kmspc()` result, optionally colouring the sites by one of the computed
+  classifications.
+
 # paar 1.0.2
 
 * Help pages were revised for clarity, consistency, and improved user guidance [#21](https://github.com/PPaccioretti/paar/issues/21)
diff --git a/R/kmspc.R b/R/kmspc.R
index 7a34bcc..e64871a 100644
--- a/R/kmspc.R
+++ b/R/kmspc.R
@@ -67,9 +67,21 @@
 #'   \item{cluster}{\code{data.frame} with cluster assignments for each evaluated number of clusters}
 #'   \item{indices}{\code{data.frame} with clustering validity indices}
 #'   \item{summaryResults}{\code{data.frame} with clustering metrics (iterations, SSDW)}
-#'   \item{pca_results}{(optional) PCA and/or spatial PCA summaries depending on arguments}
+#'   \item{pca_results}{a list with the spatial PCA results, always present:
+#'     \code{spca_summary_results} (eigenvalues, spatial variance, explained
+#'     percentage of the total inertia, cumulative percentage and Moran index
+#'     per axis), \code{eigenvectors_used}
+#'     (variable loadings on the retained axes) and \code{scores} (site
+#'     coordinates on the retained axes, the input of the clustering step).
+#'     If \code{only_spca_results = FALSE} it also contains
+#'     \code{pca_results} with the classical PCA summary, and if
+#'     \code{all_results = TRUE}, the complete \code{pca_results_all} and
+#'     \code{spca_results_all} objects.}
 #' }
 #'
+#' @seealso [spca_biplot()] to display `eigenvectors_used` and `scores`
+#'   together in a biplot.
+#'
 #' @export
 #' @example inst/examples/kmspc.R
 #'
@@ -159,7 +171,9 @@ kmspc <- function(
 
   if (!only_spca_results) {
     autov_pca <- pca$eig
-    propvar_pca <- autov_pca / sum(autov_pca)
+    # Percentages, on the same scale as the spatial PCA table below, so that
+    # both summaries can be read side by side.
+    propvar_pca <- 100 * autov_pca / sum(autov_pca)
     propvaracum_pca <- cumsum(propvar_pca)
 
     my_pca <- lapply(pca$l1, spdep::moran.mc, lw, 999)
@@ -168,7 +182,7 @@ kmspc <- function(
 
     nfila_pca <- length(variables)
     eje_pca <- seq_len(nfila_pca)
-    resultado_pca = data.frame(
+    resultado_pca <- data.frame(
       eje_pca,
       autov_pca,
       propvar_pca,
@@ -176,7 +190,7 @@ kmspc <- function(
       my_pca
     )
     resultado_pca$eje_pca <- as.factor(resultado_pca$eje_pca)
-    names(resultado_pca) = c(
+    names(resultado_pca) <- c(
       "Axis",
       "Eigenvalue",
       "Prop",
@@ -190,8 +204,14 @@ kmspc <- function(
   invisible(utils::capture.output(resms <- summary(ms)))
   var_ms <- resms[, 2, drop = F]
   nfila_ms <- length(ms$eig)
-  propvar_ms <- var_ms / nfila_ms
-  propvaracum_ms <- cumsum(propvar_ms) * 100
+
+  # Both `Prop` and `Acum. Prop.` are expressed as percentages, on the same
+  # scale as the `explainedVariance` threshold. The denominator is the total
+  # inertia of the analysis, which equals the number of variables only when
+  # the data matrix has full rank.
+  total_inertia <- sum(pca$eig)
+  propvar_ms <- 100 * var_ms / total_inertia
+  propvaracum_ms <- cumsum(propvar_ms)
 
   eje_ms <- seq_len(nfila_ms)
   resultado_ms <-
@@ -221,22 +241,26 @@ kmspc <- function(
       unlist(propvaracum_ms)
     ))
 
+  data_clust <- ms$li[num_sPC]
+
+  if (inherits(data_clust, "sf")) {
+    data_clust <- sf::st_drop_geometry(data_clust)
+  }
+
   spca_results <- list(
     spca_summary_results = resultado_ms,
-    eigenvectors_used = ms$c1[, num_sPC, drop = FALSE]
+    eigenvectors_used = ms$c1[, num_sPC, drop = FALSE],
+    scores = data_clust
   )
 
   if (all_results) {
     spca_results <- append(spca_results, list(spca_results_all = ms))
-
-    pca_results <- append(pca_results, spca_results)
   }
 
-  data_clust <- ms$li[num_sPC]
-
-  if (inherits(data_clust, "sf")) {
-    data_clust <- sf::st_drop_geometry(data_clust)
-  }
+  # The spatial PCA results are always returned: they are the output the
+  # `only_spca_results` argument refers to, and `scores` together with
+  # `eigenvectors_used` are what `spca_biplot()` needs.
+  pca_results <- append(pca_results, spca_results)
 
   my_results <- make_clasification(
     data_clust,
@@ -355,7 +379,7 @@ summarize_indices <- function(indices, number_cluster) {
   }
 
   indicesresults <- data.frame(number_cluster, indices, IndN)
-  names(indicesresults) = c(
+  names(indicesresults) <- c(
     "Num. Cluster",
     "Xie Beni",
     # "Fukuyama Sugeno",
@@ -659,6 +683,7 @@ multispati <-
 
 
 #' @noRd
+#' @exportS3Method
 summary.multispati <- function(object, ...) {
   norm.w <- function(X, w) {
     f2 <- function(v) {
@@ -731,6 +756,7 @@ summary.multispati <- function(object, ...) {
 
 
 #' @noRd
+#' @exportS3Method
 print.multispati <- function(x, ...) {
   # cat("Multispati object \n")
   # cat("class: ")
diff --git a/R/spca_biplot.R b/R/spca_biplot.R
new file mode 100644
index 0000000..185f57b
--- /dev/null
+++ b/R/spca_biplot.R
@@ -0,0 +1,223 @@
+#' Biplot of a spatial PCA
+#'
+#' @description
+#' Builds a biplot from the spatial principal component analysis computed by
+#' [kmspc()]. Sites are drawn as points using their coordinates on the
+#' retained spatial components, and variables as arrows from the origin using
+#' their loadings.
+#'
+#' The biplot complements the clustering output: the angles between arrows
+#' approximate the correlations among the original variables, the length of an
+#' arrow indicates how well the plotted plane represents that variable, and the
+#' position of a site relative to an arrow indicates its relative value for
+#' that variable. When management zones are mapped onto the points, it also
+#' shows which variables drive the separation between zones.
+#'
+#' @param x a list returned by [kmspc()].
+#'
+#' @param axes \code{numeric} vector of length two with the spatial components
+#'   to be plotted.
+#'
+#' @param cluster optional. Either a \code{character} with the name of one of
+#'   the columns of \code{x$cluster} (for example \code{"Cluster_3"}), or a
+#'   vector with one value per site. Points are coloured according to it.
+#'
+#' @param arrow_scale \code{numeric}. Factor applied to the loadings so that the
+#'   arrows are legible next to the cloud of sites. If \code{NULL} (default) it
+#'   is computed so that the longest arrow spans about 85\% of the point cloud.
+#'   Because the same factor is applied to every arrow, it does not alter their
+#'   interpretation.
+#'
+#' @param point_size,point_alpha \code{numeric}. Size and opacity of the points.
+#'
+#' @param label_size \code{numeric}. Size of the variable labels.
+#'
+#' @details
+#' Note that the sign of a component is arbitrary: the eigendecomposition fixes
+#' the direction of each axis but not its orientation, so two runs may produce
+#' mirrored biplots without any change in the conclusions. Interpretations
+#' should therefore be phrased in terms of relative positions, that is, which
+#' variables oppose each other and which go together, rather than in terms of
+#' absolute signs.
+#'
+#' The correlations suggested by the angles between arrows are computed treating
+#' the sites as independent observations. With spatially autocorrelated layers
+#' that assumption does not hold, so the magnitudes remain a valid description
+#' of the observed association but their statistical significance is
+#' misleading. See [spatial_t_test()] for a test that corrects the degrees of
+#' freedom by the effective sample size.
+#'
+#' @return A \code{ggplot} object, which can be further modified with the usual
+#'   \code{ggplot2} functions.
+#'
+#' @seealso [kmspc()], [spatial_t_test()]
+#'
+#' @export
+#' @example inst/examples/spca_biplot.R
+#'
+spca_biplot <- function(
+  x,
+  axes = c(1, 2),
+  cluster = NULL,
+  arrow_scale = NULL,
+  point_size = 1,
+  point_alpha = 0.6,
+  label_size = 3.2
+) {
+  if (!requireNamespace("ggplot2", quietly = TRUE)) {
+    stop(
+      "Package 'ggplot2' is needed for spca_biplot(). ",
+      "Please install it with install.packages('ggplot2').",
+      call. = FALSE
+    )
+  }
+
+  spca <- x[["pca_results"]]
+
+  if (is.null(spca) || is.null(spca[["scores"]])) {
+    stop(
+      "'x' does not contain spatial PCA results. ",
+      "It must be the output of kmspc().",
+      call. = FALSE
+    )
+  }
+
+  if (length(axes) != 2 || !is.numeric(axes)) {
+    stop("'axes' must be a numeric vector of length two", call. = FALSE)
+  }
+
+  scores <- spca[["scores"]]
+  loadings <- spca[["eigenvectors_used"]]
+
+  if (max(axes) > ncol(scores)) {
+    stop(
+      "Only ",
+      ncol(scores),
+      " spatial component(s) were retained, so 'axes' cannot exceed that ",
+      "number. Increase 'explainedVariance' in kmspc() to retain more ",
+      "components.",
+      call. = FALSE
+    )
+  }
+
+  sites <- data.frame(
+    axis_x = scores[[axes[1]]],
+    axis_y = scores[[axes[2]]]
+  )
+
+  arrows <- data.frame(
+    arrow_x = loadings[[axes[1]]],
+    arrow_y = loadings[[axes[2]]],
+    variable = rownames(loadings)
+  )
+
+  # The arrows are rescaled so that they are legible next to the cloud of
+  # sites; a single common factor preserves their relative lengths and angles.
+  if (is.null(arrow_scale)) {
+    max_site <- max(abs(c(sites$axis_x, sites$axis_y)), na.rm = TRUE)
+    max_arrow <- max(abs(c(arrows$arrow_x, arrows$arrow_y)), na.rm = TRUE)
+    arrow_scale <- if (max_arrow > 0) 0.85 * max_site / max_arrow else 1
+  }
+
+  arrows$arrow_x <- arrows$arrow_x * arrow_scale
+  arrows$arrow_y <- arrows$arrow_y * arrow_scale
+
+  # Optional grouping of the sites
+  zone <- NULL
+  if (!is.null(cluster)) {
+    if (is.character(cluster) && length(cluster) == 1) {
+      available <- x[["cluster"]]
+      if (is.null(available)) {
+        stop("'x' does not contain cluster assignments", call. = FALSE)
+      }
+      available <- as.data.frame(available)
+      if (!cluster %in% names(available)) {
+        stop(
+          "'",
+          cluster,
+          "' is not one of the available clusterings: ",
+          paste(names(available), collapse = ", "),
+          call. = FALSE
+        )
+      }
+      zone <- available[[cluster]]
+    } else {
+      zone <- cluster
+    }
+
+    # kmspc() keeps the rows dropped for having missing values, so the cluster
+    # vector may be longer than the matrix of scores.
+    if (length(zone) != nrow(sites)) {
+      zone <- stats::na.omit(zone)
+    }
+    if (length(zone) != nrow(sites)) {
+      stop(
+        "'cluster' has ",
+        length(zone),
+        " values but there are ",
+        nrow(sites),
+        " sites",
+        call. = FALSE
+      )
+    }
+    sites$zone <- as.factor(zone)
+  }
+
+  explained <- spca[["spca_summary_results"]]
+  axis_label <- function(i) {
+    if (!is.null(explained) && "Prop" %in% names(explained)) {
+      sprintf("%s (%.1f%%)", names(scores)[i], explained[["Prop"]][i])
+    } else {
+      names(scores)[i]
+    }
+  }
+
+  out <- ggplot2::ggplot() +
+    ggplot2::geom_hline(yintercept = 0, colour = "grey80", linewidth = 0.3) +
+    ggplot2::geom_vline(xintercept = 0, colour = "grey80", linewidth = 0.3)
+
+  if (is.null(zone)) {
+    out <- out +
+      ggplot2::geom_point(
+        data = sites,
+        ggplot2::aes(x = axis_x, y = axis_y),
+        colour = "grey60",
+        size = point_size,
+        alpha = point_alpha
+      )
+  } else {
+    out <- out +
+      ggplot2::geom_point(
+        data = sites,
+        ggplot2::aes(x = axis_x, y = axis_y, colour = zone),
+        size = point_size,
+        alpha = point_alpha
+      ) +
+      ggplot2::labs(colour = "Zone")
+  }
+
+  out +
+    ggplot2::geom_segment(
+      data = arrows,
+      ggplot2::aes(x = 0, y = 0, xend = arrow_x, yend = arrow_y),
+      arrow = ggplot2::arrow(length = ggplot2::unit(0.18, "cm")),
+      colour = "grey15",
+      linewidth = 0.5
+    ) +
+    ggplot2::geom_text(
+      data = arrows,
+      ggplot2::aes(x = arrow_x * 1.12, y = arrow_y * 1.12, label = variable),
+      size = label_size,
+      fontface = "bold"
+    ) +
+    ggplot2::coord_equal() +
+    ggplot2::labs(x = axis_label(axes[1]), y = axis_label(axes[2])) +
+    ggplot2::theme_bw()
+}
+
+
+# Columns built inside spca_biplot() and referenced within ggplot2::aes().
+# Declared here so that R CMD check does not report them as undefined globals.
+utils::globalVariables(
+  c("axis_x", "axis_y", "arrow_x", "arrow_y", "variable", "zone")
+)
diff --git a/inst/examples/spca_biplot.R b/inst/examples/spca_biplot.R
new file mode 100644
index 0000000..abdc211
--- /dev/null
+++ b/inst/examples/spca_biplot.R
@@ -0,0 +1,15 @@
+library(sf)
+data(wheat, package = 'paar')
+
+# Transform the data.frame into a sf object
+wheat_sf <- st_as_sf(wheat, coords = c('x', 'y'), crs = 32720)
+
+# Run the kmspc function
+kmspc_results <- kmspc(wheat_sf, number_cluster = 2:4)
+
+# Biplot of the first two spatial components
+spca_biplot(kmspc_results)
+
+# Colouring the sites by the three-zone classification shows which variables
+# drive the separation between management zones
+spca_biplot(kmspc_results, cluster = "Cluster_3")
diff --git a/man/kmspc.Rd b/man/kmspc.Rd
index 81ef64d..0f305a5 100644
--- a/man/kmspc.Rd
+++ b/man/kmspc.Rd
@@ -63,7 +63,16 @@ A list with the following elements:
   \item{cluster}{\code{data.frame} with cluster assignments for each evaluated number of clusters}
   \item{indices}{\code{data.frame} with clustering validity indices}
   \item{summaryResults}{\code{data.frame} with clustering metrics (iterations, SSDW)}
-  \item{pca_results}{(optional) PCA and/or spatial PCA summaries depending on arguments}
+  \item{pca_results}{a list with the spatial PCA results, always present:
+    \code{spca_summary_results} (eigenvalues, spatial variance, explained
+    percentage of the total inertia, cumulative percentage and Moran index
+    per axis), \code{eigenvectors_used}
+    (variable loadings on the retained axes) and \code{scores} (site
+    coordinates on the retained axes, the input of the clustering step).
+    If \code{only_spca_results = FALSE} it also contains
+    \code{pca_results} with the classical PCA summary, and if
+    \code{all_results = TRUE}, the complete \code{pca_results_all} and
+    \code{spca_results_all} objects.}
 }
 }
 \description{
@@ -116,3 +125,7 @@ wheat_clustered <- cbind(wheat_sf, kmspc_results$cluster)
 # Plot the results
 plot(wheat_clustered[, "Cluster_2"])
 }
+\seealso{
+[spca_biplot()] to display `eigenvectors_used` and `scores`
+  together in a biplot.
+}
diff --git a/man/spca_biplot.Rd b/man/spca_biplot.Rd
new file mode 100644
index 0000000..8c1d01d
--- /dev/null
+++ b/man/spca_biplot.Rd
@@ -0,0 +1,88 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/spca_biplot.R
+\name{spca_biplot}
+\alias{spca_biplot}
+\title{Biplot of a spatial PCA}
+\usage{
+spca_biplot(
+  x,
+  axes = c(1, 2),
+  cluster = NULL,
+  arrow_scale = NULL,
+  point_size = 1,
+  point_alpha = 0.6,
+  label_size = 3.2
+)
+}
+\arguments{
+\item{x}{a list returned by [kmspc()].}
+
+\item{axes}{\code{numeric} vector of length two with the spatial components
+to be plotted.}
+
+\item{cluster}{optional. Either a \code{character} with the name of one of
+the columns of \code{x$cluster} (for example \code{"Cluster_3"}), or a
+vector with one value per site. Points are coloured according to it.}
+
+\item{arrow_scale}{\code{numeric}. Factor applied to the loadings so that the
+arrows are legible next to the cloud of sites. If \code{NULL} (default) it
+is computed so that the longest arrow spans about 85\% of the point cloud.
+Because the same factor is applied to every arrow, it does not alter their
+interpretation.}
+
+\item{point_size, point_alpha}{\code{numeric}. Size and opacity of the points.}
+
+\item{label_size}{\code{numeric}. Size of the variable labels.}
+}
+\value{
+A \code{ggplot} object, which can be further modified with the usual
+  \code{ggplot2} functions.
+}
+\description{
+Builds a biplot from the spatial principal component analysis computed by
+[kmspc()]. Sites are drawn as points using their coordinates on the
+retained spatial components, and variables as arrows from the origin using
+their loadings.
+
+The biplot complements the clustering output: the angles between arrows
+approximate the correlations among the original variables, the length of an
+arrow indicates how well the plotted plane represents that variable, and the
+position of a site relative to an arrow indicates its relative value for
+that variable. When management zones are mapped onto the points, it also
+shows which variables drive the separation between zones.
+}
+\details{
+Note that the sign of a component is arbitrary: the eigendecomposition fixes
+the direction of each axis but not its orientation, so two runs may produce
+mirrored biplots without any change in the conclusions. Interpretations
+should therefore be phrased in terms of relative positions, that is, which
+variables oppose each other and which go together, rather than in terms of
+absolute signs.
+
+The correlations suggested by the angles between arrows are computed treating
+the sites as independent observations. With spatially autocorrelated layers
+that assumption does not hold, so the magnitudes remain a valid description
+of the observed association but their statistical significance is
+misleading. See [spatial_t_test()] for a test that corrects the degrees of
+freedom by the effective sample size.
+}
+\examples{
+library(sf)
+data(wheat, package = 'paar')
+
+# Transform the data.frame into a sf object
+wheat_sf <- st_as_sf(wheat, coords = c('x', 'y'), crs = 32720)
+
+# Run the kmspc function
+kmspc_results <- kmspc(wheat_sf, number_cluster = 2:4)
+
+# Biplot of the first two spatial components
+spca_biplot(kmspc_results)
+
+# Colouring the sites by the three-zone classification shows which variables
+# drive the separation between management zones
+spca_biplot(kmspc_results, cluster = "Cluster_3")
+}
+\seealso{
+[kmspc()], [spatial_t_test()]
+}