Breaking: get rid of index

ext/RastersArchGDALExt/RastersArchGDALExt.jl

@@ -19,7 +19,7 @@ using Rasters: GDALsource, AbstractProjected, AbstractRaster, AbstractRasterStac
     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, resample, warp, cellsize, nokw, isnokw, isnokwornothing
+import Rasters: reproject, resample, warp, cellsize, nokw, isnokw, isnokwornothing, combine
 
 import LinearAlgebra: dot, cross
 

src/Rasters.jl

@@ -69,7 +69,7 @@ export missingval, boolmask, missingmask, replace_missing, replace_missing!,
        resample, warp, zonal, crop, extend, trim, slice, combine, points,
        classify, classify!, mosaic, mosaic!, extract, rasterize, rasterize!,
        coverage, coverage!, setcrs, setmappedcrs, smapseries, cellsize, cellarea
-export crs, mappedcrs, mappedindex, mappedbounds, projectedindex, projectedbounds
+export crs, mappedcrs, mappedlookup, mappedbounds, projectedlookup, projectedbounds
 export reproject, convertlookup
 export Extent, extent
 

src/array.jl

@@ -67,7 +67,7 @@ ismem(A::AbstractRaster) = !isdisk(A)
 function Base.:(==)(A::AbstractRaster{T,N}, B::AbstractRaster{T,N}) where {T,N}
     size(A) == size(B) && all(A .== B)
 end
-for f in (:mappedbounds, :projectedbounds, :mappedindex, :projectedindex)
+for f in (:mappedbounds, :projectedbounds, :mappedlookup, :projectedlookup)
     @eval ($f)(A::AbstractRaster, dims_) = ($f)(dims(A, dims_))
     @eval ($f)(A::AbstractRaster) = ($f)(dims(A))
 end

src/lookup.jl

@@ -180,7 +180,7 @@ convertlookup(::Type, lookup::Lookup) = lookup
 convertlookup(::Type{T1}, lookup::T2) where {T1,T2<:T1} = lookup
 # Otherwise AbstractProjected needs ArchGDAL
 function convertlookup(::Type{<:Mapped}, l::Projected)
-    newindex = reproject(crs(l), mappedcrs(l), dim(l), index(l))
+    newindex = reproject(crs(l), mappedcrs(l), dim(l), val(l))
     # We use Explicit mode and make a bounds matrix
     # This way the bounds can be saved correctly to NetCDF
     span = if sampling(l) isa Points
@@ -254,30 +254,16 @@ projectedbounds(crs::GeoFormat, lookup::Mapped, dim) =
 _sort((a, b)) = a <= b ? (a, b) : (b, a)
 
 """
-    mappedindex(x)
+    mappedlookup(x)
 
 Get the index value of a dimension converted to the `mappedcrs` value.
 
 Without ArchGDAL loaded, this is just the regular dim value.
 """
-function mappedindex end
+function mappedlookup end
 
-mappedindex(dims::Tuple) = map(mappedindex, dims)
-mappedindex(dim::Dimension) = _mappedindex(parent(dim), dim)
+mappedlookup(dims::Tuple) = map(mappedlookup, dims)
+mappedlookup(dim::Dimension) = reproject(mappedcrs(dim), lookup(dim))
 
