optimize cellsize for lon-lat projections

Project.toml

@@ -16,6 +16,7 @@ FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 Flatten = "4c728ea3-d9ee-5c9a-9642-b6f7d7dc04fa"
 GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Missings = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
 Mmap = "a63ad114-7e13-5084-954f-fe012c677804"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
@@ -82,8 +83,8 @@ ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
 CFTime = "179af706-886a-5703-950a-314cd64e0468"
 CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
-GeoDataFrames = "62cb38b5-d8d2-4862-a48e-6a340996859f"
 GRIBDatasets = "82be9cdb-ee19-4151-bdb3-b400788d9abc"
+GeoDataFrames = "62cb38b5-d8d2-4862-a48e-6a340996859f"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"

ext/RastersArchGDALExt/RastersArchGDALExt.jl

@@ -21,16 +21,19 @@ using Rasters: GDALsource, AbstractProjected, RasterStackOrArray, FileArray, NoK
     GDAL_EMPTY_TRANSFORM, GDAL_TOPLEFT_X, GDAL_WE_RES, GDAL_ROT1, GDAL_TOPLEFT_Y, GDAL_ROT2, GDAL_NS_RES,
     _no_crs_error
 
-import Rasters: reproject, resample, warp, cellsize, nokw
+import Rasters: reproject, resample, warp, _spherical_cellarea, nokw
+
+import LinearAlgebra: dot, cross
 
 const RA = Rasters
+const AG = ArchGDAL
 const DD = DimensionalData
 const DA = DiskArrays
 const GI = GeoInterface
 const LA = Lookups
 
-include("cellsize.jl")
 include("gdal_source.jl")
+include("cellarea.jl")
 include("reproject.jl")
 include("resample.jl")
 include("warp.jl")

ext/RastersArchGDALExt/cellarea.jl

