[WIP] switch cellarea and reproject to use Proj

Project.toml

@@ -32,6 +32,7 @@ CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
 GRIBDatasets = "82be9cdb-ee19-4151-bdb3-b400788d9abc"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
+Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 RasterDataSources = "3cb90ccd-e1b6-4867-9617-4276c8b2ca36"
 ZarrDatasets = "519a4cdf-1362-424a-9ea1-b1d782dbb24b"
 
@@ -41,6 +42,7 @@ RastersCoordinateTransformationsExt = "CoordinateTransformations"
 RastersGRIBDatasetsExt = "GRIBDatasets"
 RastersMakieExt = "Makie"
 RastersNCDatasetsExt = "NCDatasets"
+RastersProjExt = "Proj"
 RastersRasterDataSourcesExt = "RasterDataSources"
 RastersZarrDatasetsExt = "ZarrDatasets"
 
@@ -68,6 +70,7 @@ Missings = "0.4, 1"
 NCDatasets = "0.13, 0.14"
 OffsetArrays = "1"
 ProgressMeter = "1"
+Proj = "1.3"
 RasterDataSources = "0.5.7, 0.6"
 RecipesBase = "0.7, 0.8, 1.0"
 Reexport = "0.2, 1.0"
@@ -91,11 +94,12 @@ GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 RasterDataSources = "3cb90ccd-e1b6-4867-9617-4276c8b2ca36"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "ArchGDAL", "CFTime", "CoordinateTransformations", "DataFrames", "GeoDataFrames", "GeometryBasics", "GRIBDatasets", "NCDatasets", "Plots", "RasterDataSources", "SafeTestsets", "Shapefile", "Statistics", "Test", "ZarrDatasets"]
+test = ["Aqua", "ArchGDAL", "CFTime", "CoordinateTransformations", "DataFrames", "GeoDataFrames", "GeometryBasics", "GRIBDatasets", "NCDatasets", "Plots", "Proj", "RasterDataSources", "SafeTestsets", "Shapefile", "Statistics", "Test", "ZarrDatasets"]

ext/RastersArchGDALExt/RastersArchGDALExt.jl

@@ -33,8 +33,6 @@ const GI = GeoInterface
 const LA = Lookups
 
 include("gdal_source.jl")
-include("cellarea.jl")
-include("reproject.jl")
 include("resample.jl")
 include("warp.jl")
 

ext/RastersMakieExt/plotrecipes.jl

@@ -305,9 +305,6 @@ else # surfacelike is not a thing
     Makie.convert_arguments(t::Makie.ImageLike, A::AbstractRaster{<: Any, 2}) = Makie.convert_arguments(t, _prepare_dimarray(A))
 end
 
-function Makie.convert_arguments(t::Makie.DiscreteSurface, A::AbstractRaster{<:Any,2})
-    return Makie.convert_arguments(t, _prepare_dimarray(A))
-end
 function Makie.convert_arguments(t::Makie.VolumeLike, A::AbstractRaster{<:Any,3}) 
     return Makie.convert_arguments(t, _prepare_dimarray(A))
 end

ext/RastersProjExt/RastersProjExt.jl

@@ -0,0 +1,33 @@
+module RastersProjExt
+
+using Rasters, Proj
+
+import DiskArrays,
+    Extents,
+    Missings
+
+using DimensionalData,
+    GeoFormatTypes,
+    GeoInterface
+
+
+using Rasters.Lookups
+using Rasters.Dimensions
+using Rasters: GDALsource, AbstractProjected, RasterStackOrArray, FileArray, NoKW,
+    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,
+    _no_crs_error
+
+import Rasters: reproject, _spherical_cellarea, nokw
+
+import LinearAlgebra: dot, cross
+
+const RA = Rasters
+const DD = DimensionalData
+const DA = DiskArrays
+const GI = GeoInterface
+const LA = Lookups
+
+include("cellarea.jl")
+include("reproject.jl")
+end

ext/RastersProjExt/cellarea.jl

