Add linear transformation correction when using getindex with a range…

… with a stepsize other than 1. (#97)

Project.toml

@@ -1,7 +1,7 @@
 name = "GeoArrays"
 uuid = "2fb1d81b-e6a0-5fc5-82e6-8e06903437ab"
 authors = ["Maarten Pronk <[email protected]>"]
-version = "0.7.3"
+version = "0.7.4"
 
 [deps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"

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

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