Fix _extent when called on an extent (#852)

* add a dispatch to _extent so it returns itself

* Update src/methods/burning/extents.jl

* Update src/methods/burning/extents.jl

* test extent

* uncomment a test that isn't broken anymore

* use Rasters._extent

Co-authored-by: Rafael Schouten <[email protected]>

* forward keywords

* make better use of _get_geometries

* add a test for float32

* fix test

---------

Co-authored-by: Rafael Schouten <[email protected]>

Author: tiemvanderdeure

src/methods/burning/extents.jl

@@ -1,18 +1,20 @@
 const XYExtent = Extents.Extent{(:X,:Y),Tuple{Tuple{Float64,Float64},Tuple{Float64,Float64}}}
 
 # Get the bounds of a geometry
-_extent(geom; kw...)::XYExtent = _extent(GI.trait(geom), geom; kw...)
-function _extent(::Nothing, data; geometrycolumn, kw...)::XYExtent
+function _extent(data; geometrycolumn=nothing, kw...)::XYExtent 
     geoms = _get_geometries(data, geometrycolumn)
-
+    _extent(GI.trait(geoms), geoms)
+end
+function _extent(::Nothing, geoms)::XYExtent
+    # because geoms was returned from _get_geometries, it must be an iterable of valid geometries
     g1 = first(geoms)
     if GI.trait(g1) isa GI.PointTrait 
         xs = extrema(p -> GI.x(p), geoms)
         ys = extrema(p -> GI.y(p), geoms)
         return _float64_xy_extent(Extents.Extent(X=xs, Y=ys))
     else
-        ext = reduce(geoms; init=_extent(first(geoms))) do ext, geom
-            Extents.union(ext, _extent(geom))
+        ext = reduce(geoms; init=_extent(GI.trait(g1), g1)) do ext, geom
+            Extents.union(ext, _extent(GI.trait(geom), geom))
         end
         return _float64_xy_extent(ext)
     end
@@ -33,23 +35,8 @@ function _extent(::GI.AbstractGeometryTrait, geom; kw...)::XYExtent
         return _float64_xy_extent(geomextent)
     end
 end
-_extent(::GI.AbstractFeatureTrait, feature; kw...)::XYExtent = _extent(GI.geometry(feature); kw...)
-function _extent(::GI.GeometryCollectionTrait, collection; kw...)::XYExtent
-    geometries = GI.getgeom(collection)
-    init = _float64_xy_extent(_extent(first(geometries)))
-    ext = reduce(geometries; init) do acc, g
-        Extents.union(acc, _extent(g))
-    end
-    return _float64_xy_extent(ext)
-end
-function _extent(::GI.AbstractFeatureCollectionTrait, features; kw...)::XYExtent
-    features = GI.getfeature(features)
-    init = _float64_xy_extent(_extent(first(features)))
-    ext = reduce(features; init) do acc, f
-        Extents.union(acc, _extent(f))
-    end
-    return _float64_xy_extent(ext)
-end
+
+_extent(ext::Extent; kw...) = _float64_xy_extent(ext)
 
 function _float64_xy_extent(ext::Extents.Extent)
     xbounds = map(Float64, ext.X)

test/methods.jl

@@ -762,3 +762,15 @@ test = rebuild(ga; name = :test)
     @test_throws "strictly positive" Rasters.sample(StableRNG(123), test, 3, skipmissing = true, replace = false)
     @test_throws "Cannot draw" Rasters.sample(StableRNG(123), test, 5, replace = false)
 end
+
+@testset "extent" begin
+    ga = Raster(A, (X(1.0:1:2.0), Y(1.0:1:2.0)); missingval=missing) 
+    ext = extent(ga)
+    @test ext === Extent(X=(1.0,2.0), Y=(1.0,2.0))
+    @test Rasters._extent(ext) === ext
+
+    ga2 = Raster(A, (X(Float32.(1:2)), Y(Float32.(1:2))))
+    ext2 = extent(ga2)
+    @test ext2 === Extent(X=(1.0f0,2.0f0), Y=(1.0f0,2.0f0))
+    @test Rasters._extent(ext2) === ext # currently this converts to float64!
+end

test/resample.jl

@@ -64,9 +64,7 @@ include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
         resampled = resample(cea; method)
         @test crs(cea) == crs(resampled)
         @test cea == resampled
-        # There is some floating point error here after Rasters -> GDAL -> Rasterss...
-        # Should we correct it by detecting almost identical extent and using the original?
-        # @test_broken extent(cea) == extent(resampled)
+        @test extent(cea) == extent(resampled)
     end
 
     @testset "only `res` kw changes the array size predictably" begin