diff --git a/Project.toml b/Project.toml
index 8ef210f98..d0a502536 100644
--- a/Project.toml
+++ b/Project.toml
@@ -108,7 +108,10 @@ Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+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", "StableRNGs", "StatsBase", "Test", "ZarrDatasets"]
+test = ["Aqua", "ArchGDAL", "CFTime", "CoordinateTransformations", "DataFrames", "GeoDataFrames", "GeometryBasics", "GRIBDatasets", "NCDatasets", "Plots", "Proj", "RasterDataSources", "SafeTestsets", "Shapefile", "StableRNGs", "Statistics", "StatsBase", "Test", "Unitful", "ZarrDatasets"]
+
diff --git a/test/cellarea.jl b/test/cellarea.jl
index 8f4c66ac3..7b3d7ffd5 100644
--- a/test/cellarea.jl
+++ b/test/cellarea.jl
@@ -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
@@ -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)
+         @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.
+
+         # 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