diff --git a/Project.toml b/Project.toml index e08c77bdb..97473a0d1 100644 --- a/Project.toml +++ b/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.7.2" RasterDataSources = "0.5.7, 0.6" RecipesBase = "0.7, 0.8, 1.0" Reexport = "0.2, 1.0" @@ -91,6 +94,7 @@ 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" @@ -98,4 +102,4 @@ 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"] diff --git a/ext/RastersArchGDALExt/RastersArchGDALExt.jl b/ext/RastersArchGDALExt/RastersArchGDALExt.jl index bee739cbe..dfab89cf9 100644 --- a/ext/RastersArchGDALExt/RastersArchGDALExt.jl +++ b/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") diff --git a/ext/RastersMakieExt/plotrecipes.jl b/ext/RastersMakieExt/plotrecipes.jl index 602e22921..e64883996 100644 --- a/ext/RastersMakieExt/plotrecipes.jl +++ b/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 diff --git a/ext/RastersProjExt/RastersProjExt.jl b/ext/RastersProjExt/RastersProjExt.jl new file mode 100644 index 000000000..bb9892946 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/ext/RastersProjExt/cellarea.jl b/ext/RastersProjExt/cellarea.jl new file mode 100644 index 000000000..9d9d506db --- /dev/null +++ b/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 \ No newline at end of file diff --git a/ext/RastersProjExt/reproject.jl b/ext/RastersProjExt/reproject.jl new file mode 100644 index 000000000..52f9eb471 --- /dev/null +++ b/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")) diff --git a/src/extensions.jl b/src/extensions.jl index 36fde9df3..f56f5bd5d 100644 --- a/src/extensions.jl +++ b/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) diff --git a/src/methods/reproject.jl b/src/methods/reproject.jl index a4f417edb..183729255 100644 --- a/src/methods/reproject.jl +++ b/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 diff --git a/test/cellarea.jl b/test/cellarea.jl index 22c57200f..8f4c66ac3 100644 --- a/test/cellarea.jl +++ b/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")) diff --git a/test/reproject.jl b/test/reproject.jl index e471361d8..4212e35fd 100644 --- a/test/reproject.jl +++ b/test/reproject.jl @@ -1,13 +1,13 @@ -using Rasters, Test, ArchGDAL +using Rasters, Test, Proj using Rasters.Lookups, Rasters.Dimensions using Rasters: reproject, convertlookup @testset "reproject" begin - cea = ProjString("+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0") - wktcea = convert(WellKnownText, cea) - projcea = convert(ProjString, cea) - wkt4326 = convert(WellKnownText, EPSG(4326)) - proj4326 = convert(ProjString, EPSG(4326)) + cea = ProjString("+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 +type=crs") + wktcea = convert(WellKnownText, Proj.CRS(cea)) + projcea = convert(ProjString, Proj.CRS(cea)) + wkt4326 = convert(WellKnownText, Proj.CRS(EPSG(4326))) + proj4326 = convert(ProjString, Proj.CRS(EPSG(4326))) @test reproject(proj4326, cea, Y(), -30.0) ≈ -3.658789324855012e6 @@ -37,7 +37,8 @@ using Rasters: reproject, convertlookup # Dims have changed @test dims(reprojected_raster) == (x1, y1) # Raster has not changed - @test reprojected_raster == raster + # NOTE: this test only works because both rasters are in memory. + @test parent(reprojected_raster) == parent(raster) @test_throws ArgumentError reproject(cea, EPSG(32618), Y(), [-3.658789324855012e6]) @test_throws ArgumentError reproject(cea, EPSG(32618), X(), [-3.658789324855012e6]) @@ -48,8 +49,8 @@ using Rasters: reproject, convertlookup end @testset "convertlookup" begin - projcea = ProjString("+proj=cea") - proj4326 = convert(ProjString, EPSG(4326)) + projcea = ProjString("+proj=cea +type=crs") + proj4326 = convert(ProjString, Proj.CRS(EPSG(4326))) lonstart, lonend = 0.5, 179.5 cealonstart, cealonend = reproject(proj4326, projcea, X(), [lonstart, lonend])