Spilhaus, as it flows downstream from PROJ: the R view

Also in this series: index | the Python notebook | the browser version | the other R runner | source | the story so far

The companion to the Python notebook. Same sources, same sections, same landmarks, different runtime. The point is the same too: +proj=spilhaus arrived in PROJ 9.6.0 (March 2025), and whether this R session can use it depends on which PROJ each package was linked against.

R’s version of the problem is not the Python one. On Linux, sf, terra, gdalraster and vapour all link one system libproj, so they normally agree; the question is what the system has, and on a stock Ubuntu that is often too old. On macOS and Windows, CRAN’s binaries bundle their own GDAL and PROJ, chosen by CRAN’s build recipes rather than by package authors. So the row that changes between machines is the operating system, not the wheel.

The stack here is deliberately low-level: gdalraster for I/O and the warp, wk for vertex-level geometry, ximage to draw a matrix. sf and terra appear only in the provenance table and in an appendix, as the same operation in one line, so you can see which PROJ each one silently used.

SPILHAUS <- "+proj=spilhaus"
SPILHAUS_EXPLICIT <- paste(
  "+proj=spilhaus +lon_0=66.94970198 +lat_0=-49.56371678",
  "+azi=40.17823482 +rot=45 +k_0=1 +x_0=0 +y_0=0 +R=6378137"
)
SPILHAUS_ALT <- "+proj=spilhaus +lon_0=-30 +lat_0=-20 +azi=20 +rot=45"

GEBCO_URL <- "/vsicurl/https://data.source.coop/alexgleith/gebco-2024/GEBCO_2024.tif"
RASTER_URL <- Sys.getenv("SPILHAUS_RASTER", GEBCO_URL)

NE_OCEAN_URL <- paste0(
  "/vsicurl/https://raw.githubusercontent.com/nvkelso/natural-earth-vector/",
  "master/110m_physical/ne_110m_ocean.shp"
)
OCEAN_URL <- Sys.getenv("SPILHAUS_OCEAN", NE_OCEAN_URL)

NPIX <- as.integer(Sys.getenv("SPILHAUS_NPIX", "1600"))

Sys.setenv(GDAL_DISABLE_READDIR_ON_OPEN = "EMPTY_DIR",
           GDAL_HTTP_MULTIRANGE = "YES",
           GDAL_HTTP_MERGE_CONSECUTIVE_RANGES = "YES")

1. Who has which PROJ?

One row per package. spilhaus is the result of asking that package to make a CRS from +proj=spilhaus.

HOBART <- cbind(147.33, -42.88)

## One forward-transform closure per package. The probe is "transform Hobart and get
## finite numbers", because several packages will happily build a CRS object from an
## unknown projection and only fail (or return NA) when asked to use it.
engine_forward <- list(
  gdalraster = function(lonlat, ps) gdalraster::transform_xy(lonlat, "OGC:CRS84", ps),
  PROJ       = function(lonlat, ps) PROJ::proj_trans(lonlat, ps, source = "OGC:CRS84"),
  sf         = function(lonlat, ps) sf::sf_project("OGC:CRS84", ps, lonlat, keep = TRUE),
  terra      = function(lonlat, ps) terra::crds(terra::project(terra::vect(lonlat, crs = "OGC:CRS84"), ps)),
  vapour     = function(lonlat, ps) vapour::vapour_project(lonlat, ps, "OGC:CRS84")
)
engine_inverse <- list(
  gdalraster = function(xy, ps) gdalraster::transform_xy(xy, ps, "OGC:CRS84"),
  PROJ       = function(xy, ps) PROJ::proj_trans(xy, "OGC:CRS84", source = ps),
  sf         = function(xy, ps) sf::sf_project(ps, "OGC:CRS84", xy, keep = TRUE),
  terra      = function(xy, ps) terra::crds(terra::project(terra::vect(xy, crs = ps), "OGC:CRS84")),
  vapour     = function(xy, ps) vapour::vapour_project(xy, "OGC:CRS84", ps)
)

probe <- function(pkg, ps = SPILHAUS) {
  out <- tryCatch({
    xy <- suppressWarnings(matrix(as.numeric(engine_forward[[pkg]](HOBART, ps)), ncol = 2))
    if (nrow(xy) != 1 || !all(is.finite(xy))) stop("transform returned NA/Inf") else "ok"
  }, error = function(e) paste("FAIL:", conditionMessage(e)))
  substr(out, 1, 80)
}
has <- function(pkg) requireNamespace(pkg, quietly = TRUE)
## "installed" is not "loadable": a binary built against another GDAL fails at load time
load_error <- function(pkg) tryCatch({ loadNamespace(pkg); NULL }, error = function(e) conditionMessage(e))
row_for_broken <- function(pkg, msg) data.frame(package = pkg, version = as.character(packageVersion(pkg)),
                                                gdal = "cannot load", proj = "cannot load",
                                                spilhaus = paste("FAIL:", substr(msg, 1, 60)))

rows <- list()

if (has("gdalraster") && !is.null(m <- load_error("gdalraster"))) {
  rows$gdalraster <- row_for_broken("gdalraster", m)
} else if (has("gdalraster")) {
  pv <- gdalraster::proj_version()
  rows$gdalraster <- data.frame(
    package = "gdalraster", version = as.character(packageVersion("gdalraster")),
    gdal = gdalraster::gdal_version()[4], proj = pv$name,
    spilhaus = probe("gdalraster"))
} else rows$gdalraster <- data.frame(package = "gdalraster", version = NA, gdal = NA, proj = NA, spilhaus = "not installed")

