From abee4e39d138d91cfa42234d2ec394f1d4d5a9a8 Mon Sep 17 00:00:00 2001
From: rafaqz <rafaelschouten@gmail.com>
Date: Wed, 29 Mar 2023 00:08:25 +0100
Subject: [PATCH 01/12] add geointerface reprojection

---
 Project.toml         |  2 +
 src/Proj.jl          |  2 +
 src/geointerface.jl  | 98 ++++++++++++++++++++++++++++++++++++++++++++
 test/geointerface.jl | 14 +++++++
 4 files changed, 116 insertions(+)
 create mode 100644 src/geointerface.jl
 create mode 100644 test/geointerface.jl

diff --git a/Project.toml b/Project.toml
index 923f404..d4125e4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,6 +6,7 @@ version = "1.4.0"
 CEnum = "fa961155-64e5-5f13-b03f-caf6b980ea82"
 CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
 GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
+GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 NetworkOptions = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
 PROJ_jll = "58948b4f-47e0-5654-a9ad-f609743f8632"
 
@@ -13,6 +14,7 @@ PROJ_jll = "58948b4f-47e0-5654-a9ad-f609743f8632"
 CEnum = "0.2, 0.3, 0.4"
 CoordinateTransformations = "0.6"
 GeoFormatTypes = "0.4"
+GeoInterface = "1"
 PROJ_jll = "900.100"
 julia = "1.6"
 
diff --git a/src/Proj.jl b/src/Proj.jl
index 7a8fedf..7dfbfc6 100644
--- a/src/Proj.jl
+++ b/src/Proj.jl
@@ -3,6 +3,7 @@ module Proj
 using PROJ_jll
 using CEnum
 using CoordinateTransformations
+using GeoInterface
 using NetworkOptions: ca_roots
 import GeoFormatTypes as GFT
 
@@ -92,6 +93,7 @@ include("libproj.jl")
 include("crs.jl")
 include("coord.jl")
 include("error.jl")
+include("geointerface.jl")
 
 """
     unsafe_loadstringlist(ptr::Ptr{Cstring})
diff --git a/src/geointerface.jl b/src/geointerface.jl
new file mode 100644
index 0000000..e7d013d
--- /dev/null
+++ b/src/geointerface.jl
@@ -0,0 +1,98 @@
+using GeoInterface
+
+const GI = GeoInterface
+
+# Coord is a PointTrait
+GI.isgeometry(::Type{Coord}) = true
+GI.geomtrait(::Coord) = GI.PointTrait()
+GI.x(::PointTrait, c::Coord) = c.x
+GI.y(::PointTrait, c::Coord) = c.y
+GI.z(::PointTrait, c::Coord) = c.z
+GI.getcoord(::PointTrait, c::Coord) = (c.x, c.y, c.z)
+GI.getcoord(::PointTrait, c::Coord, i::Int) = c[i]
+coordnames(::PointTrait, ::Coord) = (:X, :Y, :Z)
+
+"""
+    reproject(geometry, source_crs, target_crs)
+
+Reproject any GeoInterface.jl compatible `geometry` from `source_crs` to `target_crs`.
+
+The returned object will be constructed from `GeoInterface.WrapperGeometry`
+geometries, wrapping views of `Proj.Coord`.
+"""
+function reproject(geom, source_crs, target_crs; time=Inf)
+    source_crs1 = convert(Proj.CRS, source_crs)
+    target_crs1 = convert(Proj.CRS, target_crs)
+    trans = Proj.Transformation(source_crs1, target_crs1, always_xy = true)
+    coords = if GI.is3d(geom)
+        [Proj.Coord(GI.x(p), GI.y(p), GI.z(p), time) for p in GI.getpoint(geom)]
+    else
+        [Proj.Coord(GI.x(p), GI.y(p), 0.0, time) for p in GI.getpoint(geom)]
+    end
+    err = Proj.proj_trans_array(trans.pj, Proj.PJ_FWD, length(coords), coords)
+    err == 0 || error("Proj error $err")
+
+    crs = target_crs isa GFT.GeoFormat ? target_crs : convert(WellKnownText, crs)
+    return reconstruct(geom, coords; crs)
+end
+
+
+"""
+    reconstruct(geometry, points::AbstractVector)
+
+Reconstruct GeoInterface compatible `geometry` from the `PointTrait` geometries in `points`.
+
+The returned object will be constructed from `GeoInterface.WrapperGeometry`
+geometries, wrapping views into the `points` vector. Otherwise, the structure and
+GeoInterface traits will be the same as for `geometry`.
+
+This is a lazy operation, only allocating for outer wrappers of nested geometries.
+Later changes to the `points` vector *will* affect the returned geometry.
+
+`GeoInterface.npoint(geometry) == length(points)` must be `true`.
+"""
+function reconstruct(geom, points::AbstractVector; crs=nothing)
+    T = nonmissingtype(eltype(points))
+    isgeometry(T) || throw(ArgumentError("points of type $(T) are not GeoInterface.jl compatible objects"))
+    trait = GI.trait(first(skipmissing(points))) 
+    trait isa PointTrait || throw(ArgumentError("can only reconstruct from an array of points, got $trait"))
+    n, geom = _reconstruct(geom, points, 0, crs)
+    @assert n == GI.npoint(geom)
+    return geom
+end
+
+_reconstruct(geom, points, n::Int, crs) = _reconstruct(GI.trait(geom), geom, points, n::Int, crs)
+# Nested geometries. We need to reconstruct their child geoms.
+function _reconstruct(trait::AbstractGeometryTrait, geom, points, n, crs)
+    T = GI.geointerface_geomtype(trait)
+    geoms = map(GI.getgeom(geom)) do childgeom
+        n, reconstructed_childgeom = _reconstruct(GI.trait(childgeom), childgeom, points, n, crs)
+        reconstructed_childgeom
+    end
+    return n, T(geoms; crs)
+end
+# Bottom-level geometries that can wrap a vector of points.
+function _reconstruct(
+    trait::Union{GI.LineTrait,GI.LineStringTrait,GI.LinearRingTrait,GI.MultiPointTrait}, 
+    geom, points, n, crs
+)
+    T = GI.geointerface_geomtype(trait)
+    npoints = GI.npoint(geom)
+    v = view(points, n + 1:n + npoints)
+    geom = if GI.is3d(geom)
+        GI.ismeasured(geom) ? T{true,true}(v; crs) : T{true,false}(v; crs)
+    else
+        GI.ismeasured(geom) ? T{false,true}(v; crs) : T{false,false}(v; crs)
+    end
+    return n + npoints, geom
+end
+# Points
+function _reconstruct(::GI.PointTrait, geom, points, n, crs)
+    n += 1
+    point = if GI.is3d(geom)
+        GI.ismeasured(geom) ? GI.Point{true,true}(points[n]) : GI.Point{true,false}(points[n])
+    else
+        GI.ismeasured(geom) ? GI.Point{false,true}(points[n]) : GI.Point{false,false}(points[n])
+    end
+    return n, point
+end
diff --git a/test/geointerface.jl b/test/geointerface.jl
new file mode 100644
index 0000000..c77816a
--- /dev/null
+++ b/test/geointerface.jl
@@ -0,0 +1,14 @@
+using Proj
+using GeoInterface
+using GeoFormatTypes
+const GI = GeoInterface
+
+linearring1 = GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)])
+linearring2 = GI.LinearRing([(11, 2), (13, 4), (15, 6), (11, 2)])
+
+polygon1 = GI.Polygon([linearring1])
+polygon2 = GI.Polygon([linearring2])
+multipolygon = GI.MultiPolygon([polygon1, polygon2])
+
+Proj.reproject(multipolygon, EPSG(4326), EPSG(3857))
+

From 6966a3f46e9954fa54d6099065a71f1fe314fa56 Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Fri, 31 Mar 2023 19:39:21 +0200
Subject: [PATCH 02/12] add kwargs reproject

---
 src/geointerface.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/geointerface.jl b/src/geointerface.jl
index e7d013d..4eb10cd 100644
--- a/src/geointerface.jl
+++ b/src/geointerface.jl
@@ -20,6 +20,8 @@ Reproject any GeoInterface.jl compatible `geometry` from `source_crs` to `target
 The returned object will be constructed from `GeoInterface.WrapperGeometry`
 geometries, wrapping views of `Proj.Coord`.
 """
