diff --git a/Project.toml b/Project.toml index 54cf832..3b7d1b7 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GeoArrays" uuid = "2fb1d81b-e6a0-5fc5-82e6-8e06903437ab" authors = ["Maarten Pronk "] -version = "0.7.3" +version = "0.7.4" [deps] ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3" diff --git a/src/geoarray.jl b/src/geoarray.jl index 588bda6..23406b1 100644 --- a/src/geoarray.jl +++ b/src/geoarray.jl @@ -97,14 +97,10 @@ function Base.getindex(ga::GeoArray, i::AbstractRange, j::AbstractRange, k::Unio A = getindex(ga.A, i, j, k) x, y = first(i) - 1, first(j) - 1 t = ga.f(SVector(x, y)) - GeoArray(A, AffineMap(ga.f.linear, t), ga.crs) -end -function Base.getindex(ga::GeoArray, i::AbstractRange, j::AbstractRange) - A = getindex(ga.A, i, j, :) - x, y = first(i) - 1, first(j) - 1 - t = ga.f(SVector(x, y)) - GeoArray(A, AffineMap(ga.f.linear, t), ga.crs) + l = ga.f.linear * SMatrix{2,2}([step(i) 0; 0 step(j)]) + GeoArray(A, AffineMap(l, t), ga.crs) end +Base.getindex(ga::GeoArray, i::AbstractRange, j::AbstractRange) = Base.getindex(ga, i, j, :) Base.getindex(ga::GeoArray, I::Vararg{<:Integer,2}) = getindex(ga.A, I..., :) Base.getindex(ga::GeoArray, I::Vararg{<:Integer,3}) = getindex(ga.A, I...) @@ -169,8 +165,8 @@ See `indices` for the inverse function. function coords(ga::GeoArray, p::SVector{2,<:Integer}, strategy::AbstractStrategy) SVector{2}(ga.f(p .- strategy.offset)) end -coords(ga::GeoArray, p::Vector{<:Integer}, strategy::AbstractStrategy = Center()) = coords(ga, SVector{2}(p), strategy) -coords(ga::GeoArray, p::Tuple{<:Integer,<:Integer}, strategy::AbstractStrategy = Center()) = coords(ga, SVector{2}(p), strategy) +coords(ga::GeoArray, p::Vector{<:Integer}, strategy::AbstractStrategy=Center()) = coords(ga, SVector{2}(p), strategy) +coords(ga::GeoArray, p::Tuple{<:Integer,<:Integer}, strategy::AbstractStrategy=Center()) = coords(ga, SVector{2}(p), strategy) """ indices(ga::GeoArray, p::SVector{2,<:AbstractFloat}, strategy::AbstractStrategy) @@ -182,19 +178,19 @@ See `coords` for the inverse function. function indices(ga::GeoArray, p::SVector{2,<:AbstractFloat}, strategy::AbstractStrategy) round.(Int, inv(ga.f)(p) .+ strategy.offset)::SVector{2,Int} end -indices(ga::GeoArray, p::Vector{<:AbstractFloat}, strategy::AbstractStrategy = Center()) = indices(ga, SVector{2}(p), strategy) -indices(ga::GeoArray, p::Tuple{<:AbstractFloat,<:AbstractFloat}, strategy::AbstractStrategy = Center()) = indices(ga, SVector{2}(p), strategy) +indices(ga::GeoArray, p::Vector{<:AbstractFloat}, strategy::AbstractStrategy=Center()) = indices(ga, SVector{2}(p), strategy) +indices(ga::GeoArray, p::Tuple{<:AbstractFloat,<:AbstractFloat}, strategy::AbstractStrategy=Center()) = indices(ga, SVector{2}(p), strategy) # Generate coordinates for complete GeoArray -function coords(ga::GeoArray, strategy::AbstractStrategy = Center()) +function coords(ga::GeoArray, strategy::AbstractStrategy=Center()) (ui, uj) = size(ga)[1:2] extra = typeof(strategy) == Center ? 0 : 1 [coords(ga, SVector{2}(i, j), strategy) for i in 1:ui+extra, j in 1:uj+extra]::Matrix{SVector{2,Float64}} end # Generate coordinates for one dimension of a GeoArray -function coords(ga::GeoArray, dim::Symbol, strategy::AbstractStrategy = Center()) +function coords(ga::GeoArray, dim::Symbol, strategy::AbstractStrategy=Center()) if is_rotated(ga) error("This method cannot be used for a rotated GeoArray") end diff --git a/test/test_geoarray.jl b/test/test_geoarray.jl index 61ebeac..c7a64a6 100644 --- a/test/test_geoarray.jl +++ b/test/test_geoarray.jl @@ -6,6 +6,11 @@ using CoordinateTransformations @test length(x[1, 2]) == 5 @test size(x[1:3, 1:3]) == (3, 3, 5) @test size(x[1:3, 1:3, 1:3]) == (3, 3, 3) + + x = GeoArray(rand(10, 10, 2)) + xs = x[1:2:end, 1:2:end] + @test size(xs) == (5, 5, 2) + @test bbox(xs) == bbox(x) end @testset "Concrete" begin @@ -18,9 +23,9 @@ end @testset "Reading rasters" begin ga = GeoArrays.read(joinpath(testdatadir, "data/utmsmall.tif")) - @test bbox(ga) == (min_x = 440720.0, min_y = 3.74532e6, max_x = 446720.0, max_y = 3.75132e6) - @test bboxes(ga)[1] == (min_x = 440720.0, max_x = 440780.0, min_y = 3.75126e6, max_y = 3.75132e6) - @test bboxes(ga)[end] == (min_x = 446660.0, max_x = 446720.0, min_y = 3.74532e6, max_y = 3.74538e6) + @test bbox(ga) == (min_x=440720.0, min_y=3.74532e6, max_x=446720.0, max_y=3.75132e6) + @test bboxes(ga)[1] == (min_x=440720.0, max_x=440780.0, min_y=3.75126e6, max_y=3.75132e6) + @test bboxes(ga)[end] == (min_x=446660.0, max_x=446720.0, min_y=3.74532e6, max_y=3.74538e6) end @testset "Coords" begin @@ -34,7 +39,7 @@ end end @testset "GeoArray constructors" begin - x, y = range(4, stop = 8.0, length = 10), range(0, stop = 1, length = 9) + x, y = range(4, stop=8.0, length=10), range(0, stop=1, length=9) ga2 = GeoArray(rand(10, 9), x, y) ga2 = GeoArray(rand(10, 9), x, y, "") ga3 = GeoArray(rand(10, 9, 8), x, y) @@ -45,7 +50,7 @@ end for i in 1:length(x), j in 1:length(y) @test GeoArrays.coords(ga3, [i, j]) ≈ [x[i], y[j]] end - x, y = range(4, stop = 8.0, length = 11), range(0, stop = 1, length = 9) + x, y = range(4, stop=8.0, length=11), range(0, stop=1, length=9) @test_throws ErrorException GeoArray(rand(10, 9), x, y) end