Add support for Extent

Project.toml

@@ -5,6 +5,7 @@ version = "1.9.0"
 [deps]
 CEnum = "fa961155-64e5-5f13-b03f-caf6b980ea82"
 CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
+Extents = "411431e0-e8b7-467b-b5e0-f676ba4f2910"
 GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 NetworkOptions = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
@@ -14,6 +15,7 @@ PROJ_jll = "58948b4f-47e0-5654-a9ad-f609743f8632"
 Aqua = "0.8"
 CEnum = "0.2, 0.3, 0.4, 0.5"
 CoordinateTransformations = "0.6"
+Extents = "0.1"
 GeoFormatTypes = "0.4"
 GeoInterface = "1.3"
 NetworkOptions = "1"

src/Proj.jl

@@ -4,6 +4,7 @@ using PROJ_jll
 using CEnum
 using CoordinateTransformations
 using NetworkOptions: ca_roots
+import Extents
 import GeoFormatTypes as GFT
 import GeoInterface as GI
 

src/coord.jl

@@ -262,6 +262,64 @@ function (trans::Transformation)(coord)
     end
 end
 
+"""
+    (trans::Transformation)(extent::Extents.Extent) -> Extents.Extent
+
+Transform an `Extents.Extent` using the transformation. The extent must have X and Y
+dimensions. Extra dimensions (Z, Ti, etc.) are preserved unchanged.
+
+Throws an `ArgumentError` if the transformation is rotated or non-axis-aligned,
+since an axis-aligned bounding box cannot accurately represent a rotated extent.
+
+# Examples
+```julia
+julia> trans = Proj.Transformation("EPSG:4326", "EPSG:28992", always_xy=true)
+julia> extent = Extents.Extent(X=(5.0, 6.0), Y=(52.0, 53.0))
+julia> trans(extent)
+Extent(X = (131793.0, 199684.0), Y = (449093.0, 560431.0))
+```
+"""
+function (trans::Transformation)(extent::Extents.Extent)
+    haskey(extent, :X) && haskey(extent, :Y) ||
+        throw(ArgumentError("Extent must have X and Y dimensions"))
+
+    _check_axis_aligned(trans, extent) ||
+        throw(ArgumentError("Cannot transform Extent: CRS is rotated or non-axis-aligned"))
+
+    (xmin, xmax), (ymin, ymax) = bounds(trans, extent.X, extent.Y)
+
+    return Extents.Extent(merge(NamedTuple(extent), (X=(xmin, xmax), Y=(ymin, ymax))))
+end
+
+# Check if corners transform to an axis-aligned rectangle
+function _check_axis_aligned(trans::Transformation, extent::Extents.Extent; rtol=0.1)
+    xmin, xmax = extent.X
+    ymin, ymax = extent.Y
+
+    # Transform the 4 corners
+    c1 = trans(xmin, ymin)  # bottom-left
+    c2 = trans(xmax, ymin)  # bottom-right
+    c3 = trans(xmin, ymax)  # top-left
+    c4 = trans(xmax, ymax)  # top-right
+
+    # For axis-aligned transform, vertical edges should have similar x-coords
+    # and horizontal edges should have similar y-coords
+    # Use relative tolerance to allow for projection curvature
+    x_range = max(abs(c1[1] - c2[1]), abs(c3[1] - c4[1]), 1.0)
+    y_range = max(abs(c1[2] - c3[2]), abs(c2[2] - c4[2]), 1.0)
+
+    # Check left edge: c1 and c3 should have similar x
+    abs(c1[1] - c3[1]) < rtol * x_range || return false
+    # Check right edge: c2 and c4 should have similar x
+    abs(c2[1] - c4[1]) < rtol * x_range || return false
+    # Check bottom edge: c1 and c2 should have similar y
+    abs(c1[2] - c2[2]) < rtol * y_range || return false
+    # Check top edge: c3 and c4 should have similar y
+    abs(c3[2] - c4[2]) < rtol * y_range || return false
+
+    return true
+end
+
 """
     Proj.bounds(trans::Transformation, (xmin, xmax), (ymin,ymax); densify_pts=21) -> ((bxmin, bxmax), (bymin, bymax))
 

test/libproj.jl

@@ -1,6 +1,7 @@
 using Test
 using StaticArrays
 using Proj
+import Extents
 import GeoInterface as GI
 import PROJ_jll
 import GeoFormatTypes as GFT
@@ -253,6 +254,46 @@ end
 
 end
 
+@testset "Extent transformation" begin
+    trans = Proj.Transformation("EPSG:4326", "EPSG:28992", always_xy=true)
+
+    # Basic extent transformation
+    ext = Extents.Extent(X=(5.0, 6.0), Y=(52.0, 53.0))
+    result = trans(ext)
+    @test result isa Extents.Extent
+    @test haskey(result, :X)
+    @test haskey(result, :Y)
+    @test result.X[1] < result.X[2]
+    @test result.Y[1] < result.Y[2]
+
+    # Test that a point inside the extent transforms to inside the result
+    px, py = trans(5.5, 52.5)
+    @test result.X[1] < px < result.X[2]
+    @test result.Y[1] < py < result.Y[2]
+
+    # Extra dimensions are preserved
+    ext_with_z = Extents.Extent(X=(5.0, 6.0), Y=(52.0, 53.0), Z=(0.0, 100.0))
+    result_z = trans(ext_with_z)
+    @test haskey(result_z, :Z)
+    @test result_z.Z == (0.0, 100.0)
+
+    # Order of dimensions is preserved
+    ext_ordered = Extents.Extent(Z=(0.0, 100.0), X=(5.0, 6.0), Y=(52.0, 53.0))
+    result_ordered = trans(ext_ordered)
+    @test keys(result_ordered) == (:Z, :X, :Y)
+
+    # Error without X and Y
+    @test_throws ArgumentError trans(Extents.Extent(Z=(0.0, 100.0)))
+
+    # Error on rotated CRS (oblique mercator)
+    trans_rotated = Proj.Transformation(
+        "EPSG:4326",
+        "+proj=omerc +lat_0=45 +lonc=0 +alpha=45 +gamma=45 +k=1 +x_0=0 +y_0=0 +ellps=WGS84",
+        always_xy=true
+    )
+    @test_throws ArgumentError trans_rotated(ext)
+end
+
 @testset "dense 4D coord vector transformation" begin
     source_crs = Proj.proj_create("EPSG:4326")
     target_crs = Proj.proj_create("EPSG:28992")