+reproject(geom; source_crs, target_crs, time=Inf) =
+    reproject(geom, source_crs, target_crs; time)
 function reproject(geom, source_crs, target_crs; time=Inf)
     source_crs1 = convert(Proj.CRS, source_crs)
     target_crs1 = convert(Proj.CRS, target_crs)

From cf5d1642846268f276cf6c275952e9d91af13298 Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Fri, 31 Mar 2023 22:55:10 +0200
Subject: [PATCH 03/12] bugfix

---
 src/geointerface.jl | 57 ++++++++++++++++++++++++---------------------
 1 file changed, 30 insertions(+), 27 deletions(-)

diff --git a/src/geointerface.jl b/src/geointerface.jl
index 4eb10cd..6eb8542 100644
--- a/src/geointerface.jl
+++ b/src/geointerface.jl
@@ -7,10 +7,10 @@ GI.isgeometry(::Type{Coord}) = true
 GI.geomtrait(::Coord) = GI.PointTrait()
 GI.x(::PointTrait, c::Coord) = c.x
 GI.y(::PointTrait, c::Coord) = c.y
-GI.z(::PointTrait, c::Coord) = c.z
-GI.getcoord(::PointTrait, c::Coord) = (c.x, c.y, c.z)
+# GI.z(::PointTrait, c::Coord) = c.z
+GI.getcoord(::PointTrait, c::Coord) = (c.x, c.y)
 GI.getcoord(::PointTrait, c::Coord, i::Int) = c[i]
-coordnames(::PointTrait, ::Coord) = (:X, :Y, :Z)
+GI.coordnames(::PointTrait, ::Coord) = (:X, :Y)
 
 """
     reproject(geometry, source_crs, target_crs)
@@ -23,26 +23,25 @@ geometries, wrapping views of `Proj.Coord`.
 reproject(geom; source_crs, target_crs, time=Inf) =
     reproject(geom, source_crs, target_crs; time)
 function reproject(geom, source_crs, target_crs; time=Inf)
+    wrapped_geom, coords = reconstruct(geom) do p
+        Proj.Coord(GI.x(p), GI.y(p), (GI.is3d(geom) ? GI.z(p) : 0.0), time)
+    end
     source_crs1 = convert(Proj.CRS, source_crs)
     target_crs1 = convert(Proj.CRS, target_crs)
     trans = Proj.Transformation(source_crs1, target_crs1, always_xy = true)
-    coords = if GI.is3d(geom)
-        [Proj.Coord(GI.x(p), GI.y(p), GI.z(p), time) for p in GI.getpoint(geom)]
-    else
-        [Proj.Coord(GI.x(p), GI.y(p), 0.0, time) for p in GI.getpoint(geom)]
-    end
     err = Proj.proj_trans_array(trans.pj, Proj.PJ_FWD, length(coords), coords)
     err == 0 || error("Proj error $err")
-
-    crs = target_crs isa GFT.GeoFormat ? target_crs : convert(WellKnownText, crs)
-    return reconstruct(geom, coords; crs)
+    # crs = target_crs isa GFT.GeoFormat ? target_crs : convert(WellKnownText, crs)
+    return wrapped_geom
 end
 
 
 """
+    reconstruct(f::Function, geom)
     reconstruct(geometry, points::AbstractVector)
 
-Reconstruct GeoInterface compatible `geometry` from the `PointTrait` geometries in `points`.
+Reconstruct GeoInterface compatible `geometry` from the `PointTrait` geometries
+in `points`, or from the result of function `f` over the points in `geom`.
 
 The returned object will be constructed from `GeoInterface.WrapperGeometry`
 geometries, wrapping views into the `points` vector. Otherwise, the structure and