if (has("PROJ") && !is.null(m <- load_error("PROJ"))) {
  rows$PROJ <- row_for_broken("PROJ", m)
} else if (has("PROJ")) {
  rows$PROJ <- data.frame(
    package = "PROJ", version = as.character(packageVersion("PROJ")),
    gdal = "-", proj = PROJ::proj_version(),
    spilhaus = probe("PROJ"))
} else rows$PROJ <- data.frame(package = "PROJ", version = NA, gdal = "-", proj = NA, spilhaus = "not installed")

if (has("sf") && !is.null(m <- load_error("sf"))) {
  rows$sf <- row_for_broken("sf", m)
} else if (has("sf")) {
  v <- sf::sf_extSoftVersion()
  rows$sf <- data.frame(
    package = "sf", version = as.character(packageVersion("sf")),
    gdal = unname(v["GDAL"]), proj = unname(v["PROJ"]),
    spilhaus = probe("sf"))
} else rows$sf <- data.frame(package = "sf", version = NA, gdal = NA, proj = NA, spilhaus = "not installed")

if (has("terra") && !is.null(m <- load_error("terra"))) {
  rows$terra <- row_for_broken("terra", m)
} else if (has("terra")) {
  rows$terra <- data.frame(
    package = "terra", version = as.character(packageVersion("terra")),
    gdal = terra::gdal(lib = "gdal"), proj = terra::gdal(lib = "proj"),
    spilhaus = probe("terra"))
} else rows$terra <- data.frame(package = "terra", version = NA, gdal = NA, proj = NA, spilhaus = "not installed")

if (has("vapour") && !is.null(m <- load_error("vapour"))) {
  rows$vapour <- row_for_broken("vapour", m)
} else if (has("vapour")) {
  rows$vapour <- data.frame(
    package = "vapour", version = as.character(packageVersion("vapour")),
    gdal = vapour::vapour_gdal_version(), proj = vapour::vapour_proj_version(),
    spilhaus = probe("vapour"))
}

if (has("wk")) rows$wk <- data.frame(package = "wk", version = as.character(packageVersion("wk")),
                                     gdal = "-", proj = "- (no PROJ dependency)", spilhaus = "n/a")
if (has("ximage")) rows$ximage <- data.frame(package = "ximage", version = as.character(packageVersion("ximage")),
                                             gdal = "-", proj = "-", spilhaus = "n/a")

provenance <- do.call(rbind, rows); rownames(provenance) <- NULL
knitr::kable(provenance, caption = sprintf("R %s on %s", getRversion(), R.version$platform))

R 4.6.1 on x86_64-pc-linux-gnu

package	version	gdal	proj	spilhaus
gdalraster	2.6.1	3.8.4	9.4.0	FAIL: error importing SRS from user input
PROJ	0.7.0	-	9.4.0	FAIL: Invalid value for an argument
sf	1.1.3	3.8.4	9.4.0	FAIL: Invalid value for an argument: +proj=spilhaus
terra	1.9.50	3.8.4	9.4.0	FAIL: error in evaluating the argument ‘x’ in selecting a method for function ’c
wk	0.9.5	-	- (no PROJ dependency)	n/a

cli_version <- function(cmd, args = "--version") {
  exe <- Sys.which(cmd)
  if (!nzchar(exe)) return("not found")
  out <- tryCatch(system2(exe, args, stdout = TRUE, stderr = TRUE), error = function(e) conditionMessage(e))
  out[1]
}
cli <- c(gdalinfo = cli_version("gdalinfo"),
         gdal = cli_version("gdal"),
         proj = cli_version("proj", character()))  # bare `proj` prints "Rel. x.y.z" to stderr
knitr::kable(data.frame(tool = names(cli), version = unname(cli)), caption = "Command-line tools on PATH")

Command-line tools on PATH

tool	version
gdalinfo	GDAL 3.8.4, released 2024/02/08
gdal	not found
proj	Rel. 9.4.0, March 1st, 2024

capable <- provenance$package[provenance$spilhaus == "ok" & provenance$package %in% names(engine_forward)]
RASTER_OK <- has("gdalraster") && is.null(load_error("gdalraster"))
if (length(capable) == 0) {
  cat("No installed package can transform to +proj=spilhaus on this machine.",
      "The table above is the whole result; the document stops here.\n")
  knitr::knit_exit()
}

No installed package can transform to +proj=spilhaus on this machine. The table above is the whole result; the document stops here.

ENGINE <- capable[1]
cat("Coordinate transforms below use:", ENGINE, "\n")

Coordinate transforms below use: NA 

if (!RASTER_OK) cat("gdalraster is not loadable here, so the raster sections will be skipped.\n")

## forward/inverse constructors for any proj string, from the first capable package
make_forward <- function(proj_string) function(lonlat) engine_forward[[ENGINE]](lonlat, proj_string)
make_inverse <- function(proj_string) function(xy) engine_inverse[[ENGINE]](xy, proj_string)

## wrap: NaN where undefined, never error
safe <- function(f) function(m) {
  m <- matrix(as.numeric(m), ncol = 2)
  out <- tryCatch(f(m), error = function(e) matrix(NA_real_, nrow(m), 2))
  out <- matrix(as.numeric(out), ncol = 2)
  bad <- !is.finite(out[, 1]) | !is.finite(out[, 2]) | abs(out[, 1]) > 1e15 | abs(out[, 2]) > 1e15
  out[bad, ] <- NA_real_
  out
}
## inverse one point at a time so singular corners fail individually
pointwise <- function(f) function(m) {
  m <- matrix(as.numeric(m), ncol = 2)
  t(apply(m, 1, function(p) f(matrix(p, ncol = 2))))
}
forward <- safe(make_forward(SPILHAUS))
inverse <- pointwise(safe(make_inverse(SPILHAUS)))