diff --git a/src/fallbacks.jl b/src/fallbacks.jl
index 62dc836..c8d0119 100644
--- a/src/fallbacks.jl
+++ b/src/fallbacks.jl
@@ -126,10 +126,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)
@@ -142,6 +142,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/src/interface.jl b/src/interface.jl
index a0d98b5..b261dbd 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)
diff --git a/test/test_wrappers.jl b/test/test_wrappers.jl
index 5c343e7..4ff20bc 100644
--- a/test/test_wrappers.jl
+++ b/test/test_wrappers.jl
@@ -366,7 +366,6 @@ vecfc = GI.FeatureCollection([(geometry=(1,2), a=1, b=2)])
 @test GI.getfeature(vecfc, 1) == (geometry=(1,2), a=1, b=2)
 
 
-
 struct MaPointRappa
     x::Float64
     y::Float64
@@ -393,6 +392,21 @@ end
     )
 end
 
+@testset "Extent skips NaN (#166)" begin
+    @testset "NaNs mixed with real values" 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
+    @testset "All NaN extent should be NaN" begin
+        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
+end
+
 # TODO
 
 # # Triangle