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.

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

Fix _extent when called on an extent #852

Merged
merged 10 commits into from
Jan 17, 2025
1 change: 1 addition & 0 deletions src/methods/burning/extents.jl
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,7 @@ function _extent(::GI.AbstractFeatureCollectionTrait, features; kw...)::XYExtent
end
return _float64_xy_extent(ext)
end
_extent(ext::Extent; kw...) = _float64_xy_extent(ext)
Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I thought about this one, but do we really want to force the convert to float64? This will mean that rasterize with an extent in integers or float32 will return float64 dimensions

Copy link
Owner

Choose a reason for hiding this comment

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

We convert everything else already. I think there were 2 reasons - accuracy of ext + res/size is much better with Float64, and it reduces compilation for force the extent to always be the same object.

But we can change everything back to allow Float32 in another PR if you think its necessary

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I think we want _extent(ext::Extent; kw...) = Extents.Extent(X=ext.X, Y=ext.Y). If ext doesn't have an X or Y field it will give a reasonable error message.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Okay if it gets converted elsewhere then nevermind

Copy link
Owner

Choose a reason for hiding this comment

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

We actually want to force ::XYExtent on the method to be consistent

Copy link
Owner

Choose a reason for hiding this comment

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

We can have the Float32 converstion elsewhere but lets conform with the current setup for now

Copy link
Contributor

Choose a reason for hiding this comment

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

Yeah, but without that my actual example failed when I used a Float32 extent.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Just checked and I don't think this gets converted anywhere else, though. With the above definition of _extent I can rasterize to Float32 dimensions no problem. I haven't run into this yet, but don't see why we shouldn't allow it.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Okay I didn't realize XYExtent has the Float64 type hard-coded in. Guess that is needed for type stability?

Copy link
Owner

Choose a reason for hiding this comment

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

I think type stability was one reason, this recursive calls were not compiling away and forcing the type fixes that


function _float64_xy_extent(ext::Extents.Extent)
xbounds = map(Float64, ext.X)
Expand Down
7 changes: 7 additions & 0 deletions test/methods.jl
Original file line number Diff line number Diff line change
Expand Up @@ -762,3 +762,10 @@ test = rebuild(ga; name = :test)
@test_throws "strictly positive" Rasters.sample(StableRNG(123), test, 3, skipmissing = true, replace = false)
@test_throws "Cannot draw" Rasters.sample(StableRNG(123), test, 5, replace = false)
end

@testset "extent" begin
ga = Raster(A, (X(1.0:1:2.0), Y(1.0:1:2.0)); missingval=missing)
Copy link
Contributor

Choose a reason for hiding this comment

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

We should also test this for an Extent, that is based on float32

Copy link
Owner

Choose a reason for hiding this comment

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

@tiemvanderdeure we need this to merge

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I thought you said you wanted to fix that in a separate PR? Right now this would fail, I guess?

Copy link
Owner

Choose a reason for hiding this comment

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

Will it? I don't totally understand, won't it just convert to Float64?

We just need to check the input works

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

What I meant is the Rasters._extent(ext) === ext will fail if ext is float32-based. But I've added some tests anyway

ext = extent(ga)
@test ext === Extent(X=(1.0,2.0), Y=(1.0,2.0))
@test Rasters._extent(ext) === ext
end
4 changes: 1 addition & 3 deletions test/resample.jl
Original file line number Diff line number Diff line change
Expand Up @@ -64,9 +64,7 @@ include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
resampled = resample(cea; method)
@test crs(cea) == crs(resampled)
@test cea == resampled
# There is some floating point error here after Rasters -> GDAL -> Rasterss...
# Should we correct it by detecting almost identical extent and using the original?
# @test_broken extent(cea) == extent(resampled)
@test extent(cea) == extent(resampled)
end

@testset "only `res` kw changes the array size predictably" begin
Expand Down
Loading