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

Add tests for unitful cellarea #799

Open
wants to merge 5 commits into
base: main
Choose a base branch
from
+30 −1
Open

Add tests for unitful cellarea #799

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
a07be84
Add tests for unitful `cellarea`
asinghvi17 Oct 11, 2024
bcd7fb1
Fix cellarea tests
asinghvi17 Oct 11, 2024
f5fce8b
more test fixes
asinghvi17 Oct 12, 2024
3360f34
Merge branch 'main' into test_unitful_cellarea
rafaqz Nov 12, 2024
5a14cbe
Merge branch 'main' into test_unitful_cellarea
rafaqz Dec 13, 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
3 changes: 2 additions & 1 deletion Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,7 +102,8 @@ RasterDataSources = "3cb90ccd-e1b6-4867-9617-4276c8b2ca36"
SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[targets]
test = ["Aqua", "ArchGDAL", "CFTime", "CoordinateTransformations", "DataFrames", "GeoDataFrames", "GeometryBasics", "GRIBDatasets", "NCDatasets", "Plots", "Proj", "RasterDataSources", "SafeTestsets", "Shapefile", "Statistics", "Test", "ZarrDatasets"]
test = ["Aqua", "ArchGDAL", "CFTime", "CoordinateTransformations", "DataFrames", "GeoDataFrames", "GeometryBasics", "GRIBDatasets", "NCDatasets", "Plots", "Proj", "RasterDataSources", "SafeTestsets", "Shapefile", "Statistics", "Test", "Unitful", "ZarrDatasets"]
26 changes: 26 additions & 0 deletions test/cellarea.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
using Rasters, DimensionalData, Rasters.Lookups, Proj
using Test
using DimensionalData: @dim, YDim
import Unitful
include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))

@testset "cellarea" begin
Expand Down Expand Up @@ -70,4 +71,29 @@ include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
# test missing crs throws an error
nocrsdimz = setcrs(dimz, nothing)
@test_throws ArgumentError cellarea(nocrsdimz)

@testset "Unitful cellarea" begin
# Case 1: cellarea with unitful manifold
# This is the simplest case
unitful_manifold = Spherical(; radius = Spherical().radius * Unitful.u"m")
unitful_cellarea = cellarea(unitful_manifold, dimz_3857)
@test unitful_cellarea isa Raster{<:Unitful.Quantity}
@test Unitful.ustrip.(parent(unitful_cellarea)) == cellarea(Spherical(; radius = unitful_manifold.radius |> Unitful.ustrip), dimz_3857)
# Case 2: unitful dimensions
ux_3857 = rebuild(x_3857; val = rebuild(val(x_3857); data = val(x_3857) .* Unitful.u"m", span = Regular(val(x_3857).span.step * Unitful.u"m")))
uy_3857 = rebuild(y_3857; val = rebuild(val(y_3857); data = val(y_3857) .* Unitful.u"m", span = Regular(val(y_3857).span.step * Unitful.u"m")))
unitful_dimz_3857 = (ux_3857, uy_3857)
Comment on lines +83 to +85
Copy link
Collaborator Author

Choose a reason for hiding this comment

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

This is also super complicated, should we have DD overloads to make constructing lookups with unitful easier?

Copy link
Owner

@rafaqz rafaqz Nov 12, 2024 •
edited
Loading

Choose a reason for hiding this comment

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

I think set will just work this out for us, no need to rebuild manually?

@test cellarea(Planar(), unitful_dimz_3857) isa Raster{<:Unitful.Quantity}
@test Unitful.ustrip.(parent(cellarea(Planar(), unitful_dimz_3857))) == parent(cellarea(Planar(), dimz_3857))
# Unitful lookups shouldn't work without a unitful manifold
@test_throws Unitful.DimensionError cellarea(unitful_dimz_3857)
# The tests below fail because Unitful apparently can't convert to Float64...
# (see https://github.com/PainterQubits/Unitful.jl/issues/742)
# But we also have to re-convert back to the original unit type, which cellarea currently
# doesn't do.

Comment on lines +90 to +94
Copy link
Owner

@rafaqz rafaqz Oct 13, 2024 •
edited
Loading

Choose a reason for hiding this comment

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

Suggested change
# The tests below fail because Unitful apparently can't convert to Float64...
# (see https://github.com/PainterQubits/Unitful.jl/issues/742)
# But we also have to re-convert back to the original unit type, which cellarea currently
# doesn't do.
# This needs fixes for unit handling, probably using `ustrip(u"m", x)` in a Unitful.jl extension. That will also require checking the units of the crs.

Unitful is like a free test suite for physical models, not converting is by design!

If we are going to make units work here we need to do ustrip(u"m", x) to ensure we are actually getting metres coming in. If people pass in kilometers and we convert back and forth with Proj.jl treating them as meters in the middle, we are really breaking the contract.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

That makes sense. Maybe we have a maybe_ustrip(unit, x) defined in Rasters, its definition for Unitful unit in RastersUnitfulExt, and its definition for Proj units in RastersProjUnitfulExt?

(This is getting to be a mouthful..)

Copy link
Owner

Choose a reason for hiding this comment

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

Yeah those names! But makes sense. Somehow we need some mapping from Proj units to Unitful.

Copy link
Contributor

Choose a reason for hiding this comment

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

Why does conversion from Proj to Unitful have to live in Rasters.jl? Shouldn't this be in a ProjUnitfulExt?

Copy link
Collaborator Author

@asinghvi17 asinghvi17 Oct 28, 2024 •
edited
Loading

Choose a reason for hiding this comment

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

Hmm, we could introduce a nice interface to batch transform an array, but do it out of place in Proj.jl. Then the Unitful stuff could also live there.

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 was initially going to put that here because this calls the C-API directly, but if it doesn't have to then there's no reason to stop this from calling Proj.

# unitful_cellarea = cellarea(unitful_manifold, unitful_dimz_3857)
# @test unitful_cellarea isa Raster{<:Unitful.Quantity}
# @test Unitful.ustrip.(unitful_cellarea) == cellarea(unitful_manifold, dimz_3857)
end
end
Loading