diff --git a/ext/RastersArchGDALExt/RastersArchGDALExt.jl b/ext/RastersArchGDALExt/RastersArchGDALExt.jl
index 32e568e5e..63dce90e2 100644
--- a/ext/RastersArchGDALExt/RastersArchGDALExt.jl
+++ b/ext/RastersArchGDALExt/RastersArchGDALExt.jl
@@ -20,7 +20,7 @@ using Rasters: GDALsource, AbstractProjected, RasterStackOrArray, FileArray,
     RES_KEYWORD, SIZE_KEYWORD, CRS_KEYWORD, FILENAME_KEYWORD, SUFFIX_KEYWORD, EXPERIMENTAL,
     GDAL_EMPTY_TRANSFORM, GDAL_TOPLEFT_X, GDAL_WE_RES, GDAL_ROT1, GDAL_TOPLEFT_Y, GDAL_ROT2, GDAL_NS_RES
 
-import Rasters: reproject, resample, warp
+import Rasters: reproject, resample, warp, cellsize
 
 const RA = Rasters
 const DD = DimensionalData
@@ -28,6 +28,7 @@ const DA = DiskArrays
 const GI = GeoInterface
 const LA = LookupArrays
 
+include("cellsize.jl")
 include("gdal_source.jl")
 include("reproject.jl")
 include("resample.jl")
diff --git a/ext/RastersArchGDALExt/cellsize.jl b/ext/RastersArchGDALExt/cellsize.jl
new file mode 100644
index 000000000..833008d3c
--- /dev/null
+++ b/ext/RastersArchGDALExt/cellsize.jl
@@ -0,0 +1,106 @@
+"""
+    cellsize(x)
+
+Gives the approximate size of each cell in square km. 
+This function works for any projection, using an algorithm for polygons on a sphere. It approximates the true size to about 0.1%, depending on latitude.
+
+# Arguments
+
+- `x`: A `Raster` or a `Tuple` of `X` and `Y` dimensions. 
+
+## Example
+
+```jldoctest
+using Rasters, Rasters.LookupArrays
+dimz = X(Projected(90.0:0.1:120; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(0.1), crs=EPSG(4326))),
+       Y(Projected(0.0:0.1:50; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(0.1), crs=EPSG(4326)))
+
+cs = cellsize(dimz)
+
+```
+$EXPERIMENTAL
+"""
+
+function _great_circle_bearing(lon1::AbstractFloat, lat1::AbstractFloat, lon2::AbstractFloat, lat2::AbstractFloat)
+    dLong = lon1 - lon2
+
+    s = cos(lat2)*sin(dLong)
+    c = cos(lat1)*sin(lat2) - sin(lat1)*cos(lat2)*cos(dLong)
+
+    return atan(s, c)
+end
+
+## Get the area of a LinearRing with coordinates in radians
+# Using Gidard's theorem
+function _area_from_rads(ring; R = 6371.0088)
+    n = GI.npoint(ring)
+    area = -(n-3)*pi
+
+    prevpoint = GI.getpoint(ring, n-1)
+    point = GI.getpoint(ring, 1)
+    
+    for i in 2:n
+        nextpoint = GI.getpoint(ring, i)
+
+        beta1 = _great_circle_bearing(GI.x(point), GI.y(point), GI.x(prevpoint), GI.y(prevpoint)) 
+        beta2 = _great_circle_bearing(GI.x(point), GI.y(point), GI.x(nextpoint), GI.y(nextpoint)) 
+        angle = acos(cos(-beta1)*cos(-beta2) + sin(-beta1)*sin(-beta2))
+        area += angle
+
+        prevpoint = point
+        point = nextpoint
+    end
+    
+    return area*R^2
+end
+
+_area_from_coords(transform, geom) = _area_from_coords(transform, GI.trait(geom), geom)
+_area_from_coords(geom) = _area_from_coords(GI.trait(geom), geom)
+function _area_from_coords(::GI.LinearRingTrait, ring) # no ArchGDAL, assumes degrees
+    points = map(GI.getpoint(ring)) do p 
+        (deg2rad(GI.x(p)), deg2rad(GI.y(p)))
+    end
+
+    return _area_from_rads(GI.LinearRing(points))
+end
+function _area_from_coords(transform::ArchGDAL.CoordTransform, ::GI.LinearRingTrait, ring)
+    points = map(GI.getpoint(ring)) do p 
+        t = ArchGDAL.transform!(ArchGDAL.createpoint(p...), transform)
+        (deg2rad(GI.x(t)), deg2rad(GI.y(t)))
+    end
+    return _area_from_rads(GI.LinearRing(points))
+end
+
+function cellsize(dims::Tuple{X, Y})
+    xbnds, ybnds = DimensionalData.intervalbounds(dims)
+    if convert(CoordSys, crs(dims)) == CoordSys("Earth Projection 1, 104") # check if need to reproject
+        areas = [_area_from_coords(
+            GI.LinearRing([
+                (xb[1], yb[1]), 
+                (xb[2], yb[1]), 
+                (xb[2], yb[2]), 
+                (xb[1], yb[2]),
+                (xb[1], yb[1])
+            ]))
+            for xb in xbnds, yb in ybnds]
+    else 
+        areas = ArchGDAL.crs2transform(crs(dims), EPSG(4326)) do transform
+            [_area_from_coords(
+                transform,         
+                GI.LinearRing([
+                    (xb[1], yb[1]), 
+                    (xb[2], yb[1]), 
+                    (xb[2], yb[2]), 
+                    (xb[1], yb[2]),
+                    (xb[1], yb[1])
+                ]))
+                for xb in xbnds, yb in ybnds]
+        end
+    end
+
+    return Raster(areas, dims)
+end
+
+function cellsize(x::Raster)
+    cellsize(dims(x, (X, Y)))
+end
diff --git a/src/Rasters.jl b/src/Rasters.jl
index 0af7c82ab..a982badb0 100644
--- a/src/Rasters.jl
+++ b/src/Rasters.jl
@@ -60,7 +60,7 @@ export missingval, boolmask, missingmask, replace_missing, replace_missing!,
        aggregate, aggregate!, disaggregate, disaggregate!, mask, mask!, 
        resample, warp, zonal, crop, extend, trim, slice, combine, points,
        classify, classify!, mosaic, mosaic!, extract, rasterize, rasterize!,
