replace proj_trans with wk::xy impl.
Closed this issue · 4 comments
@anthonynorth one thing I just don't understand atm is why the calls to P = proj_create_crs_to_crs(C, *crs_in, *crs_out, NULL);
accept a string like "+proj=laea" now in C_proj_trans.c, but the same calls in wk-trans.c and C_proj_crs_text.c do not.
i. e
proj_trans(0, 0, "+proj=laea +lon_0=147", source = "OGC:CRS84")
#$x_
#[1] -12230967
#$y_
# [1] 0
proj_crs_text("+proj=laea +lon_0=147")
##[1] "CONVERSION[\"PROJ-based coordinate operation\",\n METHOD[\"PROJ-based operation method: +proj=laea +lon_0=147\"]]"
proj_trans_create("OGC:CRS84", "+proj=laea +lon_0=147")
Error in proj_trans_create("OGC:CRS84", "+proj=laea +lon_0=147") :
Unknown error (code 4096)
proj_create_operations: target_crs is not a CRS or a CoordinateMetadata
Like, why does proj_trans work now? It works with the correct addition of +type=crs
and now is probably a good time to apply that as a requirement (we could help a user with warnings).
proj_crs_text("+proj=laea +lon_0=147 +type=crs")
proj_trans_create("OGC:CRS84", "+proj=laea +lon_0=147 +type=crs")
I'll proceed with that for now, but it leaves me wondering why it's working now in proj_trans ...
PROJ::proj_version()
[1] "9.3.1"
here's a rough benchmark for reference (is this the kind of difference you were seeing?)
library(PROJ)
## dummy calls to spin up the machinery, if that matters
dummy <- local({
proj_trans_create("OGC:CRS84", "+proj=laea +type=crs")
proj_trans(cbind(0, 0), "OGC:CRS84", source = "+proj=laea +type=crs +lon_0=147")
})
## -------------- pure xy matrix proj_trans
xyn <- geosphere::randomCoordinates(1e7)
system.time(proj_trans(xyn, "+proj=laea", source = "OGC:CRS84"))
#> user system elapsed
#> 3.424 0.303 3.727
## ---------------- transform pre-existing wk/xy object
xyn <- geosphere::randomCoordinates(1e7)
xy <- wk::xy(xyn[,1], xyn[,2], crs = "OGC:CRS84")
trans <- proj_trans_create("OGC:CRS84", "+proj=laea +type=crs")
system.time(pp <- wk::wk_transform(xy, trans))
#> user system elapsed
#> 4.183 0.099 4.283
## ---------------- transform wk/xy but include packing and unpacking
xyn <- geosphere::randomCoordinates(1e7)
system.time({
xy <- wk::xy(xyn[,1], xyn[,2], crs = "OGC:CRS84")
trans <- proj_trans_create("OGC:CRS84", "+proj=laea +type=crs")
pp <- wk::wk_coords(wk::wk_transform(xy, trans))[,c("x", "y")]
})
#> user system elapsed
#> 4.118 0.400 4.518
Created on 2024-02-07 with reprex v2.0.2
- replace proj_trans with wk version
- replace all PROJ.4 strings with "+type=crs" included
- delete C_proj_trans_list and C_proj_trans_xy
- test for bad input without "type=crs" in other functions than proj_trans
- test other bad inputs (this can be pretty liberal strings like "4326" and "NAD27" that currently crash PROJ (yikes))
- pass tests and read doc for sense ....
- NEWS: we now require "+type=crs" for old-style PROJ.4 strings
all the test message changes need a bit of thought, so parking this for now: https://github.com/hypertidy/PROJ/tree/proj_trans-via-wk
It could be a bug. I'll take a look.
here's a rough benchmark for reference (is this the kind of difference you were seeing?)
Yes. I was seeing an wk_transform()
execution time overhead somewhere in the 15% ballpark, but it depends on the transform being applied. wk_transform()
allocates less than half the R memory also.
That overhead amounts to ~ 0.3s for 10 million coordinates in my environment.