Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fix stack keywords #470

Merged
merged 6 commits into from
Jul 29, 2023
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
add some tests
  • Loading branch information
rafaqz committed Jul 26, 2023
commit b518251c38cf004f5904ef97e463c5df4612c45e
89 changes: 48 additions & 41 deletions src/stack.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -197,7 +198,7 @@ function RasterStack(data::Tuple{Vararg{<:AbstractArray}}, dims::Tuple; name=not
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
Expand All @@ -224,17 +225,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;
source=nothing, name=nothing, keys=name, resize=nothing, lazy=false, dropband=true, kw...
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
Expand All @@ -247,20 +247,27 @@ function RasterStack(filename::AbstractString;
end
RasterStack(joinpath.(Ref(filename), filenames); lazy, kw...)
else
# Load as a single file
st = if haslayers(source)
_layer_raster(filename; kw...)
# 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
else
# layeresfrom (default: Band) dims acts as layers
# With bands actings as layers
RasterStack(Raster(filename; source); kw...)
end
# Maybe split the stack into separate arrays to remove extra dims.
if !(keys isa Nothing)
map(identity, st)
else
st
end
end

# Maybe read the lazy stack to memory
st1 = lazy ? st : read(st)

# 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
Expand All @@ -272,7 +279,7 @@ function RasterStack(filename::AbstractString;
end
end

function _layer_raster(filename;
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, kw...
Expand All @@ -292,6 +299,7 @@ function _layer_raster(filename;
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...
)
Expand All @@ -309,11 +317,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
Expand All @@ -327,6 +335,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
Expand Down Expand Up @@ -367,29 +397,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)...]))
Expand Down
29 changes: 28 additions & 1 deletion test/sources/gdal.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -512,6 +512,33 @@ 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 "methods" begin
@testset "mean" begin
means = map(A -> mean(parent(A); dims=2), gdalstack)
Expand Down
15 changes: 15 additions & 0 deletions test/sources/ncdatasets.jl
Original file line number Diff line number Diff line change
Expand Up @@ -364,6 +364,21 @@ 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 "load ncstack" begin
@test ncstack isa RasterStack
@test all(ismissing, missingval(ncstack))
Expand Down
Loading