-       coverage, coverage!, setcrs, setmappedcrs, smapseries
+       coverage, coverage!, setcrs, setmappedcrs, smapseries, cellsize
 export crs, mappedcrs, mappedindex, mappedbounds, projectedindex, projectedbounds
 export reproject, convertlookup
 
diff --git a/src/extensions.jl b/src/extensions.jl
index 27efd4c91..0fdf95afe 100644
--- a/src/extensions.jl
+++ b/src/extensions.jl
@@ -1,29 +1,22 @@
 # extensions
-
-# stubs that need ArchGDAL
-function resample(args...; kw...)
+function throw_extension_error(f::Function, package::String, extension::Symbol, args)
     @static if isdefined(Base, :get_extension) # julia > 1.9
-        if isnothing(Base.get_extension(Rasters, :RastersArchGDALExt))
-            throw(BackendException("ArchGDAL"))
-        else
-            throw(MethodError(resample, args))
-        end
+    if isnothing(Base.get_extension(Rasters, extension))
+        throw(BackendException(package))
     else
-        throw(BackendException("ArchGDAL"))
+        throw(MethodError(f, args))
     end
-end
-function warp(args...; kw...)
-    @static if isdefined(Base, :get_extension) # julia > 1.9
-        if isnothing(Base.get_extension(Rasters, :RastersArchGDALExt))
-            throw(BackendException("ArchGDAL"))
-        else
-            throw(MethodError(warp, args))
-        end
     else
-        throw(BackendException("ArchGDAL"))
+        throw(BackendException(package))
     end
 end
 
+
+# stubs that need ArchGDAL
+resample(args...; kw...) = throw_extension_error(resample, "ArchGDAL", :RastersArchGDALExt, args)
+warp(args...; kw...) = throw_extension_error(warp, "ArchGDAL", :RastersArchGDALExt, args)
+cellsize(args...; kw...) = throw_extension_error(cellsize, "ArchGDAL", :RastersArchGDALExt, args)
+
 # Other shared stubs
 function layerkeys end
 function smapseries end
diff --git a/test/cellsize.jl b/test/cellsize.jl
new file mode 100644
index 000000000..f5bea825c
--- /dev/null
+++ b/test/cellsize.jl
@@ -0,0 +1,25 @@
+using Rasters, Rasters.LookupArrays, ArchGDAL
+using Test
+
+include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
+
+@testset "cellsize" begin
+    dimz = X(Projected(90.0:0.1:99.9; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(0.1), crs=EPSG(4326))),
+       Y(Projected(0.0:0.1:89.9; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(0.1), crs=EPSG(4326)))
+
+    dimz_25832 = X(Projected(0.0:100:10000.0; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(100), crs=EPSG(25832))),
+       Y(Projected(0.0:100:10000.0; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(100), crs=EPSG(25832)))
+
+    cs = cellsize(dimz)
+    cs2 = cellsize(dimz_25832)
+
+    # Check the output is a raster 
+    @test cs isa Raster
+    @test cs2 isa Raster
+    # Test that the total area matches the expected area (1/72th of the Earth surface)
+    @test sum(cs) ≈ 510.1e6/72 rtol = 0.01
+    # Test all areas are about 0.01 km2
+    @test maximum(cs2) ≈ 0.01 rtol = 0.01
+    @test minimum(cs2) ≈ 0.01 rtol = 0.01
+
+end