Spatial Data in the mlr3 Ecosystem

Run a land cover classification of the city of Leipzig.

Author

Marc Becker

Published

February 27, 2023

Scope

Working with spatial data in R requires a lot of data wrangling e.g. reading from different file formats, converting between spatial formats, creating tables from point layers, and predicting spatial raster images. The goal of mlr3spatial is to simplify these workflows within the mlr3 ecosystem. As a practical example, we will perform a land cover classification for the city of Leipzig, Germany. Figure 1 illustrates the typical workflow for this type of task: Load the training data, create a spatial task, train a learner with it, and predict the final raster image.

%%{ init: { 'flowchart': { 'curve': 'bump' } } }%%

flowchart LR
    subgraph files[Files]
    vector[Vector]
    raster[Raster]
    end
    subgraph load[Load Data]
    sf
    terra
    end
    vector --> sf
    raster --> terra
    subgraph train_model[Train Model]
    task[Task]
    learner[Learner]
    end
    terra --> prediction_raster
    task --> learner
    sf --> task
    subgraph predict[Spatial Prediction]
    prediction_raster[Raster Image]
    end
    learner --> prediction_raster
Figure 1: Spatial prediction workflow in mlr3spatial.

We assume that you are familiar with the mlr3 ecosystem and know the basic concepts of remote sensing. If not, we recommend reading the mlr3book first. If you are interested in spatial resampling, check out the book chapter on spatial analysis.

Land Cover Classification

Land cover is the physical material or vegetation that covers the surface of the earth, including both natural and human-made features. Understanding land cover patterns and changes over time is critical for addressing global environmental challenges, such as climate change, land degradation, and loss of biodiversity. Land cover classification is the process of assigning land cover classes to pixels in a raster image. With mlr3spatial, we can easily perform a land cover classification within the mlr3 ecosystem.

Before we can start the land cover classification, we need to load the necessary packages. The mlr3spatial package relies on terra for processing raster data and sf for vector data. These widely used packages read all common raster and vector formats. Additionally, the stars and raster packages are supported.

library(mlr3)
library(mlr3spatial)
library(terra, exclude = "resample")
library(sf)

We will work with a Sentinel-2 scene of the city of Leipzig which consists of 7 bands with a 10 or 20m spatial resolution and an NDVI band. The data is included in the mlr3spatial package. We use the terra::rast() to load the TIFF raster file.

leipzig_raster = rast(system.file("extdata", "leipzig_raster.tif", package = "mlr3spatial"))
leipzig_raster
class       : SpatRaster 
dimensions  : 206, 154, 8  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : 731810, 733350, 5692030, 5694090  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 32N (EPSG:32632) 
source      : leipzig_raster.tif 
names       :  b02,  b03,  b04,  b06,  b07,  b08, ... 
min values  :  846,  645,  366,  375,  401,  374, ... 
max values  : 4705, 4880, 5451, 4330, 5162, 5749, ...

The training data is a GeoPackage point layer with land cover labels and spectral features. We load the file and create a simple feature point layer.

leipzig_vector = read_sf(system.file("extdata", "leipzig_points.gpkg", package = "mlr3spatial"), stringsAsFactors = TRUE)
leipzig_vector
Simple feature collection with 97 features and 9 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 731930.5 ymin: 5692136 xmax: 733220.3 ymax: 5693968
Projected CRS: WGS 84 / UTM zone 32N
# A tibble: 97 × 10
     b02   b03   b04   b06   b07   b08   b11    ndvi land_cover               geom
   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl> <fct>             <POINT [m]>
 1   903   772   426  2998  4240  4029  1816  0.809  forest     (732480.1 5693957)
 2  1270  1256  1081  1998  2493  2957  2073  0.465  urban      (732217.4 5692769)
 3  1033   996   777  2117  2748  2799  1595  0.565  urban      (732737.2 5692469)
 4   962   773   500   465   505   396   153 -0.116  water      (733169.3 5692777)
 5  1576  1527  1626  1715  1745  1768  1980  0.0418 urban      (732202.2 5692644)
 6  1125  1185   920  3058  3818  3758  2682  0.607  pasture      (732153 5693059)
 7   880   746   424  2502  3500  3397  1469  0.778  forest     (731937.9 5693722)
 8  1332  1251  1385  1663  1799  1640  1910  0.0843 urban      (732416.2 5692324)
 9   940   741   475   452   515   400   139 -0.0857 water      (732933.7 5693344)
10   902   802   454  2764  3821  3666  1567  0.780  forest     (732411.3 5693352)
# ℹ 87 more rows