@@ -0,0 +1,100 @@
+## Get the area of a LinearRing with coordinates in radians
+struct SphericalPoint{T <: Real}
+	data::NTuple{3, T}
+end
+SphericalPoint(x, y, z) = SphericalPoint((x, y, z))
+
+# define the 4 basic mathematical operators elementwise on the data tuple
+Base.:+(p::SphericalPoint, q::SphericalPoint) = SphericalPoint(p.data .+ q.data)
+Base.:-(p::SphericalPoint, q::SphericalPoint) = SphericalPoint(p.data .- q.data)
+Base.:*(p::SphericalPoint, q::SphericalPoint) = SphericalPoint(p.data .* q.data)
+Base.:/(p::SphericalPoint, q::SphericalPoint) = SphericalPoint(p.data ./ q.data)
+# Define sum on a SphericalPoint to sum across its data
+Base.sum(p::SphericalPoint) = sum(p.data)
+
+# define dot and cross products
+dot(p::SphericalPoint, q::SphericalPoint) = sum(p * q)
+function cross(a::SphericalPoint, b::SphericalPoint)
+	a1, a2, a3 = a.data
+    b1, b2, b3 = b.data
+	SphericalPoint((a2*b3-a3*b2, a3*b1-a1*b3, a1*b2-a2*b1))
+end
+
+function _spherical_quadrilateral_area(ring)
+    (p1, p2, p3, p4) = _lonlat_to_sphericalpoint.(GI.getpoint(ring))
+    area = 0.0
+    area += _spherical_triangle_area(p1, p2, p3)
+    area += _spherical_triangle_area(p3, p4, p1)
+end
+
+# Using Eriksson's formula for the area of spherical triangles: https://www.jstor.org/stable/2691141
+function _spherical_triangle_area(a, b, c)
+    #t = abs(dot(a, cross(b, c)))
+    #t /= 1 + dot(b,c) + dot(c, a) + dot(a, b)
+    t = abs(dot(a, (cross(b - a, c - a))) / dot(b + a, c + a))
+    2*atan(t)
+end
+
+_lonlat_to_sphericalpoint(args) = _lonlat_to_sphericalpoint(args...)
+function _lonlat_to_sphericalpoint(lon, lat)
+    x = cosd(lat) * cosd(lon)
+    y = cosd(lat) * sind(lon)
+    z = sind(lat)
+    return SphericalPoint(x,y,z)
+end
+
+_area_from_coords(transform, geom) = _area_from_coords(transform, GI.trait(geom), geom)
+function _area_from_coords(transform::AG.CoordTransform, trait::GI.AbstractCurveTrait, ring)
+    t = AG.transform!(GI.convert(AG.geointerface_geomtype(trait), ring), transform)
+    return _spherical_quadrilateral_area(t)
+end
+
+# For lat-lon projections. Get the area of each latitudinal band, then multiply by the width
+function _area_from_lonlat(lon::XDim, lat::YDim; radius)
+    two_pi_R2 = 2 * pi * radius * radius
+    band_area = broadcast(DD.intervalbounds(lat)) do yb
+        two_pi_R2 * (sind(yb[2]) - sind(yb[1]))
+    end
+
+    broadcast(DD.intervalbounds(lon), band_area') do xb, ba
+        abs((xb[2] - xb[1]) / 360 * ba)
+    end
+end
+
+function _spherical_cellarea(dims::Tuple{<:XDim, <:YDim}; radius = 6371008.8)
+    # check the dimensions 
+    isnothing(crs(dims)) && _no_crs_error()
+
+    areas = if _isdegrees(crs(dims)) # check if need to reproject
+        _area_from_lonlat(dims...; radius)
+    elseif !isnothing(mappedcrs(dims)) && _isdegrees(mappedcrs(dims))
+        _area_from_lonlat(reproject(dims; crs = mappedcrs(dims))...; radius)
+    else
+        xbnds, ybnds = DD.intervalbounds(dims)
+        R2 = radius * radius
+        AG.crs2transform(crs(dims), EPSG(4326); order = :trad) do transform
+            [_area_from_coords(
+                transform,         
+                GI.LinearRing([
+                    (xb[1], yb[1]), 
+                    (xb[2], yb[1]), 
+                    (xb[2], yb[2]), 
+                    (xb[1], yb[2]),
+                ])
+                ) * R2
+                for xb in xbnds, yb in ybnds]
+        end
+    end
+end
+
+# TODO: put these in ArchGDAL
+_isgeographic(crs) = _isgeographic(AG.importCRS(crs))
+_isgeographic(crs::AG.ISpatialRef) = AG.GDAL.osrisgeographic(crs) |> Bool
+
+_isdegrees(crs) = _isdegrees(AG.importCRS(crs))
+function _isdegrees(crs::AG.ISpatialRef)
+    _isgeographic(crs) || return false
+    pointer = Ref{Cstring}()
+    result = AG.GDAL.osrgetangularunits(crs, pointer)
+    return unsafe_string(pointer[]) == "degree"
+end

ext/RastersArchGDALExt/gdal_source.jl

@@ -1,5 +1,3 @@
-const AG = ArchGDAL
-
 const GDAL_LOCUS = Start()
 
 const GDAL_DIM_ORDER = (X(), Y(), Band())

src/Rasters.jl

@@ -61,7 +61,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, cellsize
+       coverage, coverage!, setcrs, setmappedcrs, smapseries, cellsize, cellarea
 export crs, mappedcrs, mappedindex, mappedbounds, projectedindex, projectedbounds
 export reproject, convertlookup
 

src/extensions.jl

@@ -157,43 +157,74 @@
 warp(args...; kw...) = throw_extension_error(warp, "ArchGDAL", :RastersArchGDALExt, args)
 
 """
-    cellsize(x)
+    cellarea(x; radius = 6371008.8, area_in_crs = false)
 
-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.
+Gives the approximate area of each gridcell of `x`.
+By assuming the earth is a sphere, it approximates the true size to about 0.1%, depending on latitude.
 
-Run `using ArchGDAL` to make this method available.
-
-# Arguments
+Run `using ArchGDAL` to make this method fully available.
 
-- `x`: A `Raster` or a `Tuple` of `X` and `Y` dimensions.
+# Keywords
+- `radius`: the radius of the sphere of the coordinates. 
+    Defaults to the arithmetic mean radius of the earth in meters.
+- `area_in_crs`: if `true`, returns the area in the units of the crs of `x`, without any reprojection. 
+    This is equal to the distance between the bounds of each gridcell. `ArchGDAL` does not to be loaded if this is set to `true`.
+    `false` by default.
 
 ## Example
 
 ```julia
 using Rasters, ArchGDAL, Rasters.Lookups
-dimz = X(Projected(90.0:10.0:120; sampling=Intervals(Start()), order=ForwardOrdered(), span=Regular(10.0), crs=EPSG(4326))),
-       Y(Projected(0.0:10.0:50; sampling=Intervals(Start()), order=ForwardOrdered(), span=Regular(10.0), crs=EPSG(4326)))
-
-cs = cellsize(dimz)
+xdim = X(Projected(90.0:10.0:120; sampling=Intervals(Start()), crs=EPSG(4326)))
+ydim = Y(Projected(0.0:10.0:50; sampling=Intervals(Start()), crs=EPSG(4326)))
+myraster = rand(xdim, ydim)
+cs = cellarea(myraster)
 
 # output
-4×6 Raster{Float64,2} with dimensions:
-  X Projected{Float64} 90.0:10.0:120.0 ForwardOrdered Regular Intervals{Start} crs: EPSG,
-  Y Projected{Float64} 0.0:10.0:50.0 ForwardOrdered Regular Intervals{Start} crs: EPSG
-extent: Extent(X = (90.0, 130.0), Y = (0.0, 60.0))
-missingval: missing
-crs: EPSG:4326
-parent:
-        0.0       10.0       20.0        30.0        40.0            50.0
-  90.0  1.2332e6   1.1952e6   1.12048e6   1.01158e6   8.72085e5  706488.0
- 100.0  1.2332e6   1.1952e6   1.12048e6   1.01158e6   8.72085e5  706488.0
- 110.0  1.2332e6   1.1952e6   1.12048e6   1.01158e6   8.72085e5  706488.0
- 120.0  1.2332e6   1.1952e6   1.12048e6   1.01158e6   8.72085e5  706488.0
+╭───────────────────────╮
+│ 4×6 Raster{Float64,2} │
+├───────────────────────┴─────────────────────────────────────────────────── dims ┐
+  ↓ X Projected{Float64} 90.0:10.0:120.0 ForwardOrdered Regular Intervals{Start},
+  → Y Projected{Float64} 0.0:10.0:50.0 ForwardOrdered Regular Intervals{Start}
+├───────────────────────────────────────────────────────────────────────── raster ┤
+  extent: Extent(X = (90.0, 130.0), Y = (0.0, 60.0))
+
+  crs: EPSG:4326
+└─────────────────────────────────────────────────────────────────────────────────┘
+   ↓ →  0.0        10.0        20.0        30.0            40.0      50.0
+  90.0  1.23017e6   1.19279e6   1.11917e6   1.01154e6  873182.0  708290.0
+ 100.0  1.23017e6   1.19279e6   1.11917e6   1.01154e6  873182.0  708290.0
+ 110.0  1.23017e6   1.19279e6   1.11917e6   1.01154e6  873182.0  708290.0
+ 120.0  1.23017e6   1.19279e6   1.11917e6   1.01154e6  873182.0  708290.0
 ```
 $EXPERIMENTAL
 """
-cellsize(args...; kw...) = throw_extension_error(cellsize, "ArchGDAL", :RastersArchGDALExt, args)
+cellarea(x; kw...) = cellarea(dims(x, (XDim, YDim)); kw...)
+function cellarea(dims::Tuple{<:XDim, <:YDim}; radius = 6371008.8, area_in_crs = false)
+    isintervals(dims) || throw(ArgumentError("Cannot calculate cell size for a `Raster` with `Points` sampling."))
+    areas = if area_in_crs
+        _planar_cellarea(dims)
+    else
+       _spherical_cellarea(dims; radius)
+    end
+    return Raster(areas; dims)
+
+end
+
+_spherical_cellarea(args...; kw...) = throw_extension_error(_spherical_cellarea, "ArchGDAL", :RastersArchGDALExt, args)
+
+function _planar_cellarea(dims::Tuple{<:XDim, <:YDim})
+    xbnds, ybnds = DD.intervalbounds(dims)
+    broadcast(xb -> xb[2] - xb[1], xbnds) .* broadcast(yb -> yb[2] - yb[1], ybnds)'
+end
+
+function cellsize(args...; kw...) 
+    @warn """
+cellsize is deprecated and will be removed in a future version, use cellarea instead. 
+Note that cellarea returns the area in square m, while cellsize still uses square km.
+"""
+    return cellarea(args...; kw..., radius = 6371.0088)
+end
 
 # Other shared stubs
 function layerkeys end