-_mappedindex(::Lookup, dim::Dimension) = index(dim)
-_mappedindex(lookup::Projected, dim::Dimension) = _mappedindex(mappedcrs(lookup), lookup, dim)
-_mappedindex(mappedcrs::Nothing, lookup::Projected, dim) =
-    error("No mappedcrs attached to $(name(dim)) dimension")
-_mappedindex(mappedcrs::GeoFormat, lookup::Projected, dim) =
-    reproject(crs(dim), mappedcrs, dim, index(dim))
-
-projectedindex(dims::Tuple) = map(projectedindex, dims)
-projectedindex(dim::Dimension) = _projectedindex(parent(dim), dim)
-
-_projectedindex(::Lookup, dim::Dimension) = index(dim)
-_projectedindex(lookup::Mapped, dim::Dimension) = _projectedindex(crs(lookup), lookup, dim)
-_projectedindex(crs::Nothing, lookup::Mapped, dim::Dimension) =
-    error("No projection crs attached to $(name(dim)) dimension")
-_projectedindex(crs::GeoFormat, lookup::Mapped, dim::Dimension) =
-    reproject(mappedcrs(dim), crs, dim, index(dim))
+projectedlookup(dims::Tuple) = map(projectedlookup, dims)
+projectedlookup(dim::Dimension) = reproject(crs(dim), lookup(dim))

src/plotrecipes.jl

@@ -73,16 +73,16 @@ end
 
 
     if eltype(A) <: ColorTypes.Colorant
-        index(x), index(y), parent(A)
+        parent(lookup(x)), parent(lookup(y)), parent(A)
     elseif get(plotattributes, :seriestype, :none) == :contourf
         A = replace_missing(A, missing)
         clims = extrema(skipmissing(A))
         :levels --> range(clims[1], clims[2], length=20)
-        index(x), index(y), clamp.(A, clims[1], clims[2])
+        parent(lookup(x)), parent(lookup(y)), clamp.(A, clims[1], clims[2])
     else
         :seriestype --> :heatmap
         A = replace_missing(A, missing)
-        index(x), index(y), parent(A)
+        parent(lookup(x)), parent(lookup(y)), parent(A)
     end
 end
 
@@ -97,7 +97,7 @@ end
     :yguide --> yguide
     :label --> ""
     z = map(_prepare_plots, dims(A))
-    parent(A), index(z)
+    parent(A), parent.(lookup(z))
 end
 
 # Plot 3d arrays as multiple tiled plots
@@ -203,7 +203,7 @@ end
         end
     else
         thinned, plotinds, nplots = _maybe_thin_plots(A)
-        titles = string.(index(A, dims(A, 1)))
+        titles = string.(lookup(A, dims(A, 1)))
         ncols, nrows = _balance_grid(nplots)
         :layout --> (ncols, nrows)
         for r in 1:nrows, c in 1:ncols

src/sources/commondatamodel.jl

@@ -604,6 +604,6 @@ function _def_dim_var!(ds::AbstractDataset, dim::Dimension)
         push!(attrib, "bounds" => boundskey)
         CDM.defVar(ds, boundskey, bounds, ("bnds", dimname))
     end
-    CDM.defVar(ds, dimname, Vector(index(dim)), (dimname,); attrib=attrib)
+    CDM.defVar(ds, dimname, Vector(lookup(dim)), (dimname,); attrib=attrib)
     return nothing
 end

test/aggregate.jl

@@ -47,17 +47,17 @@ series = RasterSeries([stack1, stack2], (Ti(dates),))
     @test step(lookup(aglon)) === 30.0
     @test val(aglon) == [30.0]
     disaglon = disaggregate(Start(), aglon, 3)
-    @test index(disaglon) == index(dimz[1])
+    @test lookup(disaglon) == lookup(dimz[1])
     @test span(disaglon) == span(dimz[1])
     @test sampling(disaglon) == sampling(dimz[1])
 
     aglat = aggregate(Start(), dimz[2], 3)
     @test step(lookup(aglat)) === 15.0
-    @test index(aglat) == LinRange(-10.0, 5.0, 2)
+    @test lookup(aglat) == LinRange(-10.0, 5.0, 2)
     disaglat = disaggregate(Start(), aglat, 3)
     # The last item is lost due to rounding in `aggregate`
-    @test index(disaglat) != index(dimz[2])
-    @test index(disaglat) === LinRange(-10.0, 15.0, 6)
+    @test lookup(disaglat) != lookup(dimz[2])
+    @test lookup(disaglat) == LinRange(-10.0, 15.0, 6)
     @test span(disaglat) == span(dimz[2])
     @test sampling(disaglat) == sampling(dimz[2])
 end