@@ -0,0 +1,185 @@
+## 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)
+    # don't assume the ring is a GI ring
+    (p1, p2, p3, p4) = _lonlat_to_sphericalpoint.(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, geom)
+function _area_from_coords(transform::Proj.Transformation, ring_points)
+    t = (transform(GI.x(p), GI.y(p)) for p in ring_points)
+    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, use_area_of_use = true)
+    # 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
+        area_of_use = if use_area_of_use
+            # Use a temporary transformation object to get the area of use,
+            # since the transformation object has to be recreated with the area of use
+            _get_area_of_use(Proj.Transformation(crs(dims), EPSG(4326); always_xy = true), Extents.Extent(X = extrema(dims[1]), Y = extrema(dims[2])))
+        else
+            C_NULL
+        end
+
+        transform = Proj.Transformation(crs(dims), EPSG(4326); always_xy = true, area = area_of_use)
+
+        result = [
+                _area_from_coords(
+                    transform,         
+                    (
+                        (xb[1], yb[1]), 
+                        (xb[2], yb[1]), 
+                        (xb[2], yb[2]), 
+                        (xb[1], yb[2]),
+                    )
+                ) * R2
+                for xb in xbnds, yb in ybnds
+            ]
+
+        if use_area_of_use
+            Proj.proj_area_destroy(area_of_use)
+        end
+
+        return result
+    end
+end
+
+function _get_area_of_use(transform::Proj.Transformation, extent::Extents.Extent; densify_pts = 21)
+    # Transform the extent using `proj_trans_bounds`
+    (xmin, xmax), (ymin, ymax) = Proj.bounds(transform, extent.X, extent.Y; densify_pts)
+
+    # Create an area of use object 
+    # This MUST be destroyed by the caller
+    area = Proj.proj_area_create()
+
+    # Set the bounding box of the area of use
+    Proj.proj_area_set_bbox(
+        area, 
+        xmin, # west_lon_degree
+        ymin, # south_lat_degree
+        xmax, # east_lon_degree
+        ymax  # north_lat_degree
+    )
+
+    return area
+end
+
+# TODO: put these in Proj (specifically the dispatches on GFT types)
+_isgeographic(crs) = _isgeographic(Proj.CRS(crs))
+_isgeographic(crs::Proj.CRS) = Proj.is_geographic(crs)
+
+_isdegrees(crs) = _isdegrees(Proj.CRS(crs))
+function _isdegrees(crs::Proj.CRS)
+    _isgeographic(crs) || return false
+    # This is a tiny bit inefficient, but it takes 500ns on my machine, 
+    # so I think we can disregard the inefficiency...
+    return axis_is_degrees(crs, 0) && axis_is_degrees(crs, 1)
+end
+
+function axis_is_degrees(crs::Proj.CRS, axis_index::Int; context::Ptr{Proj.PJ_CONTEXT} = C_NULL)
+    @assert axis_index in (0, 1)
+    # In Proj, the `CoordinateSystem` object is contained within the `CRS` object.
+    # So we need to extract it explicitly, since Proj doesn't provide utilities for this.
+    cs = Proj.proj_crs_get_coordinate_system(crs.pj, context)
+
+    # Instantiate some refs that we'll use to get information out of the PJ struct
+    auth_name = Ref{Cstring}()
+    code = Ref{Cstring}()
+    unitname = Ref{Cstring}()
+
+    # Load unit info for the given axis into the pointers
+    Proj.proj_cs_get_axis_info(
+        cs, 
+        axis_index, 
+        C_NULL, # out_name
+        C_NULL, # out_abbrev
+        C_NULL, # out_direction
+        C_NULL, # out_unit_conv_factor
+        C_NULL, # out_unit_name
+        auth_name,
+        code,
+        context
+    )
+
+    # We don't `unsafe_string` the C strings, because we're just going to pass them to Proj's unit lookup function
+    Proj.proj_uom_get_info_from_database(
+        auth_name[],
+        code[],
+        unitname,   # out_name
+        C_NULL,     # out_conv_factor
+        C_NULL,     # out_category
+        context
+    )
+    # Destroy the coordinate system object
+    Proj.proj_destroy(cs)
+
+    unit_str = unsafe_string(unitname[])
+
+    # TODO: check if this is the correct way to check if the unit is degrees
+    # We can also check if the unit category is "angular", but I chose to replicate
+    # the original test from ArchGDAL here.
+    # If the unit is not "degree", we could still have an angular unit (radians or some linearly scaled thing),
+    # in which case we should technically return true.
+    # We'd also have to return the conversion factor in this case, and maybe the category (radians or degrees)...
+    return isequal(unit_str, "degree")
+end

ext/RastersProjExt/reproject.jl