@@ -53,48 +52,52 @@ Later changes to the `points` vector *will* affect the returned geometry.
 
 `GeoInterface.npoint(geometry) == length(points)` must be `true`.
 """
-function reconstruct(geom, points::AbstractVector; crs=nothing)
-    T = nonmissingtype(eltype(points))
+function reconstruct(f, geom; crs=nothing)
+    p1 = f(first(GI.getpoint(geom)))
+    T = typeof(p1)
     isgeometry(T) || throw(ArgumentError("points of type $(T) are not GeoInterface.jl compatible objects"))
-    trait = GI.trait(first(skipmissing(points))) 
+    trait = GI.trait(p1)
     trait isa PointTrait || throw(ArgumentError("can only reconstruct from an array of points, got $trait"))
-    n, geom = _reconstruct(geom, points, 0, crs)
+    points = Vector{T}(undef, GI.npoint(geom))
+    n, geom = _reconstruct!(f, points, geom, 0, crs)
     @assert n == GI.npoint(geom)
-    return geom
+    return geom, points
 end
 
-_reconstruct(geom, points, n::Int, crs) = _reconstruct(GI.trait(geom), geom, points, n::Int, crs)
+_reconstruct!(f, points, geom, n::Int, crs) = _reconstruct!(f, points, GI.trait(geom), geom, n::Int, crs)
 # Nested geometries. We need to reconstruct their child geoms.
-function _reconstruct(trait::AbstractGeometryTrait, geom, points, n, crs)
+function _reconstruct!(f, points, trait::AbstractGeometryTrait, geom, n, crs)
     T = GI.geointerface_geomtype(trait)
     geoms = map(GI.getgeom(geom)) do childgeom
-        n, reconstructed_childgeom = _reconstruct(GI.trait(childgeom), childgeom, points, n, crs)
+        n, reconstructed_childgeom = _reconstruct!(f, points, GI.trait(childgeom), childgeom, n, crs)
         reconstructed_childgeom
     end
     return n, T(geoms; crs)
 end
 # Bottom-level geometries that can wrap a vector of points.
-function _reconstruct(
+function _reconstruct!(f, points, 
     trait::Union{GI.LineTrait,GI.LineStringTrait,GI.LinearRingTrait,GI.MultiPointTrait}, 
-    geom, points, n, crs
+    geom, n, crs
 )
     T = GI.geointerface_geomtype(trait)
     npoints = GI.npoint(geom)
     v = view(points, n + 1:n + npoints)
-    geom = if GI.is3d(geom)
+    v .= f.(GI.getpoint(geom))
+    wrapped_geom = if GI.is3d(geom)
         GI.ismeasured(geom) ? T{true,true}(v; crs) : T{true,false}(v; crs)
     else
         GI.ismeasured(geom) ? T{false,true}(v; crs) : T{false,false}(v; crs)
     end
-    return n + npoints, geom
+    return n + npoints, wrapped_geom
 end
 # Points
-function _reconstruct(::GI.PointTrait, geom, points, n, crs)
+function _reconstruct!(f, points, ::GI.PointTrait, geom, n, crs)
     n += 1
-    point = if GI.is3d(geom)
+    points[n] = f(point)
+    wrapped_point = if GI.is3d(geom)
         GI.ismeasured(geom) ? GI.Point{true,true}(points[n]) : GI.Point{true,false}(points[n])
     else
         GI.ismeasured(geom) ? GI.Point{false,true}(points[n]) : GI.Point{false,false}(points[n])
     end
-    return n, point
+    return n, wrapped_point
 end

From 2aea0b29d4231799a2043d1c294ec9614acd5c20 Mon Sep 17 00:00:00 2001
From: rafaqz <rafaelschouten@gmail.com>
Date: Mon, 3 Apr 2023 00:19:23 +0100
Subject: [PATCH 04/12] better tests

---
 test/geointerface.jl | 23 ++++++++++++++++++-----
 1 file changed, 18 insertions(+), 5 deletions(-)

diff --git a/test/geointerface.jl b/test/geointerface.jl
index c77816a..762e5ee 100644
--- a/test/geointerface.jl
+++ b/test/geointerface.jl
@@ -1,14 +1,27 @@
 using Proj
+using Test
 using GeoInterface
 using GeoFormatTypes
 const GI = GeoInterface
+using Plots
 
-linearring1 = GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)])
-linearring2 = GI.LinearRing([(11, 2), (13, 4), (15, 6), (11, 2)])
+using GeoInterfaceRecipes
+GeoInterfaceRecipes.@enable_geo_plots GeoInterface.Wrappers.WrapperGeometry
 
-polygon1 = GI.Polygon([linearring1])
-polygon2 = GI.Polygon([linearring2])
+ring1 = GI.LinearRing([(1, 2), (7, 4), (5, 6), (1, 2)])
+ring2 = GI.LinearRing([(11, 2), (20, 4), (15, 6), (11, 2)])
+hole2 = GI.LinearRing([(14, 4), (16, 4), (17, 5), (14, 4)])
+
+polygon1 = GI.Polygon([ring1])
+polygon2 = GI.Polygon([ring2, hole2])
 multipolygon = GI.MultiPolygon([polygon1, polygon2])
 
-Proj.reproject(multipolygon, EPSG(4326), EPSG(3857))
+multipolygon3857 = Proj.reproject(multipolygon, EPSG(4326), EPSG(3857))
+multipolygon4326 = Proj.reproject(multipolygon3857; target_crs=EPSG(4326))
+points1 = collect(GI.getpoint(multipolygon))
+points2 = GI.convert.(Tuple, GI.getpoint(multipolygon4326))
+# Round trip comparison
+@test all(map((p1, p2) -> all(map(isapprox, p1, p2)), points1, points2))
+@test GI.crs(multipolygon3857) == EPSG(3857)
+@test GI.crs(multipolygon4326) == EPSG(4326)
 

From 0e5f2673faa21e277a2992d77354314d23bee59b Mon Sep 17 00:00:00 2001
From: rafaqz <rafaelschouten@gmail.com>
Date: Mon, 3 Apr 2023 00:19:56 +0100
Subject: [PATCH 05/12] Coord 2 and 3

---
 src/Proj.jl         |  8 +++++--
 src/geointerface.jl | 43 +++++++++++++++++++++++++++-----------
 src/libproj.jl      | 51 ++++++++++++++++++++++++++-------------------
 3 files changed, 66 insertions(+), 36 deletions(-)

diff --git a/src/Proj.jl b/src/Proj.jl
index 7dfbfc6..91e4a7a 100644
--- a/src/Proj.jl
+++ b/src/Proj.jl
@@ -50,13 +50,17 @@ julia> Proj.Coord((1.0, 2.0, 3.0))
  Inf
 ```
 """
-struct Coord <: AbstractVector{Float64}
+struct Coord{N} <: AbstractVector{Float64}
     x::Float64
     y::Float64
     z::Float64
     t::Float64
-    Coord(x, y, z = 0.0, t = Inf) = new(x, y, z, t)
 end
+Coord{3}(x, y, z=0.0; time=Inf) = Coord{3}(x, y, z, time)
+Coord{2}(x, y; time=Inf) = Coord{2}(x, y, 0.0, time)
+Coord(x, y; time=Inf) = Coord{2}(x, y, 0.0, time)
+Coord(x, y, z; time=Inf) = Coord{3}(x, y, z, time)
+Coord(x, y, z, t) = Coord{3}(x, y, z, t)
 
 # this shields a StackOverflow from the splatting constructor
 Coord(::Real) = error("Proj.Coord takes 2 to 4 numbers, one given")
diff --git a/src/geointerface.jl b/src/geointerface.jl
index 6eb8542..3ce2877 100644
--- a/src/geointerface.jl
+++ b/src/geointerface.jl
@@ -5,30 +5,48 @@ const GI = GeoInterface
 # Coord is a PointTrait
 GI.isgeometry(::Type{Coord}) = true
 GI.geomtrait(::Coord) = GI.PointTrait()
-GI.x(::PointTrait, c::Coord) = c.x
-GI.y(::PointTrait, c::Coord) = c.y
-# GI.z(::PointTrait, c::Coord) = c.z
-GI.getcoord(::PointTrait, c::Coord) = (c.x, c.y)
-GI.getcoord(::PointTrait, c::Coord, i::Int) = c[i]
-GI.coordnames(::PointTrait, ::Coord) = (:X, :Y)
+GI.x(::PointTrait, c::Coord) = c.x # TODO: these may have the wrong order for some CRS. 
+GI.y(::PointTrait, c::Coord) = c.y 
+GI.z(::PointTrait, c::Coord{3}) = c.z
+GI.is3d(::PointTrait, c::Coord{2}) = false
+GI.is3d(::PointTrait, c::Coord{3}) = true
+GI.ismeasured(::PointTrait, c::Coord) = false
+function GI.getcoord(::PointTrait, c::Coord{N}, i::Int) where N
+    checkbounds(1:N, i)
+    @inbounds c[i]
+end
+GI.getcoord(::PointTrait, c::Coord{2}) = (c.x, c.y)
+GI.getcoord(::PointTrait, c::Coord{3}) = (c.x, c.y, c.z)
+coordnames(::PointTrait, ::Coord{2}) = (:X, :Y)
+coordnames(::PointTrait, ::Coord{3}) = (:X, :Y, :Z)
 
 """
     reproject(geometry, source_crs, target_crs)
+    reproject(geometry; [source_crs,] target_crs)
 
 Reproject any GeoInterface.jl compatible `geometry` from `source_crs` to `target_crs`.
 
 The returned object will be constructed from `GeoInterface.WrapperGeometry`
 geometries, wrapping views of `Proj.Coord`.
 """