test/array.jl

@@ -49,7 +49,7 @@ end
 end
 
 @testset "array dims have been formatted" begin
-    @test index(ga2) == (10.0:10:100, -50.0:10:50.0, [DateTime(2019)])
+    @test map(parent, lookup(ga2)) == (10.0:10:100, -50.0:10:50.0, [DateTime(2019)])
     @test dims(ga1)[1:2] == dims(ga2)[1:2]
 end
 

test/methods.jl

@@ -177,10 +177,10 @@ end
     mm_st2 = missingmask(st2)
     mm_st2_inverted = missingmask(st2; invert=true)
     @test dims(mm_st2) == dims(mm_st2_inverted) == dims(st2)
-    @test all(mm_st2 .=== [missing missing; true missing])    
-    @test all(mm_st2_inverted .=== [true true; missing true])    
-    @test all(missingmask(st2, alllayers = false) .=== [missing; true])    
-    @test all(missingmask(se) .=== missingmask(ga))
+    @test isequal(mm_st2, [missing missing; true missing])    
+    @test isequal(mm_st2_inverted, [true true; missing true])    
+    @test isequal(missingmask(st2, alllayers = false), [missing; true])    
+    @test isequal(missingmask(se), missingmask(ga))
     @test missingmask(polygon; res=1.0, boundary=:touches) == 
         fill!(Raster{Union{Missing,Bool}}(undef, X(Projected(-20:1.0:0.0; crs=nothing)), Y(Projected(10.0:1.0:30.0; crs=nothing))), true)
     x = missingmask([polygon, polygon]; collapse=false, res=1.0, boundary=:touches)
@@ -212,14 +212,14 @@ end
     @test all(mask(ga1; with=ga) .=== mask(ga2; with=ga) .=== [missing 1; 2 missing])
     @test all(mask(ga1; with=ga, invert=true) .=== mask(ga2; with=ga, invert=true) .=== [missing missing; missing 3])
     ga2 = replace_missing(ga1 .* 1.0; missingval=NaN)
-    @test all(mask(ga2; with=ga) .=== [NaN 1.0; 2.0 NaN])
-    @test all(mask(ga2; with=ga, invert=true) .=== [NaN NaN; NaN 3.0])
+    @test isequal(mask(ga2; with=ga), [NaN 1.0; 2.0 NaN])
+    @test isequal(mask(ga2; with=ga, invert=true), [NaN NaN; NaN 3.0])
     ga3 = replace_missing(ga1; missingval=-9999)
     mask!(ga3; with=ga)
-    @test all(ga3 .=== [-9999 1; 2 -9999])
+    @test isequal(ga3, [-9999 1; 2 -9999])
     ga4 = replace_missing(ga1; missingval=-9999)
     mask!(ga4; with=ga, invert=true)
-    @test all(ga4 .=== [-9999 -9999; -9999 3])
+    @test isequal(ga4, [-9999 -9999; -9999 3])
     maskfile = tempname() * ".tif"
     dmask = mask(ga3; with=ga, filename=maskfile)
     @test Rasters.isdisk(dmask)
@@ -382,9 +382,9 @@ end
 @testset "classify" begin
     A1 = [missing 1; 2 3]
     ga1 = Raster(A1, (X, Y); missingval=missing)
