From 99c2fd6dd43a8de0df59240e89f1f1abbb59daa8 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 7 Oct 2024 22:00:08 -0700 Subject: [PATCH 1/7] Switch `reproject` to use Proj instead of ArchGDAL MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This commit creates a new `reproject` method that uses Proj to transform an array of values in place using `proj_trans_generic`, only allocating a Float64 output array. It also performs the reproject check before actually transforming the array. This decreases runtime on my machine from 4.168 ms (ArchGDAL) to 665.187 μs (Proj), and allocations from 35167 allocs: 1.254 MiB (ArchGDAL) to 108 allocs: 42.859 KiB (Proj). --- Project.toml | 3 ++ ext/RastersArchGDALExt/RastersArchGDALExt.jl | 2 +- ext/RastersProjExt/RastersProjExt.jl | 32 ++++++++++++ ext/RastersProjExt/reproject.jl | 53 ++++++++++++++++++++ src/methods/reproject.jl | 2 +- 5 files changed, 90 insertions(+), 2 deletions(-) create mode 100644 ext/RastersProjExt/RastersProjExt.jl create mode 100644 ext/RastersProjExt/reproject.jl diff --git a/Project.toml b/Project.toml index e08c77bdb..fc045d451 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.3" RasterDataSources = "0.5.7, 0.6" RecipesBase = "0.7, 0.8, 1.0" Reexport = "0.2, 1.0" diff --git a/ext/RastersArchGDALExt/RastersArchGDALExt.jl b/ext/RastersArchGDALExt/RastersArchGDALExt.jl index bee739cbe..5923239ac 100644 --- a/ext/RastersArchGDALExt/RastersArchGDALExt.jl +++ b/ext/RastersArchGDALExt/RastersArchGDALExt.jl @@ -34,7 +34,7 @@ const LA = Lookups include("gdal_source.jl") include("cellarea.jl") -include("reproject.jl") +# include("reproject.jl") include("resample.jl") include("warp.jl") diff --git a/ext/RastersProjExt/RastersProjExt.jl b/ext/RastersProjExt/RastersProjExt.jl new file mode 100644 index 000000000..0d398b14a --- /dev/null +++ b/ext/RastersProjExt/RastersProjExt.jl @@ -0,0 +1,32 @@ +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("reproject.jl") +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/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 From 520723bacf875a9728c2ba12993a2905485b9c85 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 7 Oct 2024 22:00:51 -0700 Subject: [PATCH 2/7] Fix invocations to ArchGDAL GFT converts in tests --- test/reproject.jl | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/test/reproject.jl b/test/reproject.jl index e471361d8..035b0739a 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)) + 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]) @@ -49,7 +50,7 @@ end @testset "convertlookup" begin projcea = ProjString("+proj=cea") - proj4326 = convert(ProjString, EPSG(4326)) + proj4326 = convert(ProjString, Proj.CRS(EPSG(4326))) lonstart, lonend = 0.5, 179.5 cealonstart, cealonend = reproject(proj4326, projcea, X(), [lonstart, lonend]) From 11deb98b061f455c217e3a68c678a2373408902b Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 7 Oct 2024 22:01:20 -0700 Subject: [PATCH 3/7] Add `+type=crs` to all Proj-strings This is a weakness of this method compared to the old ArchGDAL method. Does this need a fix in Proj? --- ext/RastersMakieExt/plotrecipes.jl | 3 --- ext/RastersProjExt/cellarea.jl | 0 src/extensions.jl | 4 +++- test/reproject.jl | 4 ++-- 4 files changed, 5 insertions(+), 6 deletions(-) create mode 100644 ext/RastersProjExt/cellarea.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/cellarea.jl b/ext/RastersProjExt/cellarea.jl new file mode 100644 index 000000000..e69de29bb 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/test/reproject.jl b/test/reproject.jl index 035b0739a..4212e35fd 100644 --- a/test/reproject.jl +++ b/test/reproject.jl @@ -3,7 +3,7 @@ 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") + 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))) @@ -49,7 +49,7 @@ using Rasters: reproject, convertlookup end @testset "convertlookup" begin - projcea = ProjString("+proj=cea") + projcea = ProjString("+proj=cea +type=crs") proj4326 = convert(ProjString, Proj.CRS(EPSG(4326))) lonstart, lonend = 0.5, 179.5 From c1632a19c7a09dd3ba29f5d20e520216cde6b426 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 7 Oct 2024 23:05:48 -0700 Subject: [PATCH 4/7] Implement cellarea in Proj --- ext/RastersArchGDALExt/RastersArchGDALExt.jl | 2 - ext/RastersProjExt/RastersProjExt.jl | 1 + ext/RastersProjExt/cellarea.jl | 185 +++++++++++++++++++ 3 files changed, 186 insertions(+), 2 deletions(-) diff --git a/ext/RastersArchGDALExt/RastersArchGDALExt.jl b/ext/RastersArchGDALExt/RastersArchGDALExt.jl index 5923239ac..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/RastersProjExt/RastersProjExt.jl b/ext/RastersProjExt/RastersProjExt.jl index 0d398b14a..bb9892946 100644 --- a/ext/RastersProjExt/RastersProjExt.jl +++ b/ext/RastersProjExt/RastersProjExt.jl @@ -28,5 +28,6 @@ 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 index e69de29bb..9d9d506db 100644 --- a/ext/RastersProjExt/cellarea.jl +++ 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 From 276144b3ab554ac5535e71f56a6dafaa465d005a Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 7 Oct 2024 23:09:45 -0700 Subject: [PATCH 5/7] Add Proj to test project --- Project.toml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index fc045d451..d327fa190 100644 --- a/Project.toml +++ b/Project.toml @@ -94,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" @@ -101,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"] From 21e77e98989c85a561386095bf7543da2a3cb072 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 7 Oct 2024 23:15:57 -0700 Subject: [PATCH 6/7] Switch cellarea tests to Proj --- test/cellarea.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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")) From ded1ec8a1b30c70b66f6fd3b394b33814c972dbe Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Tue, 8 Oct 2024 09:27:23 -0700 Subject: [PATCH 7/7] Bump Proj compat to 1.7.2 to get the WKT support --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index d327fa190..97473a0d1 100644 --- a/Project.toml +++ b/Project.toml @@ -70,7 +70,7 @@ Missings = "0.4, 1" NCDatasets = "0.13, 0.14" OffsetArrays = "1" ProgressMeter = "1" -Proj = "1.3" +Proj = "1.7.2" RasterDataSources = "0.5.7, 0.6" RecipesBase = "0.7, 0.8, 1.0" Reexport = "0.2, 1.0"