-reproject(geom; source_crs, target_crs, time=Inf) =
-    reproject(geom, source_crs, target_crs; time)
-function reproject(geom, source_crs, target_crs; time=Inf)
-    wrapped_geom, coords = reconstruct(geom) do p
-        Proj.Coord(GI.x(p), GI.y(p), (GI.is3d(geom) ? GI.z(p) : 0.0), time)
+function reproject(geom; source_crs=nothing, target_crs, kw...)
+    source_crs = isnothing(source_crs) ? GeoInterface.crs(geom) : source_crs
+    isnothing(source_crs) && throw(ArgumentError("geom has no crs attatched. Pass a `source_crs` keyword"))
+    reproject(geom, source_crs, target_crs; kw...)
+end
+function reproject(geom, source_crs, target_crs; time=Inf, always_xy=true)
+    wrapped_geom, coords = if GI.is3d(geom)
+        reconstruct(geom; crs=target_crs) do p
+            Proj.Coord{3}(GI.x(p), GI.y(p), GI.z(p); time)
+        end
+    else
+        reconstruct(geom; crs=target_crs) do p
+            Proj.Coord{2}(GI.x(p), GI.y(p); time)
+        end
     end
     source_crs1 = convert(Proj.CRS, source_crs)
     target_crs1 = convert(Proj.CRS, target_crs)
-    trans = Proj.Transformation(source_crs1, target_crs1, always_xy = true)
+    trans = Proj.Transformation(source_crs1, target_crs1; always_xy)
     err = Proj.proj_trans_array(trans.pj, Proj.PJ_FWD, length(coords), coords)
     err == 0 || error("Proj error $err")
     # crs = target_crs isa GFT.GeoFormat ? target_crs : convert(WellKnownText, crs)
