Add ValidCoords transform (#1152)

* Add 'InDomain' transform

* Add docstring

* Add to documentation && Add to exports

* Add tests

* Update docstring

* Update example

* Add more tests

* More adjustments

* Rename InDomain to ValidCoords

* Rename helper

* Apply suggestions

* Apply suggestions

* Apply suggestions from code review

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

* Use the correct parent type

* Apply suggestions from code review

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

* Apply suggestions

* Add more tests

* Adjust code order

---------

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

docs/src/transforms.md

@@ -173,6 +173,26 @@ using Unitful: m, cm
 Point(1m, 2m, 3m) |> LengthUnit(cm)
 ```
 
+## ValidCoords
+
+```@docs
+ValidCoords
+```
+
+```@example transforms
+# load coordinate reference system
+using CoordRefSystems: LatLon
+
+# regular grid with LatLon coordinates
+grid = RegularGrid(Point(LatLon(-90, -180)), Point(LatLon(90, 180)), dims=(10, 10))
+
+# retain elements in the projection domain
+subgrid = grid |> ValidCoords(Mercator)
+
+# plot the projected grid
+viz(subgrid |> Proj(Mercator), showsegments=true)
+```
+
 ## Shadow
 
 ```@docs

src/Meshes.jl

@@ -542,6 +542,7 @@ export
   Proj,
   Morphological,
   LengthUnit,
+  ValidCoords,
   Shadow,
   Slice,
   Repair,

src/transforms.jl

@@ -96,6 +96,7 @@ include("transforms/stdcoords.jl")
 include("transforms/proj.jl")
 include("transforms/morphological.jl")
 include("transforms/lengthunit.jl")
+include("transforms/validcoords.jl")
 include("transforms/shadow.jl")
 include("transforms/slice.jl")
 include("transforms/repair.jl")

src/transforms/validcoords.jl

@@ -0,0 +1,54 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+"""
+    ValidCoords(CRS)
+    ValidCoords(code)
+
+Retain the geometries within the domain of the
+projection of type `CRS` or with EPSG/ESRI `code`.
+"""
+struct ValidCoords{CRS} <: GeometricTransform end
+
+ValidCoords(CRS) = ValidCoords{CRS}()
+
+ValidCoords(code::Type{<:EPSG}) = ValidCoords(CoordRefSystems.get(code))
+
+ValidCoords(code::Type{<:ESRI}) = ValidCoords(CoordRefSystems.get(code))
+
+parameters(::ValidCoords{CRS}) where {CRS} = (; CRS)
+
+preprocess(t::ValidCoords, d::Domain) = findall(g -> _isvalid(t, g), d)
+
+function preprocess(t::ValidCoords, d::Mesh)
+  points = vertices(d)
+  topo = topology(d)
+  findall(elements(topo)) do elem
+    is = indices(elem)
+    all(_isvalid(t, points[i]) for i in is)
+  end
+end
+
+apply(t::ValidCoords, d::Domain) = view(d, preprocess(t, d)), nothing
+
+# -----------
+# IO METHODS
+# -----------
+
+Base.show(io::IO, ::ValidCoords{CRS}) where {CRS} = print(io, "ValidCoords(CRS: $CRS)")
+
+function Base.show(io::IO, ::MIME"text/plain", t::ValidCoords{CRS}) where {CRS}
+  summary(io, t)
+  println(io)
+  print(io, "└─ CRS: $CRS")
+end
+
+# -----------------
+# HELPER FUNCTIONS
+# -----------------
+
+_isvalid(::ValidCoords{CRS}, p::Point) where {CRS} = indomain(CRS, coords(p))
+_isvalid(t::ValidCoords, g::Polytope) = all(p -> _isvalid(t, p), eachvertex(g))
+_isvalid(t::ValidCoords, g::MultiPolytope) = all(p -> _isvalid(t, p), eachvertex(g))
+_isvalid(t::ValidCoords, g::Geometry) = all(p -> _isvalid(t, p), pointify(g))

test/transforms.jl

@@ -1782,6 +1782,104 @@ end
   @test r ≈ SimpleMesh(f.(vertices(d)), topology(d))
 end
 
+@testitem "ValidCoords" setup = [Setup] begin
+  @test !isaffine(ValidCoords(Mercator))
+  @test !TB.isrevertible(ValidCoords(Mercator))
+  @test !TB.isinvertible(ValidCoords(Mercator))
+  @test TB.parameters(ValidCoords(Mercator)) == (; CRS=Mercator)
+  @test TB.parameters(ValidCoords(EPSG{3395})) == (; CRS=Mercator{WGS84Latest})
+  @test TB.parameters(ValidCoords(ESRI{54017})) == (; CRS=Behrmann{WGS84Latest})
+  f = ValidCoords(Mercator)
+  @test sprint(show, f) == "ValidCoords(CRS: CoordRefSystems.Mercator)"
+  @test sprint(show, MIME"text/plain"(), f) == """
+  ValidCoords transform
+  └─ CRS: CoordRefSystems.Mercator"""
+
+  # ---------
+  # POINTSET
+  # ---------
+
+  f = ValidCoords(Mercator)
+  d = PointSet([latlon(-90, 0), latlon(-45, 0), latlon(0, 0), latlon(45, 0), latlon(90, 0)])
+  r, c = TB.apply(f, d)
+  @test r == PointSet([latlon(-45, 0), latlon(0, 0), latlon(45, 0)])
+
+  # ------------
+  # GEOMETRYSET
+  # ------------
+
+  f = ValidCoords(Mercator)
+  t1 = Triangle(latlon(-90, 0), latlon(-45, 30), latlon(-45, -30))
+  t2 = Triangle(latlon(-45, 0), latlon(0, 30), latlon(0, -30))
+  t3 = Triangle(latlon(0, -30), latlon(0, 30), latlon(45, 0))
+  t4 = Triangle(latlon(45, -30), latlon(45, 30), latlon(90, 0))
+  d = GeometrySet([t1, t2, t3, t4])
+  r, c = TB.apply(f, d)
+  @test r == GeometrySet([t2, t3])
+
+  f = ValidCoords(Mercator)
+  t1 = Triangle(latlon(-90, 0), latlon(-45, 30), latlon(-45, -30))
+  t2 = Triangle(latlon(-45, 0), latlon(0, 30), latlon(0, -30))
+  t3 = Triangle(latlon(0, -30), latlon(0, 30), latlon(45, 0))
+  t4 = Triangle(latlon(45, -30), latlon(45, 30), latlon(90, 0))
+  m1 = Multi([t1, t4])
+  m2 = Multi([t2, t3])
+  d = GeometrySet([m1, m2])
+  r, c = TB.apply(f, d)
+  @test r == GeometrySet([m2])
+
+  f = ValidCoords(Mercator)
+  b1 = Box(latlon(-90, -30), latlon(-30, 30))
+  b2 = Box(latlon(-30, -30), latlon(30, 30))
+  b3 = Box(latlon(30, -30), latlon(90, 30))
+  d = GeometrySet([b1, b2, b3])
+  r, c = TB.apply(f, d)
+  @test r == GeometrySet([b2])
+
+  # --------------
+  # CARTESIANGRID
+  # --------------
+
+  f = ValidCoords(Mercator)
+  d = RegularGrid(latlon(-90, -180), latlon(90, 180), dims=(3, 3))
+  r, c = TB.apply(f, d)
+  linds = LinearIndices(size(d))
+  @test r == view(d, [linds[2, 1], linds[2, 2], linds[2, 3]])
+
+  # ----------------
+  # RECTILINEARGRID
+  # ----------------
+
+  f = ValidCoords(Mercator)
+  g = RegularGrid(latlon(-90, -180), latlon(90, 180), dims=(3, 3))
+  d = convert(RectilinearGrid, g)
+  r, c = TB.apply(f, d)
+  linds = LinearIndices(size(d))
+  @test r == view(d, [linds[2, 1], linds[2, 2], linds[2, 3]])
+
+  # ---------------
+  # STRUCTUREDGRID
+  # ---------------
+
+  f = ValidCoords(Mercator)
+  g = RegularGrid(latlon(-90, -180), latlon(90, 180), dims=(3, 3))
+  d = convert(StructuredGrid, g)
+  r, c = TB.apply(f, d)
+  linds = LinearIndices(size(d))
+  @test r == view(d, [linds[2, 1], linds[2, 2], linds[2, 3]])
+
+  # -----------
+  # SIMPLEMESH
+  # -----------
+
+  f = ValidCoords(Mercator)
+  g = RegularGrid(latlon(-90, -180), latlon(90, 180), dims=(3, 3))
+  d = convert(SimpleMesh, g)
+  r, c = TB.apply(f, d)
+  linds = LinearIndices(size(d))
+  @test r == view(d, [linds[2, 1], linds[2, 2], linds[2, 3]])
+end
+
 @testitem "Shadow" setup = [Setup] begin
   @test !isaffine(Shadow(:xy))
   @test !TB.isrevertible(Shadow("xy"))