We plot both layers to get an overview of the data. The training points are located in the districts of Lindenau and Zentrum West.

Code 
library(ggplot2)
library(tidyterra, exclude = "filter")

ggplot() +
  geom_spatraster_rgb(data = leipzig_raster, r = 3, g = 2, b = 1, max_col_value = 5451) +
  geom_spatvector(data = leipzig_vector, aes(color = land_cover)) +
  scale_color_viridis_d(name = "Land cover", labels = c("Forest", "Pastures", "Urban", "Water")) +
  theme_minimal()

The as_task_classif_st() function directly creates a spatial task from the point layer. This makes it unnecessary to transform the point layer to a data.frame with coordinates. Spatial tasks additionally store the coordinates of the training points. The coordinates are useful when estimating the performance of the model with spatial resampling.

task = as_task_classif_st(leipzig_vector, target = "land_cover")
task
<TaskClassifST:leipzig_vector> (97 x 9)
* Target: land_cover
* Properties: multiclass
* Features (8):
  - dbl (8): b02, b03, b04, b06, b07, b08, b11, ndvi
* Coordinates:
           X       Y
       <num>   <num>
 1: 732480.1 5693957
 2: 732217.4 5692769
 3: 732737.2 5692469
 4: 733169.3 5692777
 5: 732202.2 5692644
---                 
93: 733018.7 5692342
94: 732551.4 5692887
95: 732520.4 5692589
96: 732542.2 5692204
97: 732437.8 5692300

Now we can train a model with the task. We use a simple decision tree learner from the rpart package. The "classif_st" task is a specialization of the "classif" task and therefore works with all "classif" learners.

learner = lrn("classif.rpart")
learner$train(task)

To get a complete land cover classification of Leipzig, we have to predict on each pixel and return a raster image with these predictions. The $predict() method of the learner only works for tabular data. To predict a raster image, we use the predict_spatial() function.

# predict land cover map
land_cover = predict_spatial(leipzig_raster, learner)
Warning in warn_deprecated("DataBackend$data_formats"): DataBackend$data_formats is deprecated and will be removed in
the future.
Code 
ggplot() +
  geom_spatraster(data = land_cover) +
  scale_fill_viridis_d(name = "Land cover", labels = c("Forest", "Pastures", "Urban", "Water")) +
  theme_minimal()

Conclusion

Working with spatial data in R is very easy with the mlr3spatial package. You can quickly train a model with a point layer and predict a raster image. The mlr3spatial package is still in development and we are looking forward to your feedback and contributions.

Session Information

sessioninfo::session_info(info = "packages")
═ Session info ═══════════════════════════════════════════════════════════════════════════════════════════════════════
─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────
 ! package        * version    date (UTC) lib source
   backports        1.5.0      2024-05-23 [1] CRAN (R 4.4.1)
   checkmate        2.3.2      2024-07-29 [1] CRAN (R 4.4.1)
 P class            7.3-22     2023-05-03 [?] CRAN (R 4.4.0)
   classInt         0.4-10     2023-09-05 [1] CRAN (R 4.4.1)
   cli              3.6.3      2024-06-21 [1] CRAN (R 4.4.1)
 P codetools        0.2-20     2024-03-31 [?] CRAN (R 4.4.0)
   colorspace       2.1-1      2024-07-26 [1] CRAN (R 4.4.1)
   crayon           1.5.3      2024-06-20 [1] CRAN (R 4.4.1)
   data.table     * 1.16.2     2024-10-10 [1] CRAN (R 4.4.1)
   DBI              1.2.3      2024-06-02 [1] CRAN (R 4.4.1)
   digest           0.6.37     2024-08-19 [1] CRAN (R 4.4.1)
   dplyr            1.1.4      2023-11-17 [1] CRAN (R 4.4.1)
   e1071            1.7-14     2023-12-06 [1] CRAN (R 4.4.1)
   evaluate         0.24.0     2024-06-10 [1] CRAN (R 4.4.1)
   fansi            1.0.6      2023-12-08 [1] CRAN (R 4.4.1)
   farver           2.1.2      2024-05-13 [1] CRAN (R 4.4.1)
   fastmap          1.2.0      2024-05-15 [1] CRAN (R 4.4.1)
   future           1.34.0     2024-07-29 [1] CRAN (R 4.4.1)
   generics         0.1.3      2022-07-05 [1] CRAN (R 4.4.1)
   ggplot2        * 3.5.1      2024-04-23 [1] CRAN (R 4.4.1)
   globals          0.16.3     2024-03-08 [1] CRAN (R 4.4.1)
   glue             1.7.0      2024-01-09 [1] CRAN (R 4.4.1)
   gtable           0.3.5      2024-04-22 [1] CRAN (R 4.4.1)
   htmltools        0.5.8.1    2024-04-04 [1] CRAN (R 4.4.1)
   htmlwidgets      1.6.4      2023-12-06 [1] CRAN (R 4.4.1)
   jsonlite         1.8.9      2024-09-20 [1] CRAN (R 4.4.1)
 P KernSmooth       2.23-24    2024-05-17 [?] CRAN (R 4.4.0)
   knitr            1.48       2024-07-07 [1] CRAN (R 4.4.1)
   lgr              0.4.4      2022-09-05 [1] CRAN (R 4.4.1)
   lifecycle        1.0.4      2023-11-07 [1] CRAN (R 4.4.1)
   listenv          0.9.1      2024-01-29 [1] CRAN (R 4.4.1)
   magrittr         2.0.3      2022-03-30 [1] CRAN (R 4.4.1)
   mlr3           * 0.21.0     2024-09-24 [1] CRAN (R 4.4.1)
   mlr3misc         0.15.1     2024-06-24 [1] CRAN (R 4.4.1)
   mlr3spatial    * 0.5.0      2024-03-09 [1] CRAN (R 4.4.1)
   mlr3website    * 0.0.0.9000 2024-10-18 [1] Github (mlr-org/mlr3website@20d1ddf)
   munsell          0.5.1      2024-04-01 [1] CRAN (R 4.4.1)
   palmerpenguins   0.1.1      2022-08-15 [1] CRAN (R 4.4.1)
   paradox          1.0.1      2024-07-09 [1] CRAN (R 4.4.1)
   parallelly       1.38.0     2024-07-27 [1] CRAN (R 4.4.1)
   pillar           1.9.0      2023-03-22 [1] CRAN (R 4.4.1)
   pkgconfig        2.0.3      2019-09-22 [1] CRAN (R 4.4.1)
   proxy            0.4-27     2022-06-09 [1] CRAN (R 4.4.1)
   purrr            1.0.2      2023-08-10 [1] CRAN (R 4.4.1)
   R6               2.5.1      2021-08-19 [1] CRAN (R 4.4.1)
   Rcpp             1.0.13     2024-07-17 [1] CRAN (R 4.4.1)
   renv             1.0.3      2023-09-19 [1] CRAN (R 4.4.1)
   rlang            1.1.4      2024-06-04 [1] CRAN (R 4.4.1)
   rmarkdown        2.28       2024-08-17 [1] CRAN (R 4.4.1)
 P rpart            4.1.23     2023-12-05 [?] CRAN (R 4.4.0)
   scales           1.3.0      2023-11-28 [1] CRAN (R 4.4.1)
   sessioninfo      1.2.2      2021-12-06 [1] CRAN (R 4.4.1)
   sf             * 1.0-18     2024-10-11 [1] CRAN (R 4.4.1)
   stringi          1.8.4      2024-05-06 [1] CRAN (R 4.4.1)
   terra          * 1.7-83     2024-10-14 [1] CRAN (R 4.4.1)
   tibble           3.2.1      2023-03-20 [1] CRAN (R 4.4.1)
   tidyr            1.3.1      2024-01-24 [1] CRAN (R 4.4.1)
   tidyselect       1.2.1      2024-03-11 [1] CRAN (R 4.4.1)
   tidyterra      * 0.6.1      2024-06-08 [1] CRAN (R 4.4.1)
   units            0.8-5      2023-11-28 [1] CRAN (R 4.4.1)
   utf8             1.2.4      2023-10-22 [1] CRAN (R 4.4.1)
   uuid             1.2-1      2024-07-29 [1] CRAN (R 4.4.1)
   vctrs            0.6.5      2023-12-01 [1] CRAN (R 4.4.1)
   viridisLite      0.4.2      2023-05-02 [1] CRAN (R 4.4.1)
   withr            3.0.1      2024-07-31 [1] CRAN (R 4.4.1)
   xfun             0.48       2024-10-03 [1] CRAN (R 4.4.1)
   yaml             2.3.10     2024-07-26 [1] CRAN (R 4.4.1)

 [1] /home/marc/repositories/mlr3website/mlr-org/renv/library/R-4.4/x86_64-pc-linux-gnu
 [2] /home/marc/.cache/R/renv/sandbox/R-4.4/x86_64-pc-linux-gnu/9a444a72

 P ── Loaded and on-disk path mismatch.

──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────