Advanced Analysis with Swiss Data • mapAI

This vignette demonstrates an advanced workflow for analyzing and correcting a map with complex, non-linear distortions. While the first vignette showed a general correction, this one focuses on a more challenging real-world scenario.

We will use the swiss_cps dataset, which is included with the mapAI package. The key feature of this dataset is that the gross, linear distortions have already been removed with a Helmert transformation. The remaining errors are subtle, non-linear, and spatially variable—a perfect test case for the more advanced models in mapAI.

Our goal is to model these non-linear residuals and visualize the effect of the correction on a reference line grid.

# Load the necessary libraries for our workflow
library(mapAI)
library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.4.0; sf_use_s2() is TRUE
library(ggplot2)
library(magrittr) # For the %>% pipe operator

Step 1: Load and Explore the Data

We begin by loading the swiss_cps dataset. This sf object contains 343 control points where the source coordinates have been pre-aligned, so we can focus on the residual dx and dy errors.

# Load the built-in Swiss control points dataset
data(swiss_cps)

# Inspect the data structure
# Note that dx and dy are relatively small, representing the non-linear errors.
dplyr::glimpse(swiss_cps)
#> Rows: 343
#> Columns: 8
#> $ Index    <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18…
#> $ source_x <dbl> 612293.0, 612225.0, 615470.3, 616339.8, 616927.1, 617571.9, 6…
#> $ source_y <dbl> 267353.6, 270412.7, 270195.1, 269367.5, 269507.6, 269273.5, 2…
#> $ target_x <dbl> 611375.9, 611573.1, 615840.8, 616129.6, 617031.0, 618168.7, 6…
#> $ target_y <dbl> 267719.1, 270370.6, 270463.6, 269380.6, 268932.1, 269179.3, 2…
#> $ dx       <dbl> -917.1063, -651.9320, 370.4737, -210.1946, 103.8862, 596.7662…
#> $ dy       <dbl> 365.468249, -42.050863, 268.478200, 13.144895, -575.530782, -…
#> $ geometry <POINT [m]> POINT (612293 267353.6), POINT (612225 270412.7), POINT…

# Let's visualize the residual displacement vectors
plot_displacement(swiss_cps, title = "Residual Non-Linear Distortions")

The plot shows that even after a global alignment, there are clear spatial patterns in the remaining errors, which a simple linear model could not fix.

Step 2: Create a Reference Grid to Deform

To best visualize the distortion, we will create a regular grid of lines that covers the extent of our control points. This will serve as our “source map.” We will then apply our correction model to this grid to see how it is “un-warped.”

# Create a grid of polygons covering the data extent
grid_polygons <- sf::st_make_grid(swiss_cps, n = c(15, 15))

# Convert the polygon boundaries to a MULTILINESTRING object
# This gives us a clean line grid for visualization.
source_grid <- sf::st_cast(grid_polygons, "MULTILINESTRING") %>%
  sf::st_sf() # Convert to a full sf object

# Plot the original, regular grid
plot(sf::st_geometry(source_grid), main = "Original, Regular Reference Grid")

Step 3: Train a Non-Linear Model

Because the linear components of the distortion have already been removed from the swiss_cps data, using a helmert or lm model would be ineffective. This is the perfect scenario to use a Generalized Additive Model (gam), which is designed to capture exactly these kinds of smooth, non-linear spatial patterns.

# Train the GAM model to learn the residual distortion patterns
gam_model <- train_pai_model(swiss_cps, method = "gam")
#> Training 'gam' model...

Step 4: Apply the Correction to the Grid

Now, we use our trained gam_model to correct the regular line grid we created. The apply_pai_model function will transform every vertex in the grid according to the learned distortion.

# Apply the model to our source grid
corrected_grid <- apply_pai_model(gam_model, source_grid)
#> Applying PAI model to map features...
#> Correction complete.

Step 5: Visualize the Corrected Grid

This is the most intuitive part of the analysis. We overlay the corrected grid on the original grid. The correction will transform the regular, straight-lined grid to a warped one, in order to correct it.

# For easy plotting, add a 'status' column to each grid
source_grid$status <- "Original (Regular)"
corrected_grid$status <- "Corrected (Warped by model)"

# Combine them into a single sf object for ggplot2
comparison_data <- rbind(
  source_grid[, "status"],
  corrected_grid[, "status"]
)

ggplot(comparison_data) +
  geom_sf(aes(color = status, linetype = status), linewidth = 0.6) +
  scale_color_manual(
    name = "Grid Status",
    values = c("Original (Regular)" = "grey50", 
               "Corrected (Warped by model)" = "#e41a1c")
  ) +
  scale_linetype_manual(
    name = "Grid Status",
    values = c("Original (Regular)" = "dashed", 
               "Corrected (Warped by model)" = "solid")
  ) +
  labs(
    title = "Visualizing the Non-Linear Correction",
    subtitle = "The GAM model transforms the regular shape to warped grid"
  ) +
  theme_minimal()

Step 6: In-Depth Distortion Analysis

Finally, let’s use the advanced analysis functions to quantify the distortion that the model learned. We will analyze the distortion on a dense grid of points and visualize the results.

# 1. Create a dense grid of POINTS for analysis
analysis_points <- sf::st_make_grid(swiss_cps, n = c(25, 25)) %>%
  sf::st_centroid() %>%
  sf::st_sf()

# 2. Analyze the distortion using our trained GAM model
distortion_results <- analyze_distortion(gam_model, analysis_points)
#> Calculating distortion metrics for gam model...
#> Finalizing metrics from derivatives...
#> Distortion analysis complete.

# 3. Visualize the results for two key metrics
plot_shear <- plot_distortion_surface(
  distortion_results,
  metric = "max_shear",
  diverging = FALSE,
  palette = "magma"
) +
  labs(title = "Maximum Angular Distortion (°)")
#> Regular grid detected. Creating a surface plot with geom_raster().

plot_area <- plot_distortion_surface(
  distortion_results,
  metric = "log2_area_scale",
  diverging = TRUE
) +
  labs(title = "Log-Areal Distortion (log2)")
#> Regular grid detected. Creating a surface plot with geom_raster().

# Combine plots using the patchwork library
if (requireNamespace("patchwork", quietly = TRUE)) {
  plot_shear / plot_area
} else {
  plot_shear
}
#> Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
#> ℹ Consider using `geom_tile()` instead.
#> Raster pixels are placed at uneven horizontal intervals and will be shifted
#> ℹ Consider using `geom_tile()` instead.
#> Raster pixels are placed at uneven horizontal intervals and will be shifted
#> ℹ Consider using `geom_tile()` instead.
#> Raster pixels are placed at uneven horizontal intervals and will be shifted
#> ℹ Consider using `geom_tile()` instead.

This final analysis provides a quantitative map of the distortion, revealing complex patterns of shearing and scaling across the map extent that the gam model successfully identified and corrected.