hypertidy/PROJ

replace proj_trans with wk::xy impl.

Closed this issue · 4 comments

replace proj_trans with wk::xy impl.

@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.