Skip to content
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

Jump to bottom

Type stable calc extent #98

Merged
merged 14 commits into from
Aug 7, 2023
Merged

Type stable calc extent #98

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
f039844
fix calc extent type stability
rafaqz Apr 29, 2023
955dd3e
optimise point extent
rafaqz Apr 29, 2023
e5847cd
add ncoord to test objects
rafaqz Apr 29, 2023
8f26310
qualify xyz
rafaqz Apr 29, 2023
6cd4790
test 3d extents
rafaqz Jul 29, 2023
0761687
add missing extents point wrapper method
rafaqz Jul 29, 2023
a185bb1
test all wrapper extents
rafaqz Jul 29, 2023
e5f8fe8
update tests
rafaqz Jul 31, 2023
2824bb1
f
rafaqz Jul 31, 2023
c20cc06
show
rafaqz Jul 31, 2023
8e9ebc4
f
rafaqz Jul 31, 2023
f9f19f7
f
rafaqz Jul 31, 2023
549e0ab
f
rafaqz Jul 31, 2023
b813824
remove accidental liines
rafaqz Aug 4, 2023
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
24 changes: 20 additions & 4 deletions src/fallbacks.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -111,11 +111,27 @@ coordinates(t::AbstractFeatureCollectionTrait, fc) = [coordinates(f) for f in ge

extent(::AbstractTrait, x) = Extents.extent(x)
function calc_extent(t::AbstractPointTrait, geom)
coords = collect(getcoord(t, geom))
names = coordnames(geom)
return Extent(NamedTuple{names}(zip(coords, coords)))
x = GeoInterface.x(t, geom)
y = GeoInterface.y(t, geom)
if is3d(geom)
z = GeoInterface.z(t, geom)
return Extent(; X=(x, x), Y=(y, y), Z=(z, z))
else
return Extent(; X=(x, x), Y=(y, y))
end
end
calc_extent(t::AbstractGeometryTrait, geom) = reduce(Extents.union, (extent(f) for f in getgeom(t, geom)))
function calc_extent(t::AbstractGeometryTrait, geom)
points = getpoint(t, geom)
X = extrema(p -> x(p), points)
Y = extrema(p -> y(p), points)
if is3d(geom)
Z = extrema(p -> z(p), points)
Extent(; X, Y, Z)
else
Extent(; X, Y)
end
end
calc_extent(t::GeometryCollectionTrait, geom) = reduce(Extents.union, (extent(f) for f in getgeom(t, geom)))
Copy link
Member

Choose a reason for hiding this comment

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

pedantry: note that we call extent (and thus accept pre-calculated extents) for multigeoms instead of calling calc_extent (recalculating everything). Old behaviour and quite pragmatic, but it might be slightly unexpected.

function calc_extent(::AbstractFeatureTrait, feature)
geom = geometry(feature)
isnothing(geom) ? nothing : extent(geom)
Expand Down
1 change: 1 addition & 0 deletions src/wrappers.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -308,6 +308,7 @@ ncoord(trait::PointTrait, geom::Point) = ncoord(trait, parent(geom))
getcoord(trait::PointTrait, geom::Point, i::Integer) = getcoord(trait, parent(geom), i)
convert(::Type{Point}, ::PointTrait, geom) = Point(geom)
convert(::Type{Point}, ::PointTrait, geom::Point) = geom
extent(trait::PointTrait, geom::Point) = extent(trait, parent(geom))

x(trait::PointTrait, geom::Point) = x(trait, parent(geom))
y(trait::PointTrait, geom::Point) = y(trait, parent(geom))
Expand Down
16 changes: 15 additions & 1 deletion test/test_primitives.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,42 +35,49 @@ GeoInterface.isgeometry(::MyCurve) = true
GeoInterface.geomtrait(::MyCurve) = LineStringTrait()
GeoInterface.ngeom(::LineStringTrait, geom::MyCurve) = 2
GeoInterface.getgeom(::LineStringTrait, geom::MyCurve, i) = MyPoint()
GeoInterface.ncoord(t::LineStringTrait, geom::MyCurve) = 2

GeoInterface.ncoord(t::PolygonTrait, geom::MyPolygon) = 2
GeoInterface.isgeometry(::MyPolygon) = true
GeoInterface.geomtrait(::MyPolygon) = PolygonTrait()
GeoInterface.ngeom(::PolygonTrait, geom::MyPolygon) = 2
GeoInterface.getgeom(::PolygonTrait, geom::MyPolygon, i) = MyCurve()
GeoInterface.ncoord(t::PolygonTrait, geom::MyPolygon) = 2

GeoInterface.isgeometry(::MyTriangle) = true
GeoInterface.geomtrait(::MyTriangle) = TriangleTrait()
GeoInterface.ngeom(::TriangleTrait, geom::MyTriangle) = 3
GeoInterface.getgeom(::TriangleTrait, geom::MyTriangle, i) = MyCurve()
GeoInterface.ncoord(t::TriangleTrait, geom::MyTriangle) = 2

GeoInterface.isgeometry(::MyMultiPoint) = true
GeoInterface.geomtrait(::MyMultiPoint) = MultiPointTrait()
GeoInterface.ngeom(::MultiPointTrait, geom::MyMultiPoint) = 2
GeoInterface.getgeom(::MultiPointTrait, geom::MyMultiPoint, i) = MyPoint()
GeoInterface.ncoord(t::MultiPointTrait, geom::MyMultiPoint) = 2

GeoInterface.isgeometry(::MyMultiCurve) = true
GeoInterface.geomtrait(::MyMultiCurve) = MultiCurveTrait()
GeoInterface.ngeom(::MultiCurveTrait, geom::MyMultiCurve) = 2
GeoInterface.getgeom(::MultiCurveTrait, geom::MyMultiCurve, i) = MyCurve()
GeoInterface.ncoord(t::MultiCurveTrait, geom::MyMultiCurve) = 2

GeoInterface.isgeometry(::MyMultiPolygon) = true
GeoInterface.geomtrait(::MyMultiPolygon) = MultiPolygonTrait()
GeoInterface.ngeom(::MultiPolygonTrait, geom::MyMultiPolygon) = 2
GeoInterface.getgeom(::MultiPolygonTrait, geom::MyMultiPolygon, i) = MyPolygon()
GeoInterface.ncoord(t::MultiPolygonTrait, geom::MyMultiPolygon) = 2

GeoInterface.isgeometry(::MyTIN) = true
GeoInterface.geomtrait(::MyTIN) = PolyhedralSurfaceTrait()
GeoInterface.ngeom(::PolyhedralSurfaceTrait, geom::MyTIN) = 2
GeoInterface.getgeom(::PolyhedralSurfaceTrait, geom::MyTIN, i) = MyTriangle()
GeoInterface.ncoord(t::PolyhedralSurfaceTrait, geom::MyTIN) = 2

GeoInterface.isgeometry(::MyCollection) = true
GeoInterface.geomtrait(::MyCollection) = GeometryCollectionTrait()
GeoInterface.ngeom(::GeometryCollectionTrait, geom::MyCollection) = 2
GeoInterface.getgeom(::GeometryCollectionTrait, geom::MyCollection, i) = MyCurve()
GeoInterface.ncoord(t::GeometryCollectionTrait, geom::MyCollection) = 2

GeoInterface.isfeature(::Type{<:MyFeature}) = true
GeoInterface.trait(feature::MyFeature) = FeatureTrait()
Expand Down Expand Up @@ -117,6 +124,7 @@ GeoInterface.getfeature(::FeatureCollectionTrait, fc::MyFeatureCollection, i::In
@test testgeometry(geom)

@test GeoInterface.npoint(geom) == 2 # defaults to ngeom
@test GeoInterface.ncoord(geom) == 2
@test GeoInterface.coordinates(geom) == [[1, 2], [1, 2]]
points = GeoInterface.getpoint(geom)
point = GeoInterface.getpoint(geom, 1)
Expand All @@ -139,6 +147,7 @@ GeoInterface.getfeature(::FeatureCollectionTrait, fc::MyFeatureCollection, i::In

@test GeoInterface.nring(geom) == 2
@test GeoInterface.nhole(geom) == 1
@test GeoInterface.ncoord(geom) == 2
@test GeoInterface.coordinates(geom) == [[[1, 2], [1, 2]], [[1, 2], [1, 2]]]
lines = GeoInterface.getring(geom)
line = GeoInterface.getring(geom, 1)
Expand All @@ -162,6 +171,7 @@ GeoInterface.getfeature(::FeatureCollectionTrait, fc::MyFeatureCollection, i::In
@test testgeometry(geom)

@test GeoInterface.npoint(geom) == 2
@test GeoInterface.ncoord(geom) == 2
points = GeoInterface.getpoint(geom)
point = GeoInterface.getpoint(geom, 1)
@test GeoInterface.coordinates(geom) == [[1, 2], [1, 2]]
Expand All @@ -174,6 +184,7 @@ GeoInterface.getfeature(::FeatureCollectionTrait, fc::MyFeatureCollection, i::In
geom = MyMultiCurve()
@test testgeometry(geom)

@test GeoInterface.ncoord(geom) == 2
@test GeoInterface.nlinestring(geom) == 2
lines = GeoInterface.getlinestring(geom)
line = GeoInterface.getlinestring(geom, 1)
Expand All @@ -186,6 +197,7 @@ GeoInterface.getfeature(::FeatureCollectionTrait, fc::MyFeatureCollection, i::In
@test testgeometry(geom)

@test GeoInterface.npolygon(geom) == 2
@test GeoInterface.ncoord(geom) == 2
polygons = GeoInterface.getpolygon(geom)
polygon = GeoInterface.getpolygon(geom, 1)
@test GeoInterface.coordinates(geom) == [[[[1, 2], [1, 2]], [[1, 2], [1, 2]]], [[[1, 2], [1, 2]], [[1, 2], [1, 2]]]]
Expand All @@ -209,6 +221,7 @@ GeoInterface.getfeature(::FeatureCollectionTrait, fc::MyFeatureCollection, i::In
@test testgeometry(geom)

@test GeoInterface.ngeom(geom) == 2
@test GeoInterface.ncoord(geom) == 2
geoms = GeoInterface.getgeom(geom)
thing = GeoInterface.getgeom(geom, 1)
@test GeoInterface.coordinates(geom) == [[[1, 2], [1, 2]], [[1, 2], [1, 2]]]
Expand Down Expand Up @@ -318,6 +331,7 @@ end
@test GeoInterface.y(geom) == 2
@test_throws ArgumentError GeoInterface.z(geom)
@test_throws ArgumentError GeoInterface.m(geom)
@test GeoInterface.extent(geom) == Extents.Extent(X=(1, 1), Y=(2, 2))
geom = [1, 2, 3]
@test testgeometry(geom)
@test collect(GeoInterface.getcoord(geom)) == geom
Expand Down
30 changes: 30 additions & 0 deletions test/test_wrappers.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@ point = GI.Point(1, 2)
GI.getcoord(point, 1)
@test !GI.ismeasured(point)
@test !GI.is3d(point)
@test GI.extent(point) == Extent(X=(1, 1), Y=(2, 2))
@test point == GI.Point(point)
@test (GI.x(point), GI.y(point)) == (1, 2)
@test_throws ArgumentError GI.z(point)
Expand All @@ -25,6 +26,7 @@ pointz = GI.Point(1, 2, 3)
@test (GI.x(pointz), GI.y(pointz), GI.z(pointz)) == (1, 2, 3)
@test GI.testgeometry(pointz)
@test GI.convert(GI, pointz) === pointz
@test GI.extent(pointz) == Extents.Extent(X=(1, 1), Y=(2, 2), Z=(3, 3))

# 3D measured point
pointzm = GI.Point(; X=1, Y=2, Z=3, M=4)
Expand Down Expand Up @@ -101,6 +103,8 @@ line = GI.Line([(1, 2), (3, 4)])
@test GI.getgeom(line, 1) === (1, 2)
@test GI.getgeom(line) == [(1, 2), (3, 4)]
@test GI.testgeometry(line)
@test !GI.is3d(line)
@test GI.extent(line) == Extent(X=(1, 3), Y=(2, 4))
@test_throws ArgumentError GI.Line(point)
@test_throws ArgumentError GI.Line([(1, 2)])
@test_throws ArgumentError GI.Line([line, line])
Expand All @@ -114,6 +118,8 @@ linestring = GI.LineString([(1, 2), (3, 4)])
@test GI.getgeom(linestring, 1) === (1, 2)
@test GI.getgeom(linestring) == [(1, 2), (3, 4)]
@test GI.testgeometry(linestring)
@test !GI.is3d(linestring)
@test @inferred(GI.extent(linestring)) == Extent(X=(1, 3), Y=(2, 4))
@test_throws ArgumentError GI.LineString([(1, 2)])
linestring_crs = GI.LineString(linestring; crs=EPSG(4326))
@test parent(linestring_crs) === parent(linestring)
Expand All @@ -125,6 +131,8 @@ linearring = GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)])
@test GI.getgeom(linearring, 1) === (1, 2)
@test GI.getgeom(linearring) == [(1, 2), (3, 4), (5, 6), (1, 2)]
@test GI.testgeometry(linearring)
@test !GI.is3d(linearring)
@test @inferred(GI.extent(linearring)) == Extent(X=(1, 5), Y=(2, 6))
@test_throws ArgumentError GI.LinearRing([(1, 2)])
linearring_crs = GI.LinearRing(linearring; crs=EPSG(4326))
@test parent(linearring_crs) === parent(linearring)
Expand All @@ -137,16 +145,25 @@ polygon = GI.Polygon([linearring, linearring])
@test collect(GI.getgeom(polygon)) == [linearring, linearring]
@test collect(GI.getpoint(polygon)) == vcat(collect(GI.getpoint(linearring)), collect(GI.getpoint(linearring)))
@test GI.testgeometry(polygon)
@test !GI.is3d(polygon)
@test @inferred(GI.extent(polygon)) == Extent(X=(1, 5), Y=(2, 6))
@test GI.convert(GI, MyPolygon()) isa GI.Polygon
@test GI.convert(GI, polygon) === polygon
polygon_crs = GI.Polygon(polygon; crs=EPSG(4326))
@test parent(polygon_crs) === parent(polygon)
@test GI.crs(polygon_crs) === EPSG(4326)

linearring3d = GI.LinearRing([(1, 2, 3), (3, 4, 5), (5, 6, 7), (1, 2, 3)])
polygon3d = GI.Polygon([linearring3d, linearring3d])
@test GI.is3d(polygon3d)
@test GI.extent(polygon3d) == Extents.Extent(X=(1, 5), Y=(2, 6), Z=(3, 7))

# MultiPoint
multipoint = GI.MultiPoint([(1, 2), (3, 4), (3, 2), (1, 4), (7, 8), (9, 10)])
@test multipoint == GI.MultiPoint(multipoint)
@test GI.getgeom(multipoint, 1) === (1, 2)
@test !GI.is3d(multipoint)
@test @inferred(GI.extent(multipoint)) == Extent(X=(1, 9), Y=(2, 10))
@test_throws ArgumentError GI.MultiPoint([[(1, 2), (3, 4), (3, 2), (1, 4), (7, 8), (9, 10)]])
@test GI.testgeometry(multipoint)
multipoint_crs = GI.MultiPoint(multipoint; crs=EPSG(4326))
Expand All @@ -159,6 +176,8 @@ collection = GI.GeometryCollection(geoms)
@test collection == GI.GeometryCollection(collection)
@test GI.getgeom(collection) == geoms
@test GI.testgeometry(collection)
@test !GI.is3d(collection)
@test GI.extent(collection) == reduce(Extents.union, map(GI.extent, geoms))
collection_crs = GI.GeometryCollection(collection; crs=EPSG(4326))
@test parent(collection_crs) == parent(collection)
@test GI.crs(collection_crs) === EPSG(4326)
Expand All @@ -168,16 +187,25 @@ multicurve = GI.MultiCurve([linestring, linearring])
@test collect(GI.getpoint(multicurve)) == vcat(collect(GI.getpoint(linestring)), collect(GI.getpoint(linearring)))
@test multicurve == GI.MultiCurve(multicurve)
@test GI.getgeom(multicurve, 1) === linestring
@test !GI.is3d(multicurve)
@test GI.extent(multicurve) == Extent(X=(1, 5), Y=(2, 6))
@test_throws ArgumentError GI.MultiCurve([pointz, polygon])
@test GI.testgeometry(multicurve)
multicurve_crs = GI.MultiCurve(multicurve; crs=EPSG(4326))
@test parent(multicurve_crs) == parent(multicurve)
@test GI.crs(multicurve_crs) === EPSG(4326)


# MultiPolygon
polygon = GI.Polygon([linearring, linearring])
multipolygon = GI.MultiPolygon([polygon])
@test multipolygon == GI.MultiPolygon(multipolygon)
@test GI.getgeom(multipolygon, 1) === polygon
@test !GI.is3d(multipolygon)
@show polygon
@show GI.getgeom(polygon, 1)
# MultiPolygon extent does not infer, maybe due to nesting
@test GI.extent(multipolygon) == Extent(X=(1, 5), Y=(2, 6))
@test collect(GI.getpoint(multipolygon)) == collect(GI.getpoint(polygon))
@test_throws ArgumentError GI.MultiPolygon([[[[(1, 2), (3, 4), (3, 2), (1, 4)]]]])
@test GI.testgeometry(multipolygon)
Expand All @@ -188,6 +216,8 @@ multipolygon_crs = GI.MultiPolygon(multipolygon; crs=EPSG(4326))
# PolyhedralSurface
polyhedralsurface = GI.PolyhedralSurface([polygon, polygon])
@test polyhedralsurface == GI.PolyhedralSurface(polyhedralsurface)
@test !GI.is3d(polyhedralsurface)
@test @inferred(GI.extent(polyhedralsurface)) == Extent(X=(1, 5), Y=(2, 6))
@test GI.getgeom(polyhedralsurface, 1) === polygon
@test collect(GI.getgeom(polyhedralsurface)) == [polygon, polygon]
@test GI.getgeom(polyhedralsurface, 1) == polygon
Expand Down
Loading