-    @test all(classify(ga1, 1=>99, 2=>88, 3=>77) .=== [missing 99; 88 77])
-    @test all(classify(ga1, 1=>99, 2=>88, 3=>77; others=0) .=== [missing 99; 88 77])
-    @test all(classify(ga1, 1=>99, 2=>88; others=0) .=== [missing 99; 88 0])
+    @test isequal(classify(ga1, 1=>99, 2=>88, 3=>77), [missing 99; 88 77])
+    @test isequal(classify(ga1, 1=>99, 2=>88, 3=>77; others=0), [missing 99; 88 77])
+    @test isequal(classify(ga1, 1=>99, 2=>88; others=0), [missing 99; 88 0])
     A2 = [1.0 2.5; 3.0 4.0]
     ga2 = Raster(A2 , (X, Y))
     @test classify(ga2, (2, 3)=>:x, >(3)=>:y) == [1.0 :x; 3.0 :y]
@@ -393,7 +393,7 @@ end
     @test ga2 == [1.0 0.0; -1.0 -1.0]
     classify!(ga2, [1 2.5 0.0; 2.5 4.0 -1.0]; lower=(>), upper=(<=))
     @test ga2 == [1.0 0.0; -1.0 -1.0]
-    @test all(classify(ga1, [1 99; 2 88; 3 77]) .=== [missing 99; 88 77])
+    @test isequal(classify(ga1, [1 99; 2 88; 3 77]), [missing 99; 88 77])
     @test_throws ArgumentError classify(ga1, [1, 2, 3])
 end
 
@@ -405,9 +405,9 @@ end
     table = (geometry=[missing, (9.0, 0.1), (9.0, 0.2), (10.0, 0.3)], foo=zeros(4))
     st = RasterStack(rast, rast2)
     ga = Raster(A, (X(9.0:1.0:10.0), Y(0.1:0.1:0.2)))
