add DimensionalData extension

.github/workflows/CI.yml

@@ -18,7 +18,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.9'
+          - '1.10'
           - '1'
           # - 'nightly'
         os:

Project.toml

@@ -15,18 +15,21 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
 [weakdeps]
+DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0"
 FlexiJoins = "e37f2e79-19fa-4eb7-8510-b63b51fe0a37"
 LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb"
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 
 [extensions]
+GeometryOpsDimensionalDataExt = "DimensionalData"
 GeometryOpsFlexiJoinsExt = "FlexiJoins"
 GeometryOpsLibGEOSExt = "LibGEOS"
 GeometryOpsProjExt = "Proj"
 
 [compat]
 CoordinateTransformations = "0.5, 0.6"
 ExactPredicates = "2.2.8"
+DimensionalData = "0.27"
 FlexiJoins = "0.1.30"
 GeoInterface = "1.2"
 GeometryBasics = "0.4.7"
@@ -36,7 +39,7 @@ Proj = "1"
 SortTileRecursiveTree = "0.1"
 Statistics = "1"
 Tables = "1"
-julia = "1.9"
+julia = "1.10"
 
 [extras]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"

ext/GeometryOpsDimensionalDataExt/GeometryOpsDimensionalDataExt.jl

@@ -0,0 +1,21 @@
+module GeometryOpsDimensionalDataExt
+
+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...)
+    lookups = DD.lookup(A, dims)
+    bounds_vecs = if DD.isintervals(lookups)
+        map(DD.intervalbounds, lookups)
+    else
+        @warn "`polygonsize` is not possible for `Points` sampling, as polygons cover space by definition. Treating as `Intervals`, but this may not be appropriate"
+        map(lookups) do l
+            DD.intervalbounds(DD.set(l, DD.Intervals()))
+        end
+    end
+    GO.polygonize(bounds_vecs..., DD.AbstractDimArray; crs, kw...)
+end
+
+
+end

test/methods/polygonize.jl

@@ -1,5 +1,9 @@
-using GeometryOps, GeoInterface, Test
-import OffsetArrays, DimensionalData, Rasters
+using Test
+import GeometryOps as GO
+import GeometryOps: polygonize
+import GeoInterface as GI
+import DimensionalData as DD
+import OffsetArrays, Rasters
 
 # Missing holes throw a warning, so testing there are
 # no warnings in a range of randomisation is one way to test 
@@ -8,33 +12,33 @@ for i in (100, 300), j in (100, 300)
     @testset "bool arrays without a function return MultiPolygon" begin
         A = rand(Bool, i, j)
         @test_nowarn multipolygon = polygonize(A);
-        @test multipolygon isa GeoInterface.MultiPolygon
-        @test GeoInterface.ngeom(multipolygon) > 0
+        @test multipolygon isa GI.MultiPolygon
+        @test GI.ngeom(multipolygon) > 0
     end
 
     A = rand(i, j)
     @testset "bool functions return MultiPolygon" begin
         multipolygon = @test_nowarn polygonize(>(0.5), A);
-        @test multipolygon isa GeoInterface.MultiPolygon
-        @test GeoInterface.ngeom(multipolygon) > 0
+        @test multipolygon isa GI.MultiPolygon
+        @test GI.ngeom(multipolygon) > 0
     end
 
     @testset "other functions return FeatureCollection" begin
         fc = @test_nowarn polygonize(x -> trunc(3x), A);
-        @test fc isa GeoInterface.FeatureCollection
-        @test GeoInterface.nfeature(fc) == 3
-        @test map(GeoInterface.getfeature(fc)) do f
-            GeoInterface.properties(f).value
+        @test fc isa GI.FeatureCollection
+        @test GI.nfeature(fc) == 3
+        @test map(GI.getfeature(fc)) do f
+            GI.properties(f).value
         end == [0.0, 1.0, 2.0]
     end
 
     @testset "values are polygonized without a function" begin
         A = rand(1:3, i, j)
         fc = @test_nowarn polygonize(A)
-        fc isa GeoInterface.FeatureCollection
-        @test GeoInterface.nfeature(fc) == 3
-        @test map(GeoInterface.getfeature(fc)) do f
-            GeoInterface.properties(f).value
+        fc isa GI.FeatureCollection
+        @test GI.nfeature(fc) == 3
+        @test map(GI.getfeature(fc)) do f
+            GI.properties(f).value
         end == [1, 2, 3]
     end
 end
@@ -51,19 +55,23 @@ end
         end
         @test GO.equals(data_mp, evil_in_data_space_mp)
     end
+end
+
+@testset "Polygonize with DimensionalData compatible arrays" begin
+    data = rand(1:4, 100, 50) .== 1
+    dd = DD.DimArray(data, (DD.X(51:150), DD.Y(151:200)))
     @testset "DimensionalData" begin
-        data = rand(1:4, 100, 100) .== 1
-        evil = DimensionalData.DimArray(data, (DimensionalData.X(1:100), DimensionalData.Y(1:100)))
-        data_mp = polygonize(data)
-        evil_mp = @test_nowarn polygonize(evil)
-        @test GO.equals(data_mp, evil_mp)
+        data_mp = polygonize(51:150, 151:200, data);
+        dd_mp = polygonize(dd);
+        @test GO.equals(data_mp, dd_mp)
     end
+
     @testset "Rasters" begin
-        data = rand(1:4, 100, 100) .== 1
-        evil = Rasters.Raster(data; dims = (DimensionalData.X(1:100), DimensionalData.Y(1:100)), crs = Rasters.GeoFormatTypes.EPSG(4326))
-        data_mp = polygonize(data)
-        evil_mp = @test_nowarn polygonize(evil)
-        @test GO.equals(data_mp, evil_mp)
-        @test GI.crs(evil_mp) == GI.crs(evil)
+        data = rand(1:4, 100, 50) .== 1
+        rast = Rasters.Raster(data; dims=(DD.X(51:150), DD.Y(151:200)), crs=Rasters.GeoFormatTypes.EPSG(4326))
+        data_mp = polygonize(51:150, 151:200, data)
+        rast_mp = @test_nowarn polygonize(rast)
+        @test GO.equals(data_mp, rast_mp)
+        @test GI.crs(rast_mp) == GI.crs(evil)
     end
-end
+end