diff --git a/ext/GeometryOpsDimensionalDataExt/GeometryOpsDimensionalDataExt.jl b/ext/GeometryOpsDimensionalDataExt/GeometryOpsDimensionalDataExt.jl
index 03ad84878..81f2c9f4b 100644
--- a/ext/GeometryOpsDimensionalDataExt/GeometryOpsDimensionalDataExt.jl
+++ b/ext/GeometryOpsDimensionalDataExt/GeometryOpsDimensionalDataExt.jl
@@ -4,8 +4,10 @@ import DimensionalData as DD
 import GeometryOps as GO
 import GeoInterface as GI
 
-function GO.polygonize(A::DD.AbstractDimArray; dims=(DD.X(), DD.Y()), crs=GI.crs(A), kw...)
+function GO.polygonize(A::DD.AbstractDimArray{T, N}; dims=(DD.XDim, DD.YDim), crs=GI.crs(A), kw...) where {T, N}
+    # Extract the lookups of the dimensions we care about
     lookups = DD.lookup(A, dims)
+    # Convert them to interval-bound vectors
     bounds_vecs = if DD.isintervals(lookups)
         map(DD.intervalbounds, lookups)
     else
@@ -14,7 +16,24 @@ function GO.polygonize(A::DD.AbstractDimArray; dims=(DD.X(), DD.Y()), crs=GI.crs
             DD.intervalbounds(DD.set(l, DD.Intervals()))
         end
     end
-    GO.polygonize(bounds_vecs..., DD.AbstractDimArray; crs, kw...)
+
+    # This tree switches on the array type.
+    # If ndims < 2, then we can't polygonize anyway, so we error.
+    # If ndims > 2, then we return a DimArray across the remaining dimensions
+    # not provided in `dims`.  Each slice in `dims` (usually X and Y) 
+    if N < 2
+        error("Array had $N dimensions, but it must have at least two!")
+    elseif N == 2
+        # If ndims == 2, then we polygonize directly, and return the output as requested.
+        Ap = PermutedDimsArray(A, dims)
+        GO.polygonize(bounds_vecs..., Ap; crs, kw...)
+    else
+        map(eachslice(A; dims = DD.otherdims(A, dims))) do a
+            ap = PermutedDimsArray(a, dims)
+            GO.polygonize(bounds_vecs..., a; crs, kw...)
+        end
+    end
+            
 end