-    @test all(collect(points(ga; order=(Y, X))) .=== [missing (0.2, 9.0); (0.1, 10.0) missing])
-    @test all(collect(points(ga; order=(X, Y))) .=== [missing (9.0, 0.2); (10.0, 0.1) missing])
-    @test all(points(ga; order=(X, Y), ignore_missing=true) .===
+    @test isequal(collect(points(ga; order=(Y, X))), [missing (0.2, 9.0); (0.1, 10.0) missing])
+    @test isequal(collect(points(ga; order=(X, Y))), [missing (9.0, 0.2); (10.0, 0.1) missing])
+    @test isequal(points(ga; order=(X, Y), ignore_missing=true),
               [(9.0, 0.1) (9.0, 0.2); (10.0, 0.1) (10.0, 0.2)])
 end
 

test/reproject.jl

@@ -85,7 +85,7 @@ end
     A = DimArray(zeros(length(lon), length(lat)), (lon, lat))
     Aconv = convertlookup(Mapped, A)
 
-    @test index(Aconv) == (index(convertedlon), index(convertedlat))
+    @test lookup(Aconv) == (lookup(convertedlon), lookup(convertedlat))
     @test val.(span(Aconv)) == val.(span.((convertedlon, convertedlat)))
 end
 

test/resample.jl

@@ -30,15 +30,15 @@ raster_output = resample(cea; res=output_res, crs=output_crs, method, missingval
     @test missingval(resample(cea; res=output_res, crs=output_crs, method)) == 0x00
 end
 
-@testset "snapped size and dim index match" begin
+@testset "snapped size and dim lookup match" begin
     snaptarget = raster_output
     snapped = resample(cea; to=snaptarget)
     disk_snapped = resample(cea; to=snaptarget, filename="raster.tif")
     @test size(snapped) == size(disk_snapped) == size(snaptarget)
-    @test isapprox(index(snaptarget, Y), index(snapped, Y))
-    @test isapprox(index(snaptarget, X), index(snapped, X))
-    @test isapprox(index(snaptarget, Y), index(disk_snapped, Y))
-    @test isapprox(index(snaptarget, X), index(disk_snapped, X))
+    @test isapprox(lookup(snaptarget, Y), lookup(snapped, Y))
+    @test isapprox(lookup(snaptarget, X), lookup(snapped, X))
+    @test isapprox(lookup(snaptarget, Y), lookup(disk_snapped, Y))
+    @test isapprox(lookup(snaptarget, X), lookup(disk_snapped, X))
 end
 
 @testset "`method` only does nothing" begin

test/sources/gdal.jl

@@ -461,8 +461,8 @@ gdalpath = maybedownload(url)
             grdarray = Raster(fn)
             @test crs(grdarray) == convert(ProjString, crs(gdalarray))
             @test all(map((a, b) -> all(a .≈ b), bounds(grdarray), bounds(gdalarray)))
-            @test index(grdarray, Y) ≈ index(gdalarray, Y)
-            @test val(dims(grdarray, X)) ≈ val(dims(gdalarray, X))
+            @test lookup(grdarray, Y) ≈ lookup(gdalarray, Y)
+            @test lookup(grdarray, X) ≈ lookup(gdalarray, X)
             @test grdarray == gdalarray
         end
 
@@ -488,16 +488,16 @@ gdalpath = maybedownload(url)
             @test size(saved) == size(gdalarray)
             @test parent(saved) ≈ parent(gdalarray)
             clat, clon = DimensionalData.shiftlocus.(Ref(Center()), dims(gdalarray, (Y, X)))
-            @test index(clat) ≈ index(saved, Y)
-            @test index(clon) ≈ index(saved, X)
+            @test lookup(clat) ≈ lookup(saved, Y)
+            @test lookup(clon) ≈ lookup(saved, X)
             @test all(bounds(saved, X) .≈ bounds(clon))
             @test all(bounds(saved, Y) .≈ bounds(clat))
-            @test projectedindex(clon) ≈ projectedindex(saved, X)
+            @test projectedlookup(clon) ≈ projectedlookup(saved, X)
             @test all(projectedbounds(clon) .≈ projectedbounds(saved, X))
             # reason lat crs conversion is less accurate than lon TODO investigate further
             @test all(map((a, b) -> isapprox(a, b; rtol=1e-6),
-                projectedindex(gdalarray, Y),
-                projectedindex(DimensionalData.shiftlocus(Start(), dims(saved, Y)))
+                projectedlookup(gdalarray, Y),
+                projectedlookup(DimensionalData.shiftlocus(Start(), dims(saved, Y)))
             ))
             @test all(map((a, b) -> isapprox(a, b; rtol=1e-6), projectedbounds(saved, Y),  projectedbounds(gdalarray, Y)))
         end
@@ -616,7 +616,7 @@ gdalpath = maybedownload(url)
         @test order(dims(rast)) == (ForwardOrdered(), ForwardOrdered())
         @test span(rast) == (Regular(1.0), Regular(1.0))
         @test sampling(rast) == (Intervals(Start()), Intervals(Start()))
-        @test index(rast) == (range(; start=0.0, stop=239.0, length=240), range(start=0.0, stop=179.0, length=180))
+        @test lookup(rast) == (range(; start=0.0, stop=239.0, length=240), range(start=0.0, stop=179.0, length=180))
     end
 
 end
@@ -869,24 +869,24 @@ end
     rm("resample_a.tif")
     rm("resample_b.tif")
 
-    @testset "snapped size and dim index match" begin
+    @testset "snapped size and dim lookup match" begin
         snaptarget = aggregate(Center(), read(gdalarray), 2)
         snapped = resample(read(gdalarray); to=snaptarget)
         disk_snapped = resample(gdalarray; to=snaptarget, filename="snap_resample.tif")
         stack_snapped = resample(read(gdalstack); to=snaptarget, filename="snap_resample.tif")
         ser_snapped = resample(read(gdalser); to=snaptarget)
         extradim_snapped = resample(extradim_raster; to=snaptarget)
         @test size(snapped) == size(disk_snapped) == size(snaptarget)
-        @test isapprox(index(snaptarget, Y), index(snapped, Y))
-        @test isapprox(index(snaptarget, X), index(snapped, X))
-        @test isapprox(index(snaptarget, Y), index(stack_snapped, Y))
-        @test isapprox(index(snaptarget, X), index(stack_snapped, X))
-        @test isapprox(index(snaptarget, Y), index(first(ser_snapped), Y))
-        @test isapprox(index(snaptarget, X), index(first(ser_snapped), X))
-        @test isapprox(index(snaptarget, Y), index(disk_snapped, Y))
-        @test isapprox(index(snaptarget, X), index(disk_snapped, X))
-        @test isapprox(index(snaptarget, Y), index(extradim_snapped, Y))
-        @test isapprox(index(snaptarget, X), index(extradim_snapped, X))
+        @test isapprox(lookup(snaptarget, Y), lookup(snapped, Y))
+        @test isapprox(lookup(snaptarget, X), lookup(snapped, X))
+        @test isapprox(lookup(snaptarget, Y), lookup(stack_snapped, Y))
+        @test isapprox(lookup(snaptarget, X), lookup(stack_snapped, X))
+        @test isapprox(lookup(snaptarget, Y), lookup(first(ser_snapped), Y))
+        @test isapprox(lookup(snaptarget, X), lookup(first(ser_snapped), X))
+        @test isapprox(lookup(snaptarget, Y), lookup(disk_snapped, Y))
+        @test isapprox(lookup(snaptarget, X), lookup(disk_snapped, X))
+        @test isapprox(lookup(snaptarget, Y), lookup(extradim_snapped, Y))
+        @test isapprox(lookup(snaptarget, X), lookup(extradim_snapped, X))
         rm("snap_resample.tif")
         rm("snap_resample_a.tif")
         rm("snap_resample_b.tif")

test/sources/grd.jl

@@ -95,7 +95,7 @@ grdpath = stem * ".gri"
     @testset "custom keywords" begin
         customgrdarray = Raster(grdpath; 
             name=:test, crs=EPSG(1000), mappedcrs=EPSG(4326), refdims=(Ti(),),
-            write=true, lazy=true, dropband=false, replace_missing=true,
+            write=true, lazy=true, dropband=false,
         )
         @test name(customgrdarray) == :test
         @test refdims(customgrdarray) == (Ti(),)
@@ -259,8 +259,8 @@ grdpath = stem * ".gri"
             saved = Raster(filename2; crs=crs(grdarray))
             @test size(saved) == size(grdarray[Band(1)])
             @test all(parent(saved) .≈ parent(grdarray[Band(1)]))
-            @test index(saved, X) ≈ index(grdarray, X) .+ 0.5
-            @test index(saved, Y) ≈ index(grdarray, Y) .+ 0.5
+            @test lookup(saved, X) ≈ lookup(grdarray, X) .+ 0.5
+            @test lookup(saved, Y) ≈ lookup(grdarray, Y) .+ 0.5
             @test bounds(saved, Y) == bounds(grdarray, Y)
             @test bounds(saved, X) == bounds(grdarray, X)
             write(filename2, grdarray; force = true)

test/sources/gribdatasets.jl

@@ -148,11 +148,11 @@ v = ds[:z]
             ser = Rasters.slice(gribarray, Ti) 
             @test ser isa RasterSeries
             @test size(ser) == (4,)
-            @test index(ser, Ti) == DateTime(2017, 1, 1):Hour(12):DateTime(2017, 1, 2, 12)
+            @test lookup(ser, Ti) == DateTime(2017, 1, 1):Hour(12):DateTime(2017, 1, 2, 12)
             @test Rasters.bounds(ser) == ((DateTime(2017, 1, 1), DateTime(2017, 1, 2, 12)),)
             A = ser[1]
-            @test index(A, Y) == 90.0:-3.0:-90.0
-            @test index(A, X) == 0.0:3.0:357.0
+            @test lookup(A, Y) == 90.0:-3.0:-90.0
+            @test lookup(A, X) == 0.0:3.0:357.0
         end
     end