skip NaN in extent calculation

src/fallbacks.jl

@@ -123,10 +123,10 @@ function calc_extent(t::AbstractPointTrait, geom)
 end
 function calc_extent(t::AbstractGeometryTrait, geom)
  points = getpoint(t, geom)
- X = extrema(p -> x(p), points)
- Y = extrema(p -> y(p), points)
+ X = _nan_free_extrema(x, points)
+ Y = _nan_free_extrema(y, points)
  if is3d(geom)
- Z = extrema(p -> z(p), points)
+ Z = _nan_free_extrema(z, points)
  Extent(; X, Y, Z)
  else
  Extent(; X, Y)
@@ -139,6 +139,24 @@ function calc_extent(::AbstractFeatureTrait, feature)
 end
 calc_extent(t::AbstractFeatureCollectionTrait, fc) = reduce(Extents.union, filter(!isnothing, collect(extent(f) for f in getfeature(t, fc))))
 
+function _nan_free_extrema(f, points)
+ agg_min = agg_max = f(first(points))
+ found = false
+ for point in points
+ c = f(point)
+ isnan(c) && continue
+ if found
+ agg_min = min(agg_min, c)
+ agg_max = max(agg_max, c)
+ else
+ found = true
+ agg_min = c
+ agg_max = c
+ end
+ end
+ return agg_min, agg_max
+end
+
 # Package level `GeoInterface.convert` method
 # Packages must implement their own `traittype` method
 # that accepts a GeoInterface.jl trait and returns the

test/test_wrappers.jl

@@ -301,6 +301,18 @@ fc = GI.FeatureCollection(fc_unwrapped.parent; crs=EPSG(4326), extent=GI.extent(
 vecfc = GI.FeatureCollection([(geometry=(1,2), a=1, b=2)])
 @test GI.getfeature(vecfc, 1) == (geometry=(1,2), a=1, b=2)
 
+@testset "Extent skips NaN" begin
+ nan_multipoint = GI.MultiPoint([(1.0, NaN), (3.0, 4.0), (3.0, 2.0), (NaN, 4.0), (7.0, 8.0), (9.0, 10.0)])
+ nan_polygon = GI.Polygon(GI.LineString([(1.0, NaN), (3.0, 4.0), (3.0, 2.0), (NaN, 4.0), (7.0, 8.0), (9.0, 10.0)]))
+ @test GI.extent(nan_multipoint) == Extent(X = (1.0, 9.0), Y = (2.0, 10.0))
+ @test GI.extent(nan_polygon) == Extent(X = (1.0, 9.0), Y = (2.0, 10.0))
+
+ all_nan_multipoint = GI.MultiPoint([(NaN, NaN) for i in 1:10])
+ all_nan_ext = GI.extent(all_nan_multipoint)
+ @test all(isnan, all_nan_ext.X)
+ @test all(isnan, all_nan_ext.Y)
+end
+
 # TODO
 
 # # Triangle