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")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.
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))| package | version | gdal | proj | spilhaus |
|---|---|---|---|---|
| gdalraster | 2.6.1 | 3.13.3 | 9.9.0 | ok |
| PROJ | 0.7.0 | - | 9.9.0 | ok |
| sf | 1.1.3 | 3.13.3 | 9.9.0 | ok |
| terra | 1.9.50 | 3.13.3 | 9.9.0 | ok |
| 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")| tool | version |
|---|---|
| gdalinfo | GDAL 3.13.3 “Iowa City”, released 2026/08/13 |
| gdal | GDAL 3.13.3 “Iowa City”, released 2026/08/13 |
| proj | Rel. 9.9.0, September 15th, 2026 |
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()
}
ENGINE <- capable[1]
cat("Coordinate transforms below use:", ENGINE, "\n")Coordinate transforms below use: gdalraster 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)))2. The projection’s domain
PROJ does not publish the square’s half-width, so forward-project a dense lon/lat grid and take the extent; then inverse-project the square’s boundary to trace the cut on the globe.
square_half <- function(fwd, n = 721) {
g <- as.matrix(expand.grid(lon = seq(-180, 180, length.out = n), lat = seq(-90, 90, length.out = (n + 1) %/% 2)))
xy <- fwd(g)
ceiling(max(abs(xy), na.rm = TRUE) / 1000) * 1000
}
HALF <- square_half(forward)
n <- 300; t <- seq(-1, 1, length.out = n); edge <- 0.996 * HALF
boundary <- rbind(cbind(t * edge, -edge), cbind(edge, t * edge), cbind(rev(t) * edge, edge), cbind(-edge, rev(t) * edge))
cut <- inverse(boundary)
corners <- rbind(`corner A (Asia)` = c(115, 30), `corner B (South America)` = c(-65, -30))
cat(sprintf("Half-width %s m; %d of %d boundary samples inverted (misses cluster at the singular corners)\n",
format(HALF, big.mark = ","), sum(complete.cases(cut)), nrow(cut)))Half-width 11,803,000 m; 1172 of 1200 boundary samples inverted (misses cluster at the singular corners)op <- par(mar = c(3, 3, 2, 1))
plot(NA, xlim = c(-180, 180), ylim = c(-90, 90), asp = 1, xlab = "", ylab = "",
main = sprintf("Where the Spilhaus cut falls on the globe (half-width %.3f Mm)", HALF / 1e6))
grid(); points(cut, pch = 16, cex = 0.4, col = "#b8342c")
points(corners, cex = 1.6, lwd = 1.5); text(corners, rownames(corners), pos = 4, cex = 0.8)
par(op)3. Raster: bathymetry over /vsicurl
gdalraster::warp() is GDAL’s warper with the same option vocabulary as gdalwarp, so the target extent and size are stated explicitly rather than estimated. The result is written to /vsimem/ and read back as a matrix for ximage; nothing touches disk and only the tiles GDAL needs come over the wire.
warp_square <- function(proj_string, half, npix, url = RASTER_URL, out = "/vsimem/spilhaus.tif") {
t0 <- Sys.time()
gdalraster::warp(url, out, t_srs = proj_string,
cl_arg = c("-te", -half, -half, half, half, "-ts", npix, npix,
"-r", "average", "-dstalpha", "-overwrite", "-q"))
ds <- new(gdalraster::GDALRaster, out)
on.exit(ds$close())
v <- ds$read(band = 1, xoff = 0, yoff = 0, xsize = npix, ysize = npix, out_xsize = npix, out_ysize = npix)
a <- ds$read(band = ds$getRasterCount(), xoff = 0, yoff = 0, xsize = npix, ysize = npix, out_xsize = npix, out_ysize = npix)
m <- matrix(v, nrow = npix, ncol = npix, byrow = TRUE) # row 1 = top of the image
m[matrix(a, npix, npix, byrow = TRUE) == 0] <- NA
list(z = m, res_m = 2 * half / npix, seconds = round(as.numeric(difftime(Sys.time(), t0, units = "secs")), 1),
bands = ds$getRasterCount())
}
w <- warp_square(SPILHAUS, HALF, NPIX)
cat(sprintf("Warped %s to a %d x %d Spilhaus grid (%.1f km/px) in %s s\n",
basename(RASTER_URL), NPIX, NPIX, w$res_m / 1000, w$seconds))Warped GEBCO_2024.tif to a 1600 x 1600 Spilhaus grid (14.8 km/px) in 21.1 sbathy_cols <- grDevices::colorRampPalette(c("#08143f", "#123c7a", "#2a6fb5", "#6aa9d8", "#b9d9ea", "#e3f1f6"))(64)
render_raster <- function(z, half, add = FALSE) {
ocean <- z; ocean[!is.na(z) & z >= 0] <- NA
land <- z; land[is.na(z) | z < 0] <- NA
ximage::ximage(pmax(ocean, -7000), extent = c(-half, half, -half, half), asp = 1, col = bathy_cols,
zlim = c(-7000, 0), axes = FALSE, xlab = "", ylab = "", add = add)
ximage::ximage(matrix(as.numeric(!is.na(land)), nrow(land)), extent = c(-half, half, -half, half),
zlim = c(0, 1), col = c(NA, "#d9d2c5"), add = TRUE)
}
op <- par(mar = c(0, 0, 2, 0), bg = "#d9d2c5")
render_raster(w$z, HALF)Error in `loadNamespace()`:
! there is no package called 'ximage'title("Raster only (bathymetry), warped by GDAL via gdalraster", cex.main = 0.9)Error in `title()`:
! plot.new has not been called yetpar(op)4. Vector: Natural Earth ocean polygons, as vertices
Read remotely with gdalraster’s GDALVector, geometry returned as WKB, and handed to wk. From there the polygons are just a coordinate table with ring ids: transform the vertices, form segments between consecutive vertices of the same ring, and drop the segments that are absurdly long in projected space (those straddle the cut) or that are artefacts of the lon/lat representation (both ends on the antimeridian or at a pole). Coastlines are drawn as segments, not filled polygons.
lyr <- new(gdalraster::GDALVector, OCEAN_URL)
lyr$returnGeomAs <- "WKB"
feat <- lyr$fetch(-1)
lyr$close()
geom_col <- which(vapply(feat, is.list, logical(1)))[1] # the WKB column, whatever it is called
geom <- wk::wkb(feat[[geom_col]])
cat(sprintf("Ocean layer: %d feature(s), %s\n", length(geom), wk::wk_crs(geom)))
coast_segments <- function(geom, fwd, half, densify_deg = 0.5) {
## densify in lon/lat so the projected curves are smooth near the corners
## (wk has no densify; do it on the coordinate table)
co <- wk::wk_coords(geom)
key <- paste(co$feature_id, co$part_id, co$ring_id)
dens <- do.call(rbind, lapply(split(co, key), function(r) {
if (nrow(r) < 2) return(NULL)
steps <- pmax(1L, ceiling(sqrt(diff(r$x)^2 + diff(r$y)^2) / densify_deg))
x <- unlist(mapply(function(a, b, k) seq(a, b, length.out = k + 1)[-(k + 1)], r$x[-nrow(r)], r$x[-1], steps))
y <- unlist(mapply(function(a, b, k) seq(a, b, length.out = k + 1)[-(k + 1)], r$y[-nrow(r)], r$y[-1], steps))
data.frame(ring = r$ring_id[1] + 1e6 * r$feature_id[1], x = c(x, r$x[nrow(r)]), y = c(y, r$y[nrow(r)]))
}))
xy <- fwd(cbind(dens$x, dens$y))
same_ring <- dens$ring[-1] == dens$ring[-nrow(dens)]
step <- sqrt(diff(xy[, 1])^2 + diff(xy[, 2])^2)
synthetic <- (abs(dens$x[-1]) > 179.99 & abs(dens$x[-nrow(dens)]) > 179.99) |
(abs(dens$y[-1]) > 89.99 & abs(dens$y[-nrow(dens)]) > 89.99)
ok <- same_ring & is.finite(step) & step < half / 8 & !synthetic
i <- which(ok)
cbind(x0 = xy[i, 1], y0 = xy[i, 2], x1 = xy[i + 1, 1], y1 = xy[i + 1, 2])
}
segs <- coast_segments(geom, forward, HALF)
cat(sprintf("%d coastline segments after projecting with %s and dropping cut-crossing segments\n", nrow(segs), ENGINE))11992 coastline segments after projecting with gdalraster and dropping cut-crossing segments5. Points: landmarks, transformed directly
landmarks <- rbind(
Hobart = c(147.33, -42.88), Tokyo = c(139.69, 35.69), Reykjavik = c(-21.94, 64.15),
`Cape of Good Hope` = c(18.47, -34.36), `Bering Strait` = c(-169.0, 65.8),
`Drake Passage` = c(-63.0, -59.0), `Point Nemo` = c(-123.39, -48.88),
`Challenger Deep` = c(142.20, 11.37), `Puerto Rico Trench` = c(-66.5, 19.7),
`Mid-Atlantic Ridge (Azores)` = c(-27.0, 38.0), `Macquarie Island` = c(158.94, -54.62),
`Mawson Station` = c(62.87, -67.60), `Casey Station` = c(110.53, -66.28))
lm_xy <- forward(landmarks)
knitr::kable(data.frame(name = rownames(landmarks), lon = landmarks[, 1], lat = landmarks[, 2],
x = round(lm_xy[, 1]), y = round(lm_xy[, 2]), row.names = NULL),
caption = "Landmarks in Spilhaus metres")| name | lon | lat | x | y |
|---|---|---|---|---|
| Hobart | 147.33 | -42.88 | 2583897 | -1613781 |
| Tokyo | 139.69 | 35.69 | 9587907 | -52546 |
| Reykjavik | -21.94 | 64.15 | -5405098 | 6621233 |
| Cape of Good Hope | 18.47 | -34.36 | -2153244 | 421049 |
| Bering Strait | -169.00 | 65.80 | -4817477 | 11742840 |
| Drake Passage | -63.00 | -59.00 | -2017452 | -3916755 |
| Point Nemo | -123.39 | -48.88 | 1016488 | -5091561 |
| Challenger Deep | 142.20 | 11.37 | 6396508 | -83830 |
| Puerto Rico Trench | -66.50 | 19.70 | -9511346 | 3546491 |
| Mid-Atlantic Ridge (Azores) | -27.00 | 38.00 | -6037764 | 4511522 |
| Macquarie Island | 158.94 | -54.62 | 2012829 | -2209928 |
| Mawson Station | 62.87 | -67.60 | -168347 | -1008217 |
| Casey Station | 110.53 | -66.28 | 829890 | -1311409 |
6. Composite
composite <- function(z, half, segs, lm_xy, main) {
op <- par(mar = c(0, 0, 3, 0), bg = "white")
render_raster(z, half)
segments(segs[, 1], segs[, 2], segs[, 3], segs[, 4], col = "#3a3a3a", lwd = 0.6)
points(lm_xy, pch = 21, bg = "#f2b134", col = "black", cex = 1.1)
text(lm_xy, rownames(landmarks), pos = 4, cex = 0.75, offset = 0.4)
title(main, cex.main = 1)
par(op)
}
proj_used <- provenance$proj[provenance$package == ENGINE]
composite(w$z, HALF, segs, lm_xy,
sprintf("The world ocean in Spilhaus projection\nraster warped by GDAL (gdalraster); coastlines and points transformed by %s (PROJ %s)", ENGINE, proj_used))Error in `loadNamespace()`:
! there is no package called 'ximage'8. Another member of the family
(Skipped when the raster section was.)
Same method, different parameters: a naive recentring on the Atlantic. The square’s edges now slice the Pacific. Spilhaus’s contribution was the parameter choice, not the projection.
forward_alt <- safe(make_forward(SPILHAUS_ALT))
HALF_ALT <- square_half(forward_alt)
w_alt <- warp_square(SPILHAUS_ALT, HALF_ALT, max(400L, NPIX %/% 2L), out = "/vsimem/spilhaus_alt.tif")
segs_alt <- coast_segments(geom, forward_alt, HALF_ALT)
composite(w_alt$z, HALF_ALT, segs_alt, forward_alt(landmarks),
sprintf("%s\n(half-width %.3f Mm vs %.3f Mm for the default)", SPILHAUS_ALT, HALF_ALT / 1e6, HALF / 1e6))Error in `loadNamespace()`:
! there is no package called 'ximage'Appendix: the same thing in one line each
Convenience layers hide which PROJ did the work. Each of these either works or does not, purely according to the table in section 1.
if (requireNamespace("sf", quietly = TRUE)) {
x <- sf::st_transform(sf::st_sfc(sf::st_point(c(147.33, -42.88)), crs = "OGC:CRS84"), SPILHAUS)
print(sf::st_coordinates(x)); cat("sf used PROJ", sf::sf_extSoftVersion()[["PROJ"]], "\n")
} X Y
[1,] 2583897 -1613781
sf used PROJ 9.9.0 if (requireNamespace("terra", quietly = TRUE)) {
x <- terra::project(terra::vect(cbind(147.33, -42.88), crs = "OGC:CRS84"), SPILHAUS)
print(terra::crds(x)); cat("terra used PROJ", terra::gdal(lib = "proj"), "\n")
} x y
[1,] 2583897 -1613781
terra used PROJ 9.9.0 What to take from this
- On Linux the R packages usually agree with each other because they share one
libproj; the question is whether the system PROJ is 9.6 or later. A stock Ubuntu LTS says no. The matrix in this repo’s Actions shows what each runner image actually has. - On macOS and Windows the CRAN binaries carry their own PROJ, and those versions move on CRAN’s schedule.
- The primitives route (gdalraster for GDAL, wk for vertices, ximage for pixels) makes the dependency visible: every transform in this document is a call you can see, made by a package you can name.
- A proj-string is a recipe with hidden defaults. Cite the WKT with the parameters in it, and do not assume the one authority code you can find describes the map you meant.