Skip to content

add DimensionalData extension #178

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 6 commits into
base: main
Choose a base branch
from
Open

add DimensionalData extension #178

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
56f51ae
add DimensionalData extension
rafaqz Jul 14, 2024
8ec11d4
bugfix imports
rafaqz Jul 14, 2024
e799c08
stop testing on 1.9
rafaqz Jul 14, 2024
d040a4c
set minimum julia version to 1.9
rafaqz Jul 14, 2024
d823118
Fix capitalization
asinghvi17 Jul 27, 2024
34c6668
Merge branch 'main' into dd
asinghvi17 Sep 17, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion .github/workflows/CI.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ jobs:
fail-fast: false
matrix:
version:
- '1.9'
- '1.10'
- '1'
# - 'nightly'
os:
Expand Down
5 changes: 4 additions & 1 deletion Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,11 +16,13 @@ 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"
Expand All @@ -29,6 +31,7 @@ GeometryOpsProjExt = "Proj"
CoordinateTransformations = "0.5, 0.6"
DelaunayTriangulation = "1.0.4"
ExactPredicates = "2.2.8"
DimensionalData = "0.27"
FlexiJoins = "0.1.30"
GeoInterface = "1.2"
GeometryBasics = "0.4.7"
Expand All @@ -38,7 +41,7 @@ Proj = "1"
SortTileRecursiveTree = "0.1"
Statistics = "1"
Tables = "1"
julia = "1.9"
julia = "1.10"

[extras]
ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
Expand Down
21 changes: 21 additions & 0 deletions ext/GeometryOpsDimensionalDataExt/GeometryOpsDimensionalDataExt.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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...)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What happens with the dims that aren't given (time etc)?

Copy link
Member Author

@rafaqz rafaqz Jul 15, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah I haven't done that yet, this should be WIP.

But I image we make a feature collection with a feature for each slice, and put the lookup values as properties?

Otherwise return a DimArray of geoms of the remaining dimensions ?

Probably a single feature collection or table is most useful?

Copy link
Member Author

@rafaqz rafaqz Jul 29, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually this is a case for R style vector data cubes.

But we can just broadcast over the other dimensions and get a DimArray of X/Y geometries with time/etc dimensions.

We should be able to write that directly as a FeatureCollection with GeoJSON/ArchGDAL/Shapefile without any more work as long as we name the array :geometry.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we then want to add an assert here that the dims can only have length 2?

lookups = DD.lookup(A, dims)
bounds_vecs = if DD.isintervals(lookups)
map(DD.intervalbounds, lookups)

Check warning on line 10 in ext/GeometryOpsDimensionalDataExt/GeometryOpsDimensionalDataExt.jl

View check run for this annotation

Codecov / codecov/patch

ext/GeometryOpsDimensionalDataExt/GeometryOpsDimensionalDataExt.jl#L7-L10 

Added lines #L7 - L10 were not covered by tests
else
@warn "`polygonsize` is not possible for `Points` sampling, as polygons cover space by definition. Treating as `Intervals`, but this may not be appropriate"
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Doesn't this just mean that the ring vertices are located at the points? That seems pretty exact...

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm not sure I understand you. It's just saying that points don't cover any space so polygonizing them doesn't make sense

Copy link
Member

@asinghvi17 asinghvi17 Jan 22, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I was thinking of it more in a sense that the polygonizing operation is the concave hull of the set of points that satisfy f...but that is also probably a different algorithm!

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, somehing else! polygonize is like a groupby/combine for lots of little square polygons

map(lookups) do l
DD.intervalbounds(DD.set(l, DD.Intervals()))

Check warning on line 14 in ext/GeometryOpsDimensionalDataExt/GeometryOpsDimensionalDataExt.jl

View check run for this annotation

Codecov / codecov/patch

ext/GeometryOpsDimensionalDataExt/GeometryOpsDimensionalDataExt.jl#L12-L14 

Added lines #L12 - L14 were not covered by tests
end
end
GO.polygonize(bounds_vecs..., DD.AbstractDimArray; crs, kw...)

Check warning on line 17 in ext/GeometryOpsDimensionalDataExt/GeometryOpsDimensionalDataExt.jl

View check run for this annotation

Codecov / codecov/patch

ext/GeometryOpsDimensionalDataExt/GeometryOpsDimensionalDataExt.jl#L17 

Added line #L17 was not covered by tests
end


end
79 changes: 44 additions & 35 deletions test/methods/polygonize.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,7 @@
using Test
import GeometryOps: polygonize
import DimensionalData as DD
import Rasters, DimensionalData
using GeometryOps, GeoInterface, Test
using ..TestHelpers

Expand Down Expand Up @@ -25,7 +29,19 @@ import GeoInterface as GI
data_mp_range_floats = polygonize(range_floats, range_floats, data2)
end

import OffsetArrays, DimensionalData, Rasters
import OffsetArrays

@testset "Polygonize with xs and ys, with offsetarrays" begin
data = rand(1:4, 100, 100) .== 1
unitrange = 1:100
steprange = 1:1:100
steprangelen = range(1, 100; length = 100)
data_mp = polygonize(data)
for range in (unitrange, steprange, steprangelen)
data_mp_range = polygonize(range, range, data)
@test GO.equals(data_mp, data_mp_range)
end
end

# Missing holes throw a warning, so testing there are
# no warnings in a range of randomisation is one way to test
Expand All @@ -34,33 +50,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
Expand All @@ -77,31 +93,24 @@ 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

@testset "Polygonize with xs and ys, with offsetarrays" begin
data = rand(1:4, 100, 100) .== 1
unitrange = 1:100
steprange = 1:1:100
steprangelen = range(1, 100; length = 100)
data_mp = polygonize(data)
for range in (unitrange, steprange, steprangelen)
data_mp_range = polygonize(range, range, data)
@test GO.equals(data_mp, data_mp_range)
end
end
Loading