@@ -0,0 +1,53 @@
+
+Rasters._reproject(source::GeoFormat, target::GeoFormat, dim::XDim, vals::AbstractVector) = _reproject!(source, target, dim, Float64.(vals))
+Rasters._reproject(source::GeoFormat, target::GeoFormat, dim::YDim, vals::AbstractVector) = _reproject!(source, target, dim, Float64.(vals))
+
+function _reproject!(source::GeoFormat, target::GeoFormat, ::XDim, vals::Vector{Float64}) 
+    # First, construct a transformation from `source` to `target`.
+    # TODO: add area of use, so that the transformation is more accurate.
+    trans = Proj.Transformation(source, target; always_xy = true, direction = Proj.PJ_FWD)
+    # Check that the reprojection is valid, i.e., trans((x, 0))[1] == trans((x, 1))[1]
+    first_x = first(vals)
+    zero_y = (first_x, zero(first_x))
+    one_y = (first_x, one(first_x))
+    # Check that the reprojection is valid, i.e., trans((x, 0))[1] == trans((x, 1))[1]
+    trans(zero_y)[1] == trans(one_y)[1] || _reproject_crs_error(source, target)
+    # Now that we've proved that the reprojection is valid, we can proceed.
+    # Here, for efficiency, and to avoid allocations, we'll use `proj_trans_generic`,
+    # which mutates values in-place.
+    Proj.proj_trans_generic(
+        trans.pj,                                    # pointer to transformation
+        Proj.PJ_FWD,                                 # direction of transformation
+        vals, sizeof(typeof(first_x)), length(vals), # input x values
+        C_NULL, 0, 0,                                # default (assumed 0) y values
+        C_NULL, 0, 0,                                # default (assumed 0) z values
+        C_NULL, 0, 0                                 # default (assumed 0) t values
+    )
+    return vals
+end
+function _reproject!(source::GeoFormat, target::GeoFormat, dim::YDim, vals::AbstractVector)
+    # First, construct a transformation from `source` to `target`.
+    # TODO: add area of use, so that the transformation is more accurate.
+    trans = Proj.Transformation(source, target; always_xy = true, direction = Proj.PJ_FWD)
+    # Check that the reprojection is valid, i.e., trans((0, y))[2] == trans((1, y))[2]
+    first_y = first(vals)
+    zero_x = (zero(first_y), first_y)
+    one_x = (one(first_y), first_y)
+    # Check that the reprojection is valid, i.e., trans((0, y))[2] == trans((1, y))[2]
+    trans(zero_x)[2] == trans(one_x)[2] || _reproject_crs_error(source, target)
+    # Now that we've proved that the reprojection is valid, we can proceed.
+    # Here, for efficiency, and to avoid allocations, we'll use `proj_trans_generic`,
+    # which mutates values in-place.
+    Proj.proj_trans_generic(
+        trans.pj,                                    # pointer to transformation
+        Proj.PJ_FWD,                                 # direction of transformation
+        C_NULL, 0, 0,                                # default (assumed 0) x values
+        vals, sizeof(typeof(first_y)), length(vals), # input y values
+        C_NULL, 0, 0,                                # default (assumed 0) z values
+        C_NULL, 0, 0                                 # default (assumed 0) t values
+    )
+    return vals
+end
+
+_reproject_crs_error(source, target) = 
+    throw(ArgumentError("Cannot reproject from: \n $source \nto: \n $target, coordinate reference systems are not aligned on all axes. You may need to use `resample` instead")) 

src/extensions.jl

@@ -156,6 +156,8 @@ $EXPERIMENTAL
 """
 warp(args...; kw...) = throw_extension_error(warp, "ArchGDAL", :RastersArchGDALExt, args)
 
+# stubs that need Proj
+
 """
     cellarea([method], x)
 
@@ -211,7 +213,7 @@ function cellarea(method::GeometryOpsCore.Spherical, dims::Tuple{<:XDim, <:YDim}
     return Raster(areas; dims)
 end
 
-_spherical_cellarea(args...; kw...) = throw_extension_error(_spherical_cellarea, "ArchGDAL", :RastersArchGDALExt, args)
+_spherical_cellarea(args...; kw...) = throw_extension_error(_spherical_cellarea, "Proj", :RastersProjExt, args)
 
 function _planar_cellarea(dims::Tuple{<:XDim, <:YDim})
     xbnds, ybnds = DD.intervalbounds(dims)

src/methods/reproject.jl

@@ -74,7 +74,7 @@ function _reproject(source::Nothing, target::GeoFormat, dim::Union{XDim,YDim}, v
     reshape(reproject(source, target, dim, vec(vals)), size(vals))
 end
 function _reproject(source::GeoFormat, target::GeoFormat, dim::Union{XDim,YDim}, vals::AbstractVector) 
-    error("Rasters.jl requires backends to be loaded externally as of version 0.8. Run `using ArchGDAL` to use `reproject`")
+    throw_extension_error(_reproject, "Proj", :RastersProjExt, (source, target, dim, vals))
 end
 
 # Guess the step for arrays

test/cellarea.jl

@@ -1,4 +1,4 @@
-using Rasters, DimensionalData, Rasters.Lookups, ArchGDAL
+using Rasters, DimensionalData, Rasters.Lookups, Proj
 using Test
 using DimensionalData: @dim, YDim
 include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))