From b6ad8df9460ac29ad28be7c96cad18fe5612d0c2 Mon Sep 17 00:00:00 2001 From: rafaqz Date: Thu, 26 Sep 2024 00:03:40 +0200 Subject: [PATCH 1/3] skip nan in extent calculation --- src/fallbacks.jl | 24 +++++++++++++++++++++--- test/test_wrappers.jl | 7 +++++++ 2 files changed, 28 insertions(+), 3 deletions(-) diff --git a/src/fallbacks.jl b/src/fallbacks.jl index 9c501ad1..7af6e91f 100644 --- a/src/fallbacks.jl +++ b/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 diff --git a/test/test_wrappers.jl b/test/test_wrappers.jl index d1f35edc..8a8131ca 100644 --- a/test/test_wrappers.jl +++ b/test/test_wrappers.jl @@ -301,6 +301,13 @@ 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)) +end + # TODO # # Triangle From 8722337dd6be21ce49e66c3dc3c9dd9d2aaa933f Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Thu, 26 Sep 2024 00:32:55 +0200 Subject: [PATCH 2/3] test all nan Co-authored-by: Anshul Singhvi --- test/test_wrappers.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/test/test_wrappers.jl b/test/test_wrappers.jl index 8a8131ca..3485afb4 100644 --- a/test/test_wrappers.jl +++ b/test/test_wrappers.jl @@ -306,6 +306,11 @@ vecfc = GI.FeatureCollection([(geometry=(1,2), a=1, b=2)]) 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 From 893ce3577b5b5a3b71b6c0fa7cc309a44bd69710 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 22 Jan 2025 14:26:53 -0500 Subject: [PATCH 3/3] Add documentation on extent calculation skipping NaNs --- src/interface.jl | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/interface.jl b/src/interface.jl index 9369435a..75ed5a75 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -453,10 +453,16 @@ Retrieve the extent (bounding box) for given geom or feature. In SF this is defined as `envelope`. `Extents.extent(obj)` will be called if `extent(trait(obj), obj)`, -is not defined so it may be preferable to define `Extents.extent` directly. +is not defined, so it is preferable to define `Extents.extent` directly. -When `fallback` is true, and the obj does not have an extent, +When `fallback` is true, and the `obj` does not have an extent, an extent is calculated from the coordinates of all geometries in `obj`. + +!!! note NaN values + The GeoInterface extent calculation **ignores** all NaN values in coordinates, + but if all coordinates are NaN, it returns a NaN extent. + + For example, the extent of `[(0, 0), (NaN, 1)]` is `X = (0, 0), Y = (0, 1)`. """ function extent(obj; fallback=true) t = trait(obj)