Add boundary option to Proj and Morphological tranforms (#1117)

* Add 'boundary' option to 'Proj' and 'Morphological' tranforms

* Update docstring

* Fix typo

* Apply suggestions from code review

Co-authored-by: Júlio Hoffimann <julio.hoffimann@gmail.com>

---------

Co-authored-by: Júlio Hoffimann <julio.hoffimann@gmail.com>

src/transforms/morphological.jl

@@ -3,24 +3,34 @@
 # ------------------------------------------------------------------
 
 """
-    Morphological(fun)
+    Morphological(fun; boundary=false)
 
 Morphological transform given by a function `fun`
 that maps the coordinates of a geometry or a domain
 to new coordinates (`coords -> newcoords`).
 
+Optionally, the transform samples the `boundary` of
+polytopes, if this option is `true`, to handle distortions
+that occur in manifold conversions.
+
 # Examples
 
 ```julia
 ball = Ball((0, 0), 1)
 ball |> Morphological(c -> Cartesian(c.x + c.y, c.y, c.x - c.y))
+
+triangle = Triangle(latlon(0, 0), latlon(0, 45), latlon(45, 0))
+triangle |> Morphological(c -> LatLonAlt(c.lat, c.lon, 0.0m), boundary=true)
 ```
 """
-struct Morphological{F<:Function} <: CoordinateTransform
+struct Morphological{Boundary,F<:Function} <: CoordinateTransform
   fun::F
+  Morphological{Boundary}(fun::F) where {Boundary,F<:Function} = new{Boundary,F}(fun)
 end
 
-parameters(t::Morphological) = (; fun=t.fun)
+Morphological(fun; boundary=false) = Morphological{boundary}(fun)
+
+parameters(t::Morphological{Boundary}) where {Boundary} = (fun=t.fun, boundary=Boundary)
 
 applycoord(t::Morphological, p::Point) = Point(t.fun(coords(p)))
 
@@ -30,13 +40,26 @@ applycoord(::Morphological, v::Vec) = v
 # SPECIAL CASES
 # --------------
 
-applycoord(t::Morphological, g::Geometry) = TransformedGeometry(g, t)
+applycoord(t::Morphological, g::Primitive) = TransformedGeometry(g, t)
 
-# method to fix ambiguities
-applycoord(t::Morphological, g::TransformedGeometry) = TransformedGeometry(g, t)
+applycoord(t::Morphological{true}, g::Polytope) = TransformedGeometry(g, t)
 
 applycoord(t::Morphological, g::RegularGrid) = TransformedGrid(g, t)
 
 applycoord(t::Morphological, g::RectilinearGrid) = TransformedGrid(g, t)
 
 applycoord(t::Morphological, g::StructuredGrid) = TransformedGrid(g, t)
+
+# -----------
+# IO METHODS
+# -----------
+
+Base.show(io::IO, t::Morphological{Boundary}) where {Boundary} =
+  print(io, "Morphological(fun: $(t.fun), boundary: $Boundary)")
+
+function Base.show(io::IO, ::MIME"text/plain", t::Morphological{Boundary}) where {Boundary}
+  summary(io, t)
+  println(io)
+  println(io, "├─ fun: $(t.fun)")
+  print(io, "└─ boundary: $Boundary")
+end

src/transforms/proj.jl

@@ -3,12 +3,16 @@
 # ------------------------------------------------------------------
 
 """
-    Proj(CRS)
-    Proj(code)
+    Proj(CRS; boundary=false)
+    Proj(code; boundary=false)
 
 Convert the coordinates of geometry or domain to a given
 coordinate reference system `CRS` or EPSG/ESRI `code`.
 
+Optionally, the transform samples the `boundary` of
+polytopes, if this option is `true`, to handle distortions
+that occur in manifold conversions.
+
 ## Examples
 
 ```julia
@@ -17,17 +21,18 @@ Proj(WebMercator)
 Proj(Mercator{WGS84Latest})
 Proj(EPSG{3395})
 Proj(ESRI{54017})
+Proj(Robinson, boundary=true)
 ```
 """
-struct Proj{CRS} <: CoordinateTransform end
+struct Proj{CRS,Boundary} <: CoordinateTransform end
 
-Proj(CRS) = Proj{CRS}()
+Proj(CRS; boundary=false) = Proj{CRS,boundary}()
 
-Proj(code::Type{<:EPSG}) = Proj{CoordRefSystems.get(code)}()
+Proj(code::Type{<:EPSG}; kwargs...) = Proj(CoordRefSystems.get(code); kwargs...)
 
-Proj(code::Type{<:ESRI}) = Proj{CoordRefSystems.get(code)}()
+Proj(code::Type{<:ESRI}; kwargs...) = Proj(CoordRefSystems.get(code); kwargs...)
 
-parameters(::Proj{CRS}) where {CRS} = (; CRS)
+parameters(::Proj{CRS,Boundary}) where {CRS,Boundary} = (CRS=CRS, boundary=Boundary)
 
 # avoid constructing a new geometry or domain when the CRS is the same
 function apply(t::Proj{CRS}, g::GeometryOrDomain) where {CRS}
@@ -51,13 +56,13 @@ applycoord(::Proj, v::Vec) = v
 # SPECIAL CASES
 # --------------
 
-applycoord(t::Proj{<:Projected}, g::Geometry{<:🌐}) = TransformedGeometry(g, t)
+applycoord(t::Proj{<:Projected}, g::Primitive{<:🌐}) = TransformedGeometry(g, t)
+
+applycoord(t::Proj{<:Geographic}, g::Primitive{<:𝔼}) = TransformedGeometry(g, t)
 
-applycoord(t::Proj{<:Geographic}, g::Geometry{<:𝔼}) = TransformedGeometry(g, t)
+applycoord(t::Proj{<:Projected,true}, g::Polytope{K,<:🌐}) where {K} = TransformedGeometry(g, t)
 
-# methods to fix ambiguities
-applycoord(t::Proj{<:Projected}, g::TransformedGeometry{<:🌐}) = TransformedGeometry(g, t)
-applycoord(t::Proj{<:Geographic}, g::TransformedGeometry{<:𝔼}) = TransformedGeometry(g, t)
+applycoord(t::Proj{<:Geographic,true}, g::Polytope{K,<:𝔼}) where {K} = TransformedGeometry(g, t)
 
 applycoord(t::Proj, g::RegularGrid) = TransformedGrid(g, t)
 
@@ -69,12 +74,13 @@ applycoord(t::Proj, g::StructuredGrid) = TransformedGrid(g, t)
 # IO METHODS
 # -----------
 
-Base.show(io::IO, ::Proj{CRS}) where {CRS} = print(io, "Proj(CRS: $CRS)")
+Base.show(io::IO, ::Proj{CRS,Boundary}) where {CRS,Boundary} = print(io, "Proj(CRS: $CRS, boundary: $Boundary)")
 
-function Base.show(io::IO, ::MIME"text/plain", t::Proj{CRS}) where {CRS}
+function Base.show(io::IO, ::MIME"text/plain", t::Proj{CRS,Boundary}) where {CRS,Boundary}
   summary(io, t)
   println(io)
-  print(io, "└─ CRS: $CRS")
+  println(io, "├─ CRS: $CRS")
+  print(io, "└─ boundary: $Boundary")
 end
 
 # -----------------

test/sets.jl

@@ -30,22 +30,24 @@
   q = Quadrangle(latlon(0, 0), latlon(0, 1), latlon(1, 1), latlon(1, 0)) |> Proj(WebMercator)
   geoms = [s, t, q]
   gset = GeometrySet(geoms)
+  @test eltype(gset) <: Polytope
   @test crs(gset) <: LatLon
   gset = GeometrySet(g for g in geoms)
+  @test eltype(gset) <: Polytope
   @test crs(gset) <: LatLon
   geoms = [t, s, q]
   gset = GeometrySet(geoms)
-  @test eltype(gset) <: TransformedGeometry
+  @test eltype(gset) <: Polytope
   @test crs(gset) <: PlateCarree
   gset = GeometrySet(g for g in geoms)
-  @test eltype(gset) <: TransformedGeometry
+  @test eltype(gset) <: Polytope
   @test crs(gset) <: PlateCarree
   geoms = [q, s, t]
   gset = GeometrySet(geoms)
-  @test eltype(gset) <: TransformedGeometry
+  @test eltype(gset) <: Polytope
   @test crs(gset) <: WebMercator
   gset = GeometrySet(g for g in geoms)
-  @test eltype(gset) <: TransformedGeometry
+  @test eltype(gset) <: Polytope
   @test crs(gset) <: WebMercator
 
   # conversion

test/transforms.jl

@@ -1251,14 +1251,15 @@ end
   @test !isaffine(Proj(Polar))
   @test !TB.isrevertible(Proj(Polar))
   @test !TB.isinvertible(Proj(Polar))
-  @test TB.parameters(Proj(Polar)) == (; CRS=Polar)
-  @test TB.parameters(Proj(EPSG{3395})) == (; CRS=Mercator{WGS84Latest})
-  @test TB.parameters(Proj(ESRI{54017})) == (; CRS=Behrmann{WGS84Latest})
+  @test TB.parameters(Proj(Polar)) == (CRS=Polar, boundary=false)
+  @test TB.parameters(Proj(EPSG{3395})) == (CRS=Mercator{WGS84Latest}, boundary=false)
+  @test TB.parameters(Proj(ESRI{54017})) == (CRS=Behrmann{WGS84Latest}, boundary=false)
   f = Proj(Mercator)
-  @test sprint(show, f) == "Proj(CRS: CoordRefSystems.Mercator)"
+  @test sprint(show, f) == "Proj(CRS: CoordRefSystems.Mercator, boundary: false)"
   @test sprint(show, MIME"text/plain"(), f) == """
   Proj transform
-  └─ CRS: CoordRefSystems.Mercator"""
+  ├─ CRS: CoordRefSystems.Mercator
+  └─ boundary: false"""
 
   # ----
   # VEC
@@ -1457,14 +1458,39 @@ end
   f = Proj(crs(merc(0, 0)))
   r, c = TB.apply(f, d)
   @test r === d
+
+  # ----------------
+  # BOUNDARY OPTION
+  # ----------------
+
+  f = Proj(Mercator, boundary=false)
+  g = Triangle(latlon(0, 0), latlon(0, 45), latlon(45, 0))
+  r, c = TB.apply(f, g)
+  @test r isa Triangle
+  f = Proj(Mercator, boundary=true)
+  r, c = TB.apply(f, g)
+  @test r isa TransformedGeometry
+
+  f = Proj(LatLon, boundary=false)
+  g = Triangle(merc(0, 0), merc(1, 0), merc(1, 1))
+  r, c = TB.apply(f, g)
+  @test r isa Triangle
+  f = Proj(LatLon, boundary=true)
+  r, c = TB.apply(f, g)
+  @test r isa TransformedGeometry
 end
 
 @testitem "Morphological" setup = [Setup] begin
   f = Morphological(c -> c)
   @test !isaffine(f)
   @test !TB.isrevertible(f)
   @test !TB.isinvertible(f)
-  @test TB.parameters(f) == (; fun=f.fun)
+  @test TB.parameters(f) == (fun=f.fun, boundary=false)
+  @test sprint(show, f) == "Morphological(fun: $(f.fun), boundary: false)"
+  @test sprint(show, MIME"text/plain"(), f) == """
+  Morphological transform
+  ├─ fun: $(f.fun)
+  └─ boundary: false"""
 
   # ----
   # VEC
@@ -1590,6 +1616,18 @@ end
   d = SimpleMesh(p, c)
   r, c = TB.apply(f, d)
   @test r ≈ SimpleMesh(f.(vertices(d)), topology(d))
+
+  # ----------------
+  # BOUNDARY OPTION
+  # ----------------
+
+  f = Morphological(c -> Cartesian(c.x, c.y, zero(c.x)), boundary=false)
+  g = Triangle(cart(0, 0), cart(1, 0), cart(1, 1))
+  r, c = TB.apply(f, g)
+  @test r isa Triangle
+  f = Morphological(c -> Cartesian(c.x, c.y, zero(c.x)), boundary=true)
+  r, c = TB.apply(f, g)
+  @test r isa TransformedGeometry
 end
 
 @testitem "LengthUnit" setup = [Setup] begin