Skip to content

Instantly share code, notes, and snippets.

@asinghvi17
Created April 14, 2024 14:18
Show Gist options
  • Save asinghvi17/6553bec4cb18acf254d3a714142899bf to your computer and use it in GitHub Desktop.
Save asinghvi17/6553bec4cb18acf254d3a714142899bf to your computer and use it in GitHub Desktop.
How to use `proj_trans_generic`
using Proj, BenchmarkTools
import GeometryOps as GO
# using proj_trans_generic on a vector of tuples - or could be GB.Point or SVector or GI.Point with arbitrary metadata, doesn't matter!
points = tuple.(rand(10000), rand(10000))
t = Proj.Transformation("+proj=longlat", "+proj=natearth"; always_xy = true)
function batch_project_numsfirst!(points, projection)
# GC.@preserve points begin # not technically necessary
ptr = pointer(points) # we use `pointer` and not `Ref` here because we need to do pointer arithmetic.
x_ptr = ptr
y_ptr = ptr + sizeof(Float64) # see above, we do pointer arithmetic here to offset y by 1 Float64
stride = sizeof(first(points)) # it's assumed this is equal for all points
n_points = length(points) # regular C pointer-as-array semantics
z_ptr, z_stride, z_num = if GI.is3d(first(points)) # only populate Z if 3D.
ptr + 2*sizeof(Float64), stride, length(points)
else
typeof(ptr)(C_NULL), 0, 0
end
@ccall Proj.libproj.proj_trans_generic(
projection.pj::Ptr{Cvoid}, # P
projection.direction::Proj.PJ_DIRECTION,# direction
x_ptr::Ptr{Cdouble}, # x
stride::Csize_t, # sx
n_points::Csize_t, # nx
y_ptr::Ptr{Cdouble}, # y
stride::Csize_t, # sy
n_points::Csize_t, # ny
z_ptr::Ptr{Cdouble}, # z
z_stride::Csize_t, # sz
z_num::Csize_t, # nz
C_NULL::Ptr{Cdouble}, # t - we ignore it here, but could potentially include
0::Csize_t, # st
0::Csize_t, # nt
)::Csize_t
# end
return points
end
batch_project_numsfirst!(tuple.(rand(10000), rand(10000)), t)
@benchmark batch_project_numsfirst!(dt, $t) setup=(dt=tuple.(rand(10000), rand(10000))) evals=1
@benchmark GO.reproject(points, t)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment