fix stack keywords (#470)

* fix stack keywords

* add some tests

* fix resize

* fix cat test

* fix ncdatasets tests

* bugfix gdal tests

src/stack.jl

@@ -138,6 +138,7 @@ Load a file path or a `NamedTuple` of paths as a `RasterStack`, or convert argum
 # Keywords
 
 - `name`: Used as stack layer names when a `Tuple`, `Vector` or splat of `Raster` is passed in.
+ Has no effect when `NameTuple` is used - the `NamedTuple` keys are the layer names.
 - `metadata`: A `Dict` or `DimensionalData.Metadata` object.
 - `refdims`: `Tuple` of `Dimension` that the stack was sliced from.
 - `layersfrom`: `Dimension` to source stack layers from if the file is not already multi-layered.
@@ -166,7 +167,7 @@ function RasterStack(
  filenames::Union{AbstractArray{<:AbstractString},Tuple{<:AbstractString,Vararg}};
  name=map(filekey, filenames), keys=name, kw...
 )
- RasterStack(NamedTuple{Tuple(keys)}(Tuple(filenames)); kw...)
+ RasterStack(NamedTuple{cleankeys(Tuple(keys))}(Tuple(filenames)); kw...)
 end
 function RasterStack(filenames::NamedTuple{K,<:Tuple{<:AbstractString,Vararg}};
  crs=nothing, mappedcrs=nothing, source=nothing, lazy=false, dropband=true, kw...
@@ -192,12 +193,11 @@ function RasterStack(filenames::NamedTuple{K,<:Tuple{<:AbstractString,Vararg}};
 end
 # Multi Raster stack from Tuple of AbstractArray
 function RasterStack(data::Tuple{Vararg{<:AbstractArray}}, dims::Tuple; name=nothing, keys=name, kw...)
- isnothing(keys) && throw(ArgumentError("`name` or `keys` keyword must be a tuple of `Symbol`"))
  return RasterStack(NamedTuple{cleankeys(keys)}(data), dims; kw...)
 end
 # Multi Raster stack from NamedTuple of AbstractArray
 function RasterStack(data::NamedTuple{<:Any,<:Tuple{Vararg{<:AbstractArray}}}, dims::Tuple; kw...)
- # TODO: make this more sophisticated an match dimension length to axes?
+ # TODO: make this more sophisticated and match dimension length to axes?
  layers = map(data) do A
  Raster(A, dims[1:ndims(A)])
  end
@@ -224,20 +224,16 @@ function RasterStack(layers::NamedTuple{<:Any,<:Tuple{Vararg{<:AbstractRaster}}}
  layermetadata=map(DD.metadata, layers)
  missingval=map(Rasters.missingval, layers)
  return RasterStack(
- data, dims, refdims, layerdims, metadata,
- layermetadata, missingval
+ data, dims, refdims, layerdims, metadata, layermetadata, missingval
  )
 end
-# Single-file stack from a string
+# Stack from a String
 function RasterStack(filename::AbstractString;
- dims=nothing, refdims=(), metadata=nothing, crs=nothing, mappedcrs=nothing,
- layerdims=nothing, layermetadata=nothing, missingval=nothing,
- source=nothing, name=nothing, keys=name, layersfrom=nothing,
- resize=nothing, lazy=false, ext=nothing, dropband=true,
+ source=nothing, name=nothing, keys=name, lazy=false, dropband=true, kw...
 )
  source = isnothing(source) ? _sourcetype(filename) : _sourcetype(source)
  st = if isdir(filename)
- # Load a whole directory
+ # Load as a whole directory
  filenames = readdir(filename)
  length(filenames) > 0 || throw(ArgumentError("No files in directory $filename"))
  # Detect keys from names
@@ -248,35 +244,30 @@ function RasterStack(filename::AbstractString;
  else
  keys
  end
- RasterStack(joinpath.(Ref(filename), filenames); keys)
+ RasterStack(joinpath.(Ref(filename), filenames); lazy, kw...)
  else
+ # Load as a single file
  st = if haslayers(source)
- crs = defaultcrs(source, crs)
- mappedcrs = defaultmappedcrs(source, mappedcrs)
- data, field_kw = _open(filename; source) do ds
- dims = dims isa Nothing ? DD.dims(ds, crs, mappedcrs) : dims
- refdims = refdims == () || refdims isa Nothing ? () : refdims
- layerdims = layerdims isa Nothing ? DD.layerdims(ds) : layerdims
- metadata = metadata isa Nothing ? DD.metadata(ds) : metadata
- layermetadata = layermetadata isa Nothing ? DD.layermetadata(ds) : layermetadata
- missingval = missingval isa Nothing ? Rasters.missingval(ds) : missingval
- data = FileStack{source}(ds, filename; keys)
- data, (; dims, refdims, layerdims, metadata, layermetadata, missingval)
+ # With multiple named layers
+ l_st = _layer_stack(filename; source, name, keys, kw...)
+
+ # Maybe split the stack into separate arrays to remove extra dims.
+ if !(keys isa Nothing)
+ map(identity, l_st)
+ else
+ l_st
  end
- RasterStack(data; field_kw...)
- else
- # Band dims acts as layers
- RasterStack(Raster(filename; lazy); layersfrom)
- end
- # Maybe split the stack into separate arrays to remove extra dims.
- if !(keys isa Nothing)
- map(identity, st)
  else
- st
+ # With bands actings as layers
+ RasterStack(Raster(filename; source); kw...)
  end
  end
+
+ # Maybe read the lazy stack to memory
  st1 = lazy ? st : read(st)
- if hasdim(st1, Band()) && size(st1, Band()) < 2
+
+ # Maybe drop the Band dimension
+ if dropband && hasdim(st1, Band()) && size(st1, Band()) == 1
  if lazy
  return view(st1, Band(1)) # TODO fix dropdims in DiskArrays
  else
@@ -286,10 +277,31 @@ function RasterStack(filename::AbstractString;
  return st1
  end
 end
+
+function _layer_stack(filename;
+ dims=nothing, refdims=(), metadata=nothing, crs=nothing, mappedcrs=nothing,
+ layerdims=nothing, layermetadata=nothing, missingval=nothing,
+ source=nothing, name=nothing, keys=name, resize=nothing, kw...
+)
+ crs = defaultcrs(source, crs)
+ mappedcrs = defaultmappedcrs(source, mappedcrs)
+ data, field_kw = _open(filename; source) do ds
+ dims = dims isa Nothing ? DD.dims(ds, crs, mappedcrs) : dims
+ refdims = refdims == () || refdims isa Nothing ? () : refdims
+ layerdims = layerdims isa Nothing ? DD.layerdims(ds) : layerdims
+ metadata = metadata isa Nothing ? DD.metadata(ds) : metadata
+ layermetadata = layermetadata isa Nothing ? DD.layermetadata(ds) : layermetadata
+ missingval = missingval isa Nothing ? Rasters.missingval(ds) : missingval
+ data = FileStack{source}(ds, filename; keys)
+ data, (; dims, refdims, layerdims, metadata, layermetadata, missingval)
+ end
+ return RasterStack(data; field_kw..., kw...)
+end
+
+# Stack from a Raster
 function RasterStack(A::Raster;
  layersfrom=nothing, name=nothing, keys=name, metadata=metadata(A), refdims=refdims(A), kw...
 )
-
  keys = keys isa Union{AbstractString,Symbol,Name} ? (keys,) : keys
  layers = if isnothing(layersfrom)
  keys = if keys isa Nothing
@@ -303,11 +315,11 @@ function RasterStack(A::Raster;
  slices = slice(A, layersfrom)
  NamedTuple{cleankeys(keys)}(Tuple(slices))
  end
- RasterStack(layers; refdims=refdims, metadata=metadata, kw...)
+ return RasterStack(layers; refdims=refdims, metadata=metadata, kw...)
 end
-# Stack from stack, dims args
+# Stack from stack and dims args
 RasterStack(st::AbstractRasterStack, dims::Tuple; kw...) = RasterStack(st; dims, kw...)
-# Stack from table, dims args
+# Stack from table and dims args
 function RasterStack(table, dims::Tuple; name=_not_a_dimcol(table, dims), keys=name, kw...)
  # TODO use `name` everywhere, not keys
  if keys isa Symbol
@@ -321,6 +333,28 @@ function RasterStack(table, dims::Tuple; name=_not_a_dimcol(table, dims), keys=n
  end
  return RasterStack(layers, dims; kw...)
 end
+# Rebuild from internals
+function RasterStack(
+ data::Union{FileStack,OpenStack,NamedTuple{<:Any,<:Tuple{Vararg{<:AbstractArray}}}};
+ dims, refdims=(), layerdims, metadata=NoMetadata(), layermetadata, missingval)
+ return RasterStack(
+ data, dims, refdims, layerdims, metadata, layermetadata, missingval
+ )
+end
+# RasterStack from another stack
+function RasterStack(s::AbstractDimStack; name=cleankeys(Base.keys(s)), keys=name,
+ data=NamedTuple{keys}(s[key] for key in keys),
+ dims=dims(s), refdims=refdims(s), layerdims=DD.layerdims(s),
+ metadata=metadata(s), layermetadata=DD.layermetadata(s),
+ missingval=missingval(s)
+)
+ st = RasterStack(
+ data, DD.dims(s), refdims, layerdims, metadata, layermetadata, missingval
+ )
+
+ # TODO This is a bit of a hack, it should use `formatdims`.
+ return set(st, dims...)
+end
 
 function DD.modify(f, s::AbstractRasterStack{<:FileStack})
  open(s) do o
@@ -361,29 +395,6 @@ function _layerkeysfromdim(A, dim)
  end
 end
 
-# Rebuild from internals
-function RasterStack(
- data::Union{FileStack,OpenStack,NamedTuple{<:Any,<:Tuple{Vararg{<:AbstractArray}}}};
- dims, refdims=(), layerdims, metadata=NoMetadata(), layermetadata, missingval)
- st = RasterStack(
- data, dims, refdims, layerdims, metadata, layermetadata, missingval
- )
-end
-# RasterStack from another stack
-function RasterStack(s::AbstractDimStack; name=cleankeys(Base.keys(s)), keys=name,
- data=NamedTuple{keys}(s[key] for key in keys),
- dims=dims(s), refdims=refdims(s), layerdims=DD.layerdims(s),
- metadata=metadata(s), layermetadata=DD.layermetadata(s),
- missingval=missingval(s)
-)
- st = RasterStack(
- data, DD.dims(s), refdims, layerdims, metadata, layermetadata, missingval
- )
-
- # TODO This is a bit of a hack, it should use `formatdims`.
- return set(st, dims...)
-end
-
 Base.convert(::Type{RasterStack}, src::AbstractDimStack) = RasterStack(src)
 
 Raster(stack::RasterStack) = cat(values(stack)...; dims=Band([keys(stack)...]))

src/utils.jl

@@ -322,4 +322,3 @@ function _float64_xy_extent(ext::Extents.Extent)
  ybounds = map(Float64, ext.Y)
  return Extents.Extent(X=xbounds, Y=ybounds)
 end
-

test/sources/gdal.jl

@@ -72,7 +72,7 @@ gdalpath = maybedownload(url)
  end
 
  @testset "read and write band names" begin
- A = set(cat(gdalarray, gdalarray; dims=Band), Band=>1:2)
+ A = cat(gdalarray, gdalarray; dims=Band(1:2))
  named = set(A, Band => string.(Ref("layer_"), dims(A, Band)))
  tempfile = tempname() * ".tif"
  write(tempfile, named)
@@ -275,7 +275,7 @@ gdalpath = maybedownload(url)
  end
 
  @testset "3d, with subsetting" begin
- geoA2 = cat(gdalarray, gdalarray; dims=Band)[Y(4.224e6..4.226e6), X(-28492..0)]
+ geoA2 = cat(gdalarray, gdalarray; dims=Band(1:2))[Y(4.224e6..4.226e6), X(-28492..0)]
  geoA2 = set(geoA2, Band => Band(1:2))
  filename2 = tempname() * ".tif"
  write(filename2, geoA2)
@@ -433,7 +433,7 @@ gdalpath = maybedownload(url)
  am = AffineMap([60.0 20; 40 60], [first.(bounds(gdalarray, (X, Y)))...])
  xap = Rasters.AffineProjected(am; crs=crs(gdalarray), paired_lookup=parent(lookup(gdalarray, X)))
  yap = Rasters.AffineProjected(am; crs=crs(gdalarray), paired_lookup=parent(lookup(gdalarray, Y)))
- twoband = cat(gdalarray, gdalarray; dims=Band)
+ twoband = cat(gdalarray, gdalarray; dims=Band(1:2))
  affine_dims = DimensionalData.format((X(xap), Y(yap), Band(1:2)), twoband)
  rotated = rebuild(twoband; dims=affine_dims);
  @test occursin("Extent", sprint(show, MIME"text/plain"(), rotated))
@@ -512,6 +512,38 @@ end
  @test view(gdalstack, Y(2:3), X(1))[:a] == [0x00, 0x6b]
  end
 
+ @testset "lazy" begin
+ gdalstack_lazy = RasterStack((a=gdalpath, b=gdalpath); lazy=true)
+ @test Rasters.isdisk(gdalstack_lazy.a)
+ @test Rasters.isdisk(gdalstack_lazy.b)
+ gdalstack_eager = RasterStack((a=gdalpath, b=gdalpath); lazy=false)
+ @test !Rasters.isdisk(gdalstack_eager.a)
+ @test !Rasters.isdisk(gdalstack_eager.b)
+ end
+
+ @testset "dropband" begin
+ gdalstack_band = RasterStack((a=gdalpath, b=gdalpath); dropband=false)
+ @test hasdim(gdalstack_band, Band)
+ gdalstack_noband = RasterStack((a=gdalpath, b=gdalpath); dropband=true)
+ @test !hasdim(gdalstack_noband, Band)
+ end
+
+ @testset "source" begin
+ no_ext = tempname()
+ cp(gdalpath, no_ext)
+ a = RasterStack((a=no_ext, b=no_ext); source=:gdal)
+ b = RasterStack((a=no_ext, b=no_ext); source=Rasters.GDALsource())
+ # GDAL is the fallback anyway
+ c = RasterStack((a=no_ext, b=no_ext))
+ @test a == b == c == gdalstack
+ rm(no_ext)
+ end
+
+ @testset "name" begin
+ gdalstack_names = RasterStack((gdalpath, gdalpath); name=(:c, :d))
+ @test keys(gdalstack_names) == (:c, :d)
+ end
+
  @testset "methods" begin
  @testset "mean" begin
  means = map(A -> mean(parent(A); dims=2), gdalstack)

test/sources/ncdatasets.jl

@@ -364,6 +364,33 @@ end
  @test parent(eagerstack[:xi]) isa Array
  end
 
+ @testset "source" begin
+ no_ext = tempname()
+ cp(ncmulti, no_ext)
+ a = RasterStack(no_ext; source=:netcdf)
+ b = RasterStack(no_ext; source=Rasters.NCDsource())
+ @test a == b == ncstack
+ rm(no_ext)
+ end
+
+ @testset "crs" begin
+ st = RasterStack(ncmulti; crs=EPSG(3857), mappedcrs=EPSG(3857))
+ @test crs(st) == EPSG(3857)
+ @test mappedcrs(st) == EPSG(3857)
+ end
+
+ @testset "name" begin
+ @testset "multi name from single file" begin
+ @time small_stack = RasterStack(ncmulti; name=(:sofllac, :xlvi))
+ @test keys(small_stack) == (:sofllac, :xlvi)
+ end
+ @testset "multi file with single name" begin
+ tempnc = tempname() * ".nc"
+ write(tempnc, rebuild(Raster(ncsingle); name=:tos2))
+ @time small_stack = RasterStack((ncsingle, tempnc); name=(:tos, :tos2))
+ end
+ end
+
  @testset "load ncstack" begin
  @test ncstack isa RasterStack
  @test all(ismissing, missingval(ncstack))

test/stack.jl

@@ -132,13 +132,21 @@ end
  stack_b = RasterStack((ga1=Raster(data1 .+ 10, dims1b), ga2=Raster(data2 .+ 20, dims2b)))
  catstack = cat(stack_a, stack_b; dims=X())
  @test size(first(catstack)) == (20, 11)
+ @test size(last(catstack)) == (20, 11, 1)
  @test val(dims(catstack, X)) ≈ 10.0:10.0:200.0
  #@test step(dims(first(catstack), X())) == 10.0
  @test DimensionalData.bounds(dims(first(catstack), X)) == (10.0, 200.0)
  @test catstack[:ga1][Y(1)] == 1.0:20.0
  @test catstack[:ga2][Y(1), Ti(1)] == 2.0:2.0:40.0
- catstack = cat(stack_a, stack_b; dims=(X(), Y()))
+ dims2c = (dims1b..., Ti([DateTime(2019)]))
+ stack_c = set(stack_b, X=>110:10:200, Y=>60:10:160)
+ catstack = cat(stack_a, stack_c; dims=(X(), Y()))
  @test size(first(catstack)) == (20, 22)
+ @test size(last(catstack)) == (20, 22, 1)
+ @test dims(catstack) == 
+ (X(Sampled(10:10.0:200, ForwardOrdered(), Regular(10.0), Points(), NoMetadata())),
+ Y(Sampled(-50:10.0:160, ForwardOrdered(), Regular(10.0), Points(), NoMetadata())),
+ Ti(Sampled([DateTime(2019)], ForwardOrdered(), Irregular(nothing, nothing), Points(), NoMetadata())))
 end
 
 @testset "copy" begin