@@ -69,6 +87,7 @@ _reconstruct!(f, points, geom, n::Int, crs) = _reconstruct!(f, points, GI.trait(
 function _reconstruct!(f, points, trait::AbstractGeometryTrait, geom, n, crs)
     T = GI.geointerface_geomtype(trait)
     geoms = map(GI.getgeom(geom)) do childgeom
+        childcrs = nothing # We dont pass the CRS down, the type gets too complicated
         n, reconstructed_childgeom = _reconstruct!(f, points, GI.trait(childgeom), childgeom, n, crs)
         reconstructed_childgeom
     end
diff --git a/src/libproj.jl b/src/libproj.jl
index 11caf55..89bb96a 100644
--- a/src/libproj.jl
+++ b/src/libproj.jl
@@ -490,20 +490,20 @@ function proj_degree_output(P, dir)
     @ccall libproj.proj_degree_output(P::Ptr{PJ}, dir::PJ_DIRECTION)::Cint
 end
 
-function proj_trans(P, direction, coord)
-    @ccall libproj.proj_trans(P::Ptr{PJ}, direction::PJ_DIRECTION, coord::Coord)::Coord
+function proj_trans(P, direction, coord::Coord{N}) where N
+    @ccall libproj.proj_trans(P::Ptr{PJ}, direction::PJ_DIRECTION, coord::Coord{N})::Coord{N}
 end
 
 function proj_trans_get_last_used_operation(P)
     @ccall libproj.proj_trans_get_last_used_operation(P::Ptr{PJ})::Ptr{PJ}
 end
 
-function proj_trans_array(P, direction, n, coord)
+function proj_trans_array(P, direction, n, coord::Vector{Coord{N}}) where N
     @ccall libproj.proj_trans_array(
         P::Ptr{PJ},
         direction::PJ_DIRECTION,
         n::Csize_t,
-        coord::Ptr{Coord},
+        coord::Ptr{Coord{N}},
     )::Cint
 end
 
@@ -566,37 +566,44 @@ end
 
 ` Doxygen_Suppress `
 """
-function proj_coord(x, y, z = 0.0, t = Inf)
-    @ccall libproj.proj_coord(x::Cdouble, y::Cdouble, z::Cdouble, t::Cdouble)::Coord
+function proj_coord(x, y; time=Inf)
+    z = 0.0 
+    @ccall libproj.proj_coord(x::Cdouble, y::Cdouble, z::Cdouble, t::Cdouble)::Coord{2}
+end
+function proj_coord(x, y, z; time=Inf)
+    @ccall libproj.proj_coord(x::Cdouble, y::Cdouble, z::Cdouble, t::Cdouble)::Coord{3}
+end
+function proj_coord(x, y, z, t)
+    @ccall libproj.proj_coord(x::Cdouble, y::Cdouble, z::Cdouble, t::Cdouble)::Coord{3}
 end
 
-function proj_roundtrip(P, direction, n, coord)
+function proj_roundtrip(P, direction, n, coord::Coord{N}) where N
     @ccall libproj.proj_roundtrip(
         P::Ptr{PJ},
         direction::PJ_DIRECTION,
         n::Cint,
-        coord::Ptr{Coord},
+        coord::Ptr{Coord{N}},
     )::Cdouble
 end
 
-function proj_lp_dist(P, a, b)
-    @ccall libproj.proj_lp_dist(P::Ptr{PJ}, a::Coord, b::Coord)::Cdouble
+function proj_lp_dist(P, a::Coord{Na}, b::Coord{Nb}) where {Na,Nb}
+    @ccall libproj.proj_lp_dist(P::Ptr{PJ}, a::Coord{Na}, b::Coord{Nb})::Cdouble
 end
 
-function proj_lpz_dist(P, a, b)
-    @ccall libproj.proj_lpz_dist(P::Ptr{PJ}, a::Coord, b::Coord)::Cdouble
+function proj_lpz_dist(P, a::Coord{Na}, b::Coord{Nb}) where {Na,Nb}
+    @ccall libproj.proj_lpz_dist(P::Ptr{PJ}, a::Coord{Na}, b::Coord{Nb})::Cdouble
 end
 
-function proj_xy_dist(a, b)
-    @ccall libproj.proj_xy_dist(a::Coord, b::Coord)::Cdouble
+function proj_xy_dist(a::Coord{Na}, b::Coord{Nb}) where {Na,Nb}
+    @ccall libproj.proj_xy_dist(a::Coord{Na}, b::Coord{Nb})::Cdouble
 end
 
-function proj_xyz_dist(a, b)
-    @ccall libproj.proj_xyz_dist(a::Coord, b::Coord)::Cdouble
+function proj_xyz_dist(a::Coord{Na}, b::Coord{Nb}) where {Na,Nb}
+    @ccall libproj.proj_xyz_dist(a::Coord{Na}, b::Coord{Nb})::Cdouble
 end
 
-function proj_geod(P, a, b)
-    @ccall libproj.proj_geod(P::Ptr{PJ}, a::Coord, b::Coord)::Coord
+function proj_geod(P, a::Coord{N}, b::Coord{N}) where {N}
+    @ccall libproj.proj_geod(P::Ptr{PJ}, a::Coord{N}, b::Coord{N})::Coord{N}
 end
 
 function proj_context_errno(ctx = C_NULL)
@@ -644,8 +651,8 @@ function proj_log_func(app_data, logf, ctx = C_NULL)
     )::Cvoid
 end
 
-function proj_factors(P, lp)
-    @ccall libproj.proj_factors(P::Ptr{PJ}, lp::Coord)::PJ_FACTORS
+function proj_factors(P, lp::Coord{N}) where N
+    @ccall libproj.proj_factors(P::Ptr{PJ}, lp::Coord{N})::PJ_FACTORS
 end
 
 function proj_info()
@@ -1699,12 +1706,12 @@ function proj_list_destroy(result)
     @ccall libproj.proj_list_destroy(result::Ptr{PJ_OBJ_LIST})::Cvoid
 end
 
-function proj_get_suggested_operation(operations, direction, coord, ctx = C_NULL)
+function proj_get_suggested_operation(operations, direction, coord::Coord{N}, ctx = C_NULL) where N
     @ccall libproj.proj_get_suggested_operation(
         ctx::Ptr{PJ_CONTEXT},
         operations::Ptr{PJ_OBJ_LIST},
         direction::PJ_DIRECTION,
-        coord::Coord,
+        coord::Coord{N},
     )::Cint
 end
 

From 1c8bf0f1441eec7dc11c59fa65ce606cd543e967 Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Thu, 27 Apr 2023 22:46:14 +0200
Subject: [PATCH 06/12] dont change proj

---
 src/libproj.jl | 51 ++++++++++++++++++++++----------------------------
 1 file changed, 22 insertions(+), 29 deletions(-)

diff --git a/src/libproj.jl b/src/libproj.jl
index 89bb96a..11caf55 100644
--- a/src/libproj.jl
+++ b/src/libproj.jl
@@ -490,20 +490,20 @@ function proj_degree_output(P, dir)
     @ccall libproj.proj_degree_output(P::Ptr{PJ}, dir::PJ_DIRECTION)::Cint
 end
 
-function proj_trans(P, direction, coord::Coord{N}) where N
-    @ccall libproj.proj_trans(P::Ptr{PJ}, direction::PJ_DIRECTION, coord::Coord{N})::Coord{N}
+function proj_trans(P, direction, coord)
+    @ccall libproj.proj_trans(P::Ptr{PJ}, direction::PJ_DIRECTION, coord::Coord)::Coord
 end
 
 function proj_trans_get_last_used_operation(P)
     @ccall libproj.proj_trans_get_last_used_operation(P::Ptr{PJ})::Ptr{PJ}
 end
 
-function proj_trans_array(P, direction, n, coord::Vector{Coord{N}}) where N
+function proj_trans_array(P, direction, n, coord)
     @ccall libproj.proj_trans_array(
         P::Ptr{PJ},
         direction::PJ_DIRECTION,
         n::Csize_t,
-        coord::Ptr{Coord{N}},
+        coord::Ptr{Coord},
     )::Cint
 end
 
@@ -566,44 +566,37 @@ end
 
 ` Doxygen_Suppress `
 """
-function proj_coord(x, y; time=Inf)
-    z = 0.0 
-    @ccall libproj.proj_coord(x::Cdouble, y::Cdouble, z::Cdouble, t::Cdouble)::Coord{2}
-end
-function proj_coord(x, y, z; time=Inf)
-    @ccall libproj.proj_coord(x::Cdouble, y::Cdouble, z::Cdouble, t::Cdouble)::Coord{3}
-end
-function proj_coord(x, y, z, t)
-    @ccall libproj.proj_coord(x::Cdouble, y::Cdouble, z::Cdouble, t::Cdouble)::Coord{3}
+function proj_coord(x, y, z = 0.0, t = Inf)
+    @ccall libproj.proj_coord(x::Cdouble, y::Cdouble, z::Cdouble, t::Cdouble)::Coord
 end
 
-function proj_roundtrip(P, direction, n, coord::Coord{N}) where N
+function proj_roundtrip(P, direction, n, coord)
     @ccall libproj.proj_roundtrip(
         P::Ptr{PJ},
         direction::PJ_DIRECTION,
         n::Cint,
-        coord::Ptr{Coord{N}},
+        coord::Ptr{Coord},
     )::Cdouble
 end
 
-function proj_lp_dist(P, a::Coord{Na}, b::Coord{Nb}) where {Na,Nb}
-    @ccall libproj.proj_lp_dist(P::Ptr{PJ}, a::Coord{Na}, b::Coord{Nb})::Cdouble
+function proj_lp_dist(P, a, b)
+    @ccall libproj.proj_lp_dist(P::Ptr{PJ}, a::Coord, b::Coord)::Cdouble
 end
 
-function proj_lpz_dist(P, a::Coord{Na}, b::Coord{Nb}) where {Na,Nb}
-    @ccall libproj.proj_lpz_dist(P::Ptr{PJ}, a::Coord{Na}, b::Coord{Nb})::Cdouble
+function proj_lpz_dist(P, a, b)
+    @ccall libproj.proj_lpz_dist(P::Ptr{PJ}, a::Coord, b::Coord)::Cdouble
 end
 
-function proj_xy_dist(a::Coord{Na}, b::Coord{Nb}) where {Na,Nb}
-    @ccall libproj.proj_xy_dist(a::Coord{Na}, b::Coord{Nb})::Cdouble
+function proj_xy_dist(a, b)
+    @ccall libproj.proj_xy_dist(a::Coord, b::Coord)::Cdouble
 end
 
-function proj_xyz_dist(a::Coord{Na}, b::Coord{Nb}) where {Na,Nb}
-    @ccall libproj.proj_xyz_dist(a::Coord{Na}, b::Coord{Nb})::Cdouble
+function proj_xyz_dist(a, b)
+    @ccall libproj.proj_xyz_dist(a::Coord, b::Coord)::Cdouble
 end
 
-function proj_geod(P, a::Coord{N}, b::Coord{N}) where {N}
-    @ccall libproj.proj_geod(P::Ptr{PJ}, a::Coord{N}, b::Coord{N})::Coord{N}
+function proj_geod(P, a, b)
+    @ccall libproj.proj_geod(P::Ptr{PJ}, a::Coord, b::Coord)::Coord
 end
 
 function proj_context_errno(ctx = C_NULL)
@@ -651,8 +644,8 @@ function proj_log_func(app_data, logf, ctx = C_NULL)
     )::Cvoid
 end
 
-function proj_factors(P, lp::Coord{N}) where N
-    @ccall libproj.proj_factors(P::Ptr{PJ}, lp::Coord{N})::PJ_FACTORS
+function proj_factors(P, lp)
+    @ccall libproj.proj_factors(P::Ptr{PJ}, lp::Coord)::PJ_FACTORS
 end
 
 function proj_info()
@@ -1706,12 +1699,12 @@ function proj_list_destroy(result)
     @ccall libproj.proj_list_destroy(result::Ptr{PJ_OBJ_LIST})::Cvoid
 end
 
-function proj_get_suggested_operation(operations, direction, coord::Coord{N}, ctx = C_NULL) where N
+function proj_get_suggested_operation(operations, direction, coord, ctx = C_NULL)
     @ccall libproj.proj_get_suggested_operation(
         ctx::Ptr{PJ_CONTEXT},
         operations::Ptr{PJ_OBJ_LIST},
         direction::PJ_DIRECTION,
-        coord::Coord{N},
+        coord::Coord,
     )::Cint
 end
 

From c65bdc91c5960125c14b8731581a69faea1ad4ca Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Thu, 27 Apr 2023 22:46:43 +0200
Subject: [PATCH 07/12] use Point rather than Coord

---
 src/geointerface.jl | 60 ++++++++++++++++++++++++++++++++-------------
 1 file changed, 43 insertions(+), 17 deletions(-)

diff --git a/src/geointerface.jl b/src/geointerface.jl
index 3ce2877..809367a 100644
--- a/src/geointerface.jl
+++ b/src/geointerface.jl
@@ -2,23 +2,49 @@ using GeoInterface
 
 const GI = GeoInterface
 
-# Coord is a PointTrait
-GI.isgeometry(::Type{Coord}) = true
-GI.geomtrait(::Coord) = GI.PointTrait()
-GI.x(::PointTrait, c::Coord) = c.x # TODO: these may have the wrong order for some CRS. 
-GI.y(::PointTrait, c::Coord) = c.y 
-GI.z(::PointTrait, c::Coord{3}) = c.z
-GI.is3d(::PointTrait, c::Coord{2}) = false
-GI.is3d(::PointTrait, c::Coord{3}) = true
-GI.ismeasured(::PointTrait, c::Coord) = false
-function GI.getcoord(::PointTrait, c::Coord{N}, i::Int) where N
+struct Point{D}
+    c::Coord
+end
+Point{D}(x, y, z, t) where N = Point{N}(Coord(x, y, z, y))
+
+Point{3}(x, y, z=0.0; time=Inf) = Point{3}(x, y, z, time)
+Point{2}(x, y; time=Inf) = Point{2}(x, y, 0.0, time)
+Point(x, y; time=Inf) = Point{2}(x, y, 0.0, time)
+Point(x, y, z; time=Inf) = Point{3}(x, y, z, time)
+Point(x, y, z, t) = Point{3}(x, y, z, t)
+
+# this shields a StackOverflow from the splatting constructor
+Point(::Real) = error("Proj.Point takes 2 to 4 numbers, one given")
+Point(v) = Point(v...)
+function Point(v::AbstractVector{<:Real})
+    n = length(v)
+    return if n == 2
+        Point{2}(v[begin], v[begin+1])
+    elseif n == 3
+        Point{3}(v[begin], v[begin+1], v[begin+2])
+    elseif n == 4
+        Point{3}(v[begin], v[begin+1], v[begin+2], v[begin+3])
+    else
+        error("Proj.Point takes 2 to 4 numbers")
+    end
+end
+
+GI.isgeometry(::Type{Point}) = true
+GI.geomtrait(::Point) = GI.PointTrait()
+GI.x(::PointTrait, c::Point) = c.x # TODO: these may have the wrong order for some CRS. 
+GI.y(::PointTrait, c::Point) = c.y 
+GI.z(::PointTrait, c::Point{3}) = c.z
+GI.is3d(::PointTrait, c::Point{2}) = false
+GI.is3d(::PointTrait, c::Point{3}) = true
+GI.ismeasured(::PointTrait, c::Point) = false
+function GI.getcoord(::PointTrait, c::Point{N}, i::Int) where N
     checkbounds(1:N, i)
     @inbounds c[i]
 end
-GI.getcoord(::PointTrait, c::Coord{2}) = (c.x, c.y)
-GI.getcoord(::PointTrait, c::Coord{3}) = (c.x, c.y, c.z)
-coordnames(::PointTrait, ::Coord{2}) = (:X, :Y)
-coordnames(::PointTrait, ::Coord{3}) = (:X, :Y, :Z)
+GI.getcoord(::PointTrait, c::Point{2}) = (c.x, c.y)
+GI.getcoord(::PointTrait, c::Point{3}) = (c.x, c.y, c.z)
+coordnames(::PointTrait, ::Point{2}) = (:X, :Y)
+coordnames(::PointTrait, ::Point{3}) = (:X, :Y, :Z)
 
 """
     reproject(geometry, source_crs, target_crs)
@@ -27,7 +53,7 @@ coordnames(::PointTrait, ::Coord{3}) = (:X, :Y, :Z)
 Reproject any GeoInterface.jl compatible `geometry` from `source_crs` to `target_crs`.
 
 The returned object will be constructed from `GeoInterface.WrapperGeometry`
-geometries, wrapping views of `Proj.Coord`.
+geometries, wrapping views of `Proj.Point`.
 """
 function reproject(geom; source_crs=nothing, target_crs, kw...)
     source_crs = isnothing(source_crs) ? GeoInterface.crs(geom) : source_crs
@@ -37,11 +63,11 @@ end
 function reproject(geom, source_crs, target_crs; time=Inf, always_xy=true)
     wrapped_geom, coords = if GI.is3d(geom)
         reconstruct(geom; crs=target_crs) do p
-            Proj.Coord{3}(GI.x(p), GI.y(p), GI.z(p); time)
+            Proj.Point{3}(GI.x(p), GI.y(p), GI.z(p); time)
         end
     else
         reconstruct(geom; crs=target_crs) do p
-            Proj.Coord{2}(GI.x(p), GI.y(p); time)
+            Proj.Point{2}(GI.x(p), GI.y(p); time)
         end
     end
     source_crs1 = convert(Proj.CRS, source_crs)

From be9849bbf449067c530efc68aadaf132604fdb9d Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Thu, 27 Apr 2023 22:53:29 +0200
Subject: [PATCH 08/12] refactor

---
 src/Proj.jl          |  8 ++------
 src/geointerface.jl  | 20 ++++++++++----------
 test/geointerface.jl |  1 -
 3 files changed, 12 insertions(+), 17 deletions(-)

diff --git a/src/Proj.jl b/src/Proj.jl
index 91e4a7a..7dfbfc6 100644
--- a/src/Proj.jl
+++ b/src/Proj.jl
@@ -50,17 +50,13 @@ julia> Proj.Coord((1.0, 2.0, 3.0))
  Inf
 ```
 """
-struct Coord{N} <: AbstractVector{Float64}
+struct Coord <: AbstractVector{Float64}
     x::Float64
     y::Float64
     z::Float64
     t::Float64
+    Coord(x, y, z = 0.0, t = Inf) = new(x, y, z, t)
 end
-Coord{3}(x, y, z=0.0; time=Inf) = Coord{3}(x, y, z, time)
-Coord{2}(x, y; time=Inf) = Coord{2}(x, y, 0.0, time)
-Coord(x, y; time=Inf) = Coord{2}(x, y, 0.0, time)
-Coord(x, y, z; time=Inf) = Coord{3}(x, y, z, time)
-Coord(x, y, z, t) = Coord{3}(x, y, z, t)
 
 # this shields a StackOverflow from the splatting constructor
 Coord(::Real) = error("Proj.Coord takes 2 to 4 numbers, one given")
diff --git a/src/geointerface.jl b/src/geointerface.jl
index 809367a..c465d71 100644
--- a/src/geointerface.jl
+++ b/src/geointerface.jl
@@ -3,9 +3,9 @@ using GeoInterface
 const GI = GeoInterface
 
 struct Point{D}
-    c::Coord
+    coord::Coord
 end
-Point{D}(x, y, z, t) where N = Point{N}(Coord(x, y, z, y))
+Point{D}(x, y, z, t) where D = Point{D}(Coord(x, y, z, y))
 
 Point{3}(x, y, z=0.0; time=Inf) = Point{3}(x, y, z, time)
 Point{2}(x, y; time=Inf) = Point{2}(x, y, 0.0, time)
@@ -29,20 +29,20 @@ function Point(v::AbstractVector{<:Real})
     end
 end
 
-GI.isgeometry(::Type{Point}) = true
+GI.isgeometry(::Type{<:Point}) = true
 GI.geomtrait(::Point) = GI.PointTrait()
-GI.x(::PointTrait, c::Point) = c.x # TODO: these may have the wrong order for some CRS. 
-GI.y(::PointTrait, c::Point) = c.y 
-GI.z(::PointTrait, c::Point{3}) = c.z
+GI.x(::PointTrait, c::Point) = c.coord.x # TODO: these may have the wrong order for some CRS. 
+GI.y(::PointTrait, c::Point) = c.coord.y 
+GI.z(::PointTrait, c::Point{3}) = c.coord.z
 GI.is3d(::PointTrait, c::Point{2}) = false
 GI.is3d(::PointTrait, c::Point{3}) = true
 GI.ismeasured(::PointTrait, c::Point) = false
-function GI.getcoord(::PointTrait, c::Point{N}, i::Int) where N
+function GI.getcoord(::PointTrait, p::Point{N}, i::Int) where N
     checkbounds(1:N, i)
-    @inbounds c[i]
+    @inbounds p.coord[i]
 end
-GI.getcoord(::PointTrait, c::Point{2}) = (c.x, c.y)
-GI.getcoord(::PointTrait, c::Point{3}) = (c.x, c.y, c.z)
+GI.getcoord(::PointTrait, p::Point{2}) = (c.x, c.y)
+GI.getcoord(::PointTrait, p::Point{3}) = (c.x, c.y, c.z)
 coordnames(::PointTrait, ::Point{2}) = (:X, :Y)
 coordnames(::PointTrait, ::Point{3}) = (:X, :Y, :Z)
 
diff --git a/test/geointerface.jl b/test/geointerface.jl
index 762e5ee..118c257 100644
--- a/test/geointerface.jl
+++ b/test/geointerface.jl
@@ -3,7 +3,6 @@ using Test
 using GeoInterface
 using GeoFormatTypes
 const GI = GeoInterface
-using Plots
 
 using GeoInterfaceRecipes
 GeoInterfaceRecipes.@enable_geo_plots GeoInterface.Wrappers.WrapperGeometry

From 9b4e3c86759ecd5b336b69baed3e4e7657fc6a7a Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Thu, 27 Apr 2023 22:56:46 +0200
Subject: [PATCH 09/12] test properly

---
 test/geointerface.jl | 22 +++++++++++++++++++---
 1 file changed, 19 insertions(+), 3 deletions(-)

diff --git a/test/geointerface.jl b/test/geointerface.jl
index 118c257..106c9c7 100644
--- a/test/geointerface.jl
+++ b/test/geointerface.jl
@@ -11,16 +11,32 @@ ring1 = GI.LinearRing([(1, 2), (7, 4), (5, 6), (1, 2)])
 ring2 = GI.LinearRing([(11, 2), (20, 4), (15, 6), (11, 2)])
 hole2 = GI.LinearRing([(14, 4), (16, 4), (17, 5), (14, 4)])
 
+ref_points3857 = map(GeoInterface.getpoint(multipolygon)) do p
+    trans([GeoInterface.x(p), GeoInterface.y(p)])
+end
+
 polygon1 = GI.Polygon([ring1])
 polygon2 = GI.Polygon([ring2, hole2])
 multipolygon = GI.MultiPolygon([polygon1, polygon2])
 
 multipolygon3857 = Proj.reproject(multipolygon, EPSG(4326), EPSG(3857))
 multipolygon4326 = Proj.reproject(multipolygon3857; target_crs=EPSG(4326))
-points1 = collect(GI.getpoint(multipolygon))
-points2 = GI.convert.(Tuple, GI.getpoint(multipolygon4326))
+points4326_1 = collect(GI.getpoint(multipolygon))
+points4326_2 = GI.convert.(Tuple, GI.getpoint(multipolygon4326))
+points3857 = GI.convert.(Tuple, GI.getpoint(multipolygon3857))
+
+# Comparison to regular `trans` on points
+@test all(map((p1, p2) -> all(map(isapprox, p1, p2)), ref_points3857, points3857))
+
 # Round trip comparison
-@test all(map((p1, p2) -> all(map(isapprox, p1, p2)), points1, points2))
+@test all(map((p1, p2) -> all(map(isapprox, p1, p2)), points4326_1, points4326_2))
+
+# Embedded crs check
 @test GI.crs(multipolygon3857) == EPSG(3857)
 @test GI.crs(multipolygon4326) == EPSG(4326)
 
+source_crs = convert(Proj.CRS, EPSG(4326))
+target_crs = convert(Proj.CRS, EPSG(3857))
+trans = Proj.Transformation(source_crs, target_crs; always_xy=true)
+
+

From e1d8de9ae08d8cfbfcc95cc215cb7ec20d4f9909 Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Thu, 27 Apr 2023 23:07:32 +0200
Subject: [PATCH 10/12] allow passing in transform

---
 src/geointerface.jl | 24 ++++++++++++++++++------
 1 file changed, 18 insertions(+), 6 deletions(-)

diff --git a/src/geointerface.jl b/src/geointerface.jl
index c465d71..f153bef 100644
--- a/src/geointerface.jl
+++ b/src/geointerface.jl
@@ -55,12 +55,25 @@ Reproject any GeoInterface.jl compatible `geometry` from `source_crs` to `target
 The returned object will be constructed from `GeoInterface.WrapperGeometry`
 geometries, wrapping views of `Proj.Point`.
 """
-function reproject(geom; source_crs=nothing, target_crs, kw...)
-    source_crs = isnothing(source_crs) ? GeoInterface.crs(geom) : source_crs
-    isnothing(source_crs) && throw(ArgumentError("geom has no crs attatched. Pass a `source_crs` keyword"))
-    reproject(geom, source_crs, target_crs; kw...)
+function reproject(geom; 
+    source_crs=nothing, target_crs=nothing, transform=nothing, kw...
+)
+    if isnothing(transform)
+        source_crs = isnothing(source_crs) ? GeoInterface.crs(geom) : source_crs
+        isnothing(source_crs) && throw(ArgumentError("geom has no crs attatched. Pass a `source_crs` keyword"))
+        reproject(geom, source_crs, target_crs; transform, kw...)
+    else
+        reproject(geom, transform)
+    end
+end
+function reproject(geom, source_crs, target_crs; 
+    time=Inf, 
+    always_xy=true,
+    transform=Proj.Transformation(source_crs1, target_crs1; always_xy),
+)
+    reproject(geom, transform; time)
 end
-function reproject(geom, source_crs, target_crs; time=Inf, always_xy=true)
+function reproject(geom, transform::Transformation; time)
     wrapped_geom, coords = if GI.is3d(geom)
         reconstruct(geom; crs=target_crs) do p
             Proj.Point{3}(GI.x(p), GI.y(p), GI.z(p); time)
@@ -72,7 +85,6 @@ function reproject(geom, source_crs, target_crs; time=Inf, always_xy=true)
     end
     source_crs1 = convert(Proj.CRS, source_crs)
     target_crs1 = convert(Proj.CRS, target_crs)
-    trans = Proj.Transformation(source_crs1, target_crs1; always_xy)
     err = Proj.proj_trans_array(trans.pj, Proj.PJ_FWD, length(coords), coords)
     err == 0 || error("Proj error $err")
     # crs = target_crs isa GFT.GeoFormat ? target_crs : convert(WellKnownText, crs)

From 020c54dc8a5b28233af1b8549963976a3bfed6ef Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Fri, 28 Apr 2023 19:45:40 +0200
Subject: [PATCH 11/12] actually run the tests

---
 test/geointerface.jl | 10 ++++------
 test/runtests.jl     |  1 +
 2 files changed, 5 insertions(+), 6 deletions(-)

diff --git a/test/geointerface.jl b/test/geointerface.jl
index 106c9c7..68e8616 100644
--- a/test/geointerface.jl
+++ b/test/geointerface.jl
@@ -11,6 +11,10 @@ ring1 = GI.LinearRing([(1, 2), (7, 4), (5, 6), (1, 2)])
 ring2 = GI.LinearRing([(11, 2), (20, 4), (15, 6), (11, 2)])
 hole2 = GI.LinearRing([(14, 4), (16, 4), (17, 5), (14, 4)])
 
+# Set up a regular tranformation of the points for reference
+source_crs = convert(Proj.CRS, EPSG(4326))
+target_crs = convert(Proj.CRS, EPSG(3857))
+trans = Proj.Transformation(source_crs, target_crs; always_xy=true)
 ref_points3857 = map(GeoInterface.getpoint(multipolygon)) do p
     trans([GeoInterface.x(p), GeoInterface.y(p)])
 end
@@ -34,9 +38,3 @@ points3857 = GI.convert.(Tuple, GI.getpoint(multipolygon3857))
 # Embedded crs check
 @test GI.crs(multipolygon3857) == EPSG(3857)
 @test GI.crs(multipolygon4326) == EPSG(4326)
-
-source_crs = convert(Proj.CRS, EPSG(4326))
-target_crs = convert(Proj.CRS, EPSG(3857))
-trans = Proj.Transformation(source_crs, target_crs; always_xy=true)
-
-
diff --git a/test/runtests.jl b/test/runtests.jl
index 1597bff..28d01bc 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -8,5 +8,6 @@ is_approx(a, b) = all(isapprox(c[1], c[2]) for c in zip(a, b))
 @testset "Proj" begin
     include("libproj.jl")
     include("applications.jl")
+    include("geointerface.jl")
 
 end # testset "Proj"

From c1ed76cac6df0fddbb16f61fd478b95df5f2f366 Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Fri, 28 Apr 2023 19:57:53 +0200
Subject: [PATCH 12/12] fix docs

---
 src/geointerface.jl | 44 ++++++++++++++++++++++++++++----------------
 1 file changed, 28 insertions(+), 16 deletions(-)

diff --git a/src/geointerface.jl b/src/geointerface.jl
index f153bef..49b461a 100644
--- a/src/geointerface.jl
+++ b/src/geointerface.jl
@@ -47,13 +47,30 @@ coordnames(::PointTrait, ::Point{2}) = (:X, :Y)
 coordnames(::PointTrait, ::Point{3}) = (:X, :Y, :Z)
 
 """
-    reproject(geometry, source_crs, target_crs)
-    reproject(geometry; [source_crs,] target_crs)
+    reproject(geometry; source_crs, target_crs, transform, always_xy, time)
+    reproject(geometry, source_crs, target_crs; always_xy, time)
+    reproject(geometry, transform; always_xy, time)
 
 Reproject any GeoInterface.jl compatible `geometry` from `source_crs` to `target_crs`.
 
 The returned object will be constructed from `GeoInterface.WrapperGeometry`
-geometries, wrapping views of `Proj.Point`.
+geometries, wrapping views of a `Vector{Proj.Point{D}}`, where `D` is the dimension.
+
+# Arguments
+
+- `geometry`: Any GeoInterface.jl compatible geometries.
+- `source_crs`: the source coordinate referece system, as a GeoFormatTypes.jl object or a string.
+- `target_crs`: the target coordinate referece system, as a GeoFormatTypes.jl object or a string.
+
+If these a passed as keywords, `transform` will take priority. 
+Without it `target_crs` is always needed, and `source_crs` is
+needed if it is not retreivable from the geometry with `GeoInterface.crs(geometry)`.
+
+# Keywords
+
+-`always_xy`: force x, y coordinate order, `true` by default. 
+    `false` will expect and return points in the crs coordinate order.
+-`time`: the time for the coordinates. `Inf` by default.
 """
 function reproject(geom; 
     source_crs=nothing, target_crs=nothing, transform=nothing, kw...
@@ -92,22 +109,17 @@ function reproject(geom, transform::Transformation; time)
 end
 
 
-"""
-    reconstruct(f::Function, geom)
-    reconstruct(geometry, points::AbstractVector)
-
-Reconstruct GeoInterface compatible `geometry` from the `PointTrait` geometries
-in `points`, or from the result of function `f` over the points in `geom`.
+    # reconstruct(f::Function, geom)
 
-The returned object will be constructed from `GeoInterface.WrapperGeometry`
-geometries, wrapping views into the `points` vector. Otherwise, the structure and
-GeoInterface traits will be the same as for `geometry`.
+# Reconstruct GeoInterface compatible `geometry` from the `PointTrait` geometries
+# in `points`, or from the result of function `f` over the points in `geom`.
 
-This is a lazy operation, only allocating for outer wrappers of nested geometries.
-Later changes to the `points` vector *will* affect the returned geometry.
+# The returned object will be constructed from `GeoInterface.WrapperGeometry`
+# geometries, wrapping views into the `points` vector. Otherwise, the structure and
+# GeoInterface traits will be the same as for `geometry`.
 
-`GeoInterface.npoint(geometry) == length(points)` must be `true`.
-"""
+# This is a lazy operation, only allocating for outer wrappers of nested geometries.
+# Later changes to the `points` vector *will* affect the returned geometry.
 function reconstruct(f, geom; crs=nothing)
     p1 = f(first(GI.getpoint(geom)))
     T = typeof(p1)