Running R models from the GeoDMS

The GeoDMS can hand a data matrix to an R script, let R apply a fitted model to it, and read the result back. This is the route asked for in issue 659: instead of hard-coding coefficients in an Expression, such as

attribute<float32> suitability (domain) := logit(b1 * hoogte + b2 * afstandTot);

a model fitted in R - lm, glm, randomForest, xgboost, anything with a predict method - is saved with saveRDS and applied to the GeoDMS data.

Nothing R-specific is needed in the GeoDMS for this: the exchange runs over files and Exec_ec, exactly as it does for Python. Rscript simply takes the place of python.

Calling R in-process, through Rcpp, was explored and rejected. Rcpp - and R’s own C API - requires the toolchain that R itself was built with, which on Windows is the MinGW-w64 gcc/gfortran of Rtools, whereas the GeoDMS is built with MSVC. Bridging that would mean a separate gcc-built DLL behind a flat C ABI, and then hosting a single-threaded, non-reentrant interpreter with longjmp-based error handling inside an engine that schedules its work across worker threads. The file route costs one process start and one file round trip, and needs no such machinery.

the exchange format

Use Parquet. It carries the Value type of every column in the file itself, so no .csvt is needed and R gets correctly typed columns. On the R side either

  • arrow::read_parquet() / arrow::write_parquet(), or
  • nanoparquet::read_parquet() / nanoparquet::write_parquet(), which has no further dependencies and is therefore the lighter installation.

Arrow is the equally good sibling, written and read through the same GDAL storage managers, and R takes it with arrow::read_ipc_file() / arrow::write_ipc_file(). Pick Arrow when the same files are also read by Julia, which takes it directly with Arrow.Table - see Julia for a production pipeline built on it.

Csv also works and is a reasonable fallback, but then the types are guessed on at least one of the two sides.

the R script

The script takes the input file, the saved model and the output file as arguments, so that the GeoDMS decides where everything lives:

# predict.R  --  usage: Rscript predict.R <input.parquet> <model.rds> <output.parquet>
args  <- commandArgs(trailingOnly = TRUE)
stopifnot(length(args) == 3L)

X     <- nanoparquet::read_parquet(args[1])
model <- readRDS(args[2])

out <- data.frame(suitability = as.numeric(predict(model, newdata = X)))
nanoparquet::write_parquet(out, args[3])

Let the script fail loudly: an R stop() writes to stderr and makes Rscript exit non-zero, and since GeoDMS 20.16.0 both the message and the ExitCode reach the GeoDMS message log (see Exec_ec). On older versions a failing script under GeoDMS GUI leaves nothing but the ExitCode, so there it is worth appending > log.txt 2>&1 to the command and reading log.txt back as a str storage.

the GeoDMS configuration 

container suitability_model
{
   parameter<string> Dir      := Expand(., '%LocalDataProjDir%/R');
   parameter<string> RScript  := Expand(., '%env:ProgramFiles%/R/R-4.3.3/bin/Rscript.exe');

   // step 1: the predictors, written by the GeoDMS
   unit<uint32> X := domain
   ,  StorageName     = "= Dir + '/temp/X.parquet'"
   ,  StorageType     = "gdalwrite.vect"
   ,  StorageReadOnly = "False"
   {
      attribute<float32> hoogte      := domain/hoogte;
      attribute<float32> afstandTot  := domain/afstandTot;
   }

   // step 2: apply the saved model
   parameter<uint32> exitcode := exec_ec(RScript
      ,  DQuote(Dir + '/predict.R') + ' ' + DQuote(Dir + '/temp/X.parquet')
                                    + ' ' + DQuote(Dir + '/model.rds')
                                    + ' ' + DQuote(Dir + '/temp/y.parquet')
      ,  Dir);

   // step 3: read the predictions back. The ExitCode in the StorageName orders the read
   // after step 2; ':= domain' relates the imported rows to the configured domain, which
   // is safe here because the script preserves row order (see [Parquet](parquet.html) for the
   // alternatives when it does not).
   unit<uint32> y := domain
   ,  StorageName     = "= exitcode == 0 ? Dir + '/temp/y.parquet' : ''"
   ,  StorageType     = "gdal.vect"
   ,  StorageReadOnly = "True"
   {
      attribute<float32> suitability;
   }

   attribute<float32> suitability (domain) := y/suitability;
}

running it

Step 1 has to be a GeoDmsRun call of its own: a storage the GeoDMS writes cannot be ordered before an exec_ec inside one run, because exec_ec runs at meta-information time and the write happens later, when data is calculated. See Exec_ec for the mechanism and for the constructions that do not work. Steps 2 and 3 may share a call:

@echo off
set GEODMS="C:\Program Files\ObjectVision\GeoDms20.16.0\GeoDmsRun.exe"
set CFG=%~dp0cfg\main.dms

%GEODMS% %CFG% /suitability_model/X
if errorlevel 1 exit /b 1

%GEODMS% %CFG% /suitability_model/suitability

fitting the model

Fitting is a separate script over the same mechanism: write [y:X] from the GeoDMS, fit in R, saveRDS the model, and let the configuration above pick it up on the next run.

# fit.R  --  usage: Rscript fit.R <yX.parquet> <model.rds>
args <- commandArgs(trailingOnly = TRUE)
df   <- nanoparquet::read_parquet(args[1])
saveRDS(randomForest::randomForest(y ~ ., df), args[2])

notes

  • The GeoDMS does not install or manage R; the modeller installs it and configures the path to Rscript.exe. Keeping that path in one Parameter, as above, keeps the configuration portable.
  • R start-up plus package loading costs on the order of a second per call, so call the script once for the whole matrix, never per row or per region.
  • exec_ec occupies the main thread while R runs; the user interface does not repaint in the meantime. The exec: still waiting for ... after n seconds notices in the message log (20.16.0 and later) are what shows that the run is progressing.

Export an explicit id column with the predictors and join the results back on it, rather than relying on the script preserving row order; Julia shows that pattern with Rjoin.

see also