projr

Prototype for exploring what a full PROJ R package feels like. Three transformation interfaces are provided to compare — use whichever feels right, or all of them.

Setup

# After cloning, from inside the package directory:
cpp11::cpp_register()   # generates src/cpp11.cpp and R/cpp11.R
devtools::load_all()

Interface comparison

Option A: typed variants (proj_trans_xy, _xyz, _xyzt)

The caller chooses the function, the output is always exactly that shape.

# 2D — the common case
proj_trans_xy(147:148, c(-42, -45), "OGC:CRS84", "EPSG:32755")
#> $x
#> [1] 321932.5 ...
#> $y
#> [1] 5345455 ...

# 3D — geocentric
proj_trans_xyz(147, -42, 100, "OGC:CRS84", "EPSG:4978")

# 4D — time-dependent
proj_trans_xyzt(147, -42, 100, 2025.0,
                "EPSG:7844", "EPSG:9309")  # GDA2020 -> ITRF2020

Pros: No ambiguity about what z/t mean. Fast — no optional-arg dispatch. No surprise columns in output. Easy to grep/autocomplete.

Cons: 3 functions for one concept. Users who just want "transform" have to pick. The _xyzt form is rare enough to feel odd.

Option B: single function with optional z/t (proj_trans)

Output shape mirrors input. Give it xy, get xy back. Give it xyz, get xyz.

# 2D
proj_trans(147:148, c(-42, -45), source = "OGC:CRS84", target = "EPSG:32755")
#> $x  $y    (2 elements)

# 3D — just add z
proj_trans(147, -42, z = 100, source = "OGC:CRS84", target = "EPSG:4978")
#> $x  $y  $z  (3 elements)

# 4D
proj_trans(147, -42, z = 100, t = 2025.0,
           source = "EPSG:7844", target = "EPSG:9309")
#> $x  $y  $z  $t  (4 elements)

Pros: One function to learn. Output is minimal — no z column when you didn't ask for z. Natural for 2D (the 90% case).

Cons: z and t are keyword-only which is correct but slightly more typing. The variable-shape output might surprise programmatic users who expect a fixed schema. Slight overhead from NULL checks (negligible in practice).

Option C: with direction (proj_transform)

Like Option B but adds a direction argument for inverse transforms.

# Forward
pt <- proj_transform(147, -42, source = "OGC:CRS84", target = "EPSG:32755")

# Inverse — same source/target, flip direction
proj_transform(pt$x, pt$y, source = "OGC:CRS84", target = "EPSG:32755",
               direction = "inverse")

Pros: Clean inverse without swapping source/target strings. Matches the PROJ C API's PJ_FWD/PJ_INV model. Useful for pipelines.

Cons: Another argument to learn. For CRS-to-CRS transforms, just swapping source/target is arguably clearer. Direction is more meaningful for pipeline strings like +proj=pipeline +step ....

What about matrix input?

All three options work with vectors. For matrix/data.frame input, a thin wrapper in R is cleaner than putting it in C++:

# Not exported yet — just showing the idea
proj_trans_matrix <- function(m, source, target) {
  stopifnot(ncol(m) >= 2)
  result <- if (ncol(m) == 2) {
    proj_trans_xy(m[,1], m[,2], source, target)
  } else if (ncol(m) == 3) {
    proj_trans_xyz(m[,1], m[,2], m[,3], source, target)
  } else {
    proj_trans_xyzt(m[,1], m[,2], m[,3], m[,4], source, target)
  }
  do.call(cbind, result)
}

But this is a UX question for later — the C layer just works with vectors.

CRS objects

crs <- proj_crs("EPSG:32755")
crs
#> <projr_crs> WGS 84 / UTM zone 55S
#>   Type: ProjectedCRS
#>   ID:   EPSG:32755

proj_get_type(crs)
#> [1] "ProjectedCRS"

proj_as_wkt(crs, "WKT1_GDAL")
proj_as_projjson(crs)

# Decompose
gcrs <- proj_crs_get_geodetic_crs(crs)
gcrs
#> <projr_crs> WGS 84
#>   Type: GeographicCRS
#>   ID:   EPSG:4326

proj_get_ellipsoid(crs)
#> $name [1] "WGS 84"
#> $semi_major [1] 6378137
#> $semi_minor [1] 6356752.3
#> $inv_flattening [1] 298.2572

# Identify an unknown CRS
mystery <- proj_crs("+proj=longlat +datum=WGS84 +type=crs")
proj_identify(mystery)
#>   authority code confidence
#> 1      EPSG 4326        70

Projection factors (for tissot)

proj_factors(seq(140, 155, by = 5), rep(-42, 4),
             "+proj=utm +zone=55 +datum=WGS84")
#>   meridional_scale parallel_scale areal_scale angular_distortion ...

Database queries

proj_get_authorities()
#> [1] "EPSG" "ESRI" "IAU_2015" "IGNF" "NKG" "OGC"

head(proj_get_codes("EPSG", "Ellipsoid"))
#> [1] "7001" "7002" "7003" ...

proj_db_path()
#> [1] "/usr/share/proj/proj.db"

Bounds

proj_trans_bounds(140, -50, 155, -35, "OGC:CRS84", "EPSG:32755")
#>       xmin       ymin       xmax       ymax
#> -611021.4 -5559753.0  945396.1 -3804479.3