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

calling read on 0 dimensional raster fails #473

Closed
maxfreu opened this issue Jul 5, 2023 · 7 comments · Fixed by meggart/DiskArrays.jl#100
Closed

calling read on 0 dimensional raster fails #473

maxfreu opened this issue Jul 5, 2023 · 7 comments · Fixed by meggart/DiskArrays.jl#100

Comments

@maxfreu
Copy link
Contributor

maxfreu commented Jul 5, 2023

Hi! I would have expected read to work in the following case:

someraster = "..."
r = Raster(someraster; lazy=true)
r[X(Near(0)), Y(Near(0)), Band(1)]  # returns value
v = @view r[X(Near(0)), Y(Near(0)), Band(1)]  # returns view and prints nicely
read(v)  # errors with no method matching for copyto! Int16 something deep down in diskarrays

I just dug out old code and I think it used to work. Maybe it can be fixed with some code along the lines of:

function Rasters.read(r::Raster)
    if ndims(r) == 0
        return rebuild(r, first.(r))
    else
        return modify(Array, x)
    end
end

function Rasters.read(r::RasterSeries)
    if ndims(first(r)) == 0
        return rebuild(r, read.(r))
    else
        return modify(Array, x)
    end
end

# sth similar for rasterstack?
@rafaqz
Copy link
Owner

rafaqz commented Jul 5, 2023

Yeah that should work. Can you paste the error?

MWE is always nice too ;)

@maxfreu
Copy link
Contributor Author

maxfreu commented Jul 6, 2023

Sorry, I was quite tired yesterday and missed the essential. A minimum failing example is provided in the initial post. Here is the stack trace:

ERROR: MethodError: no method matching copyto!(::Int16, ::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}, Tuple{}, typeof(identity), Tuple{Int16}})

Closest candidates are:
  copyto!(::AbstractArray, ::Base.Broadcast.Broadcasted{<:Base.Broadcast.AbstractArrayStyle{0}})
   @ Base broadcast.jl:929
  copyto!(::AbstractArray, ::Base.Broadcast.Broadcasted)
   @ Base broadcast.jl:926
  copyto!(::AbstractArray, ::Any)
   @ Base abstractarray.jl:943

Stacktrace:
  [1] materialize!
    @ ./broadcast.jl:884 [inlined]
  [2] materialize!(dest::Int16, bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}, Nothing, typeof(identity), Tuple{Int16}})
    @ Base.Broadcast ./broadcast.jl:881
  [3] (::DiskArrays.var"#59#61"{Array{Int16, 0}, Base.Broadcast.Broadcasted{DiskArrays.ChunkStyle{0}, Tuple{}, typeof(identity), Tuple{DiskArrays.SubDiskArray{Int16, 0}}}})(cnow::Tuple{})
    @ DiskArrays ~/.julia/packages/DiskArrays/4dn4s/src/broadcast.jl:42
  [4] foreach(f::DiskArrays.var"#59#61"{Array{Int16, 0}, Base.Broadcast.Broadcasted{DiskArrays.ChunkStyle{0}, Tuple{}, typeof(identity), Tuple{DiskArrays.SubDiskArray{Int16, 0}}}}, itr::DiskArrays.GridChunks{0})
    @ Base ./abstractarray.jl:3073
  [5] copyto!(dest::Array{Int16, 0}, bc::Base.Broadcast.Broadcasted{DiskArrays.ChunkStyle{0}, Tuple{}, typeof(identity), Tuple{DiskArrays.SubDiskArray{Int16, 0}}})
    @ DiskArrays ~/.julia/packages/DiskArrays/4dn4s/src/broadcast.jl:38
  [6] materialize!
    @ ./broadcast.jl:884 [inlined]
  [7] materialize!
    @ ./broadcast.jl:881 [inlined]
  [8] _Array(a::DiskArrays.SubDiskArray{Int16, 0})
    @ DiskArrays ~/.julia/packages/DiskArrays/4dn4s/src/array.jl:42
  [9] Array(a::DiskArrays.SubDiskArray{Int16, 0})
    @ DiskArrays ~/.julia/packages/DiskArrays/4dn4s/src/array.jl:5
 [10] (::Rasters.var"#31#32"{UnionAll})(O::Raster{Int16, 0, Tuple{}, Tuple{X{Projected{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Float64}, DimensionalData.Dimensions.LookupArrays.Intervals{DimensionalData.Dimensions.LookupArrays.Start}, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, WellKnownText{GeoFormatTypes.CRS}, Nothing, X{Colon}}}, Y{Projected{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.LookupArrays.ReverseOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Float64}, DimensionalData.Dimensions.LookupArrays.Intervals{DimensionalData.Dimensions.LookupArrays.Start}, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, WellKnownText{GeoFormatTypes.CRS}, Nothing, Y{Colon}}}, Band{DimensionalData.Dimensions.LookupArrays.Categorical{Int64, UnitRange{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}, DiskArrays.SubDiskArray{Int16, 0}, Symbol, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, Int16})
    @ Rasters ~/.julia/dev/Rasters/src/array.jl:128
 [11] (::Rasters.var"#39#40"{Rasters.var"#31#32"{UnionAll}, Raster{Int16, 0, Tuple{}, Tuple{X{Projected{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Float64}, DimensionalData.Dimensions.LookupArrays.Intervals{DimensionalData.Dimensions.LookupArrays.Start}, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, WellKnownText{GeoFormatTypes.CRS}, Nothing, X{Colon}}}, Y{Projected{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.LookupArrays.ReverseOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Float64}, DimensionalData.Dimensions.LookupArrays.Intervals{DimensionalData.Dimensions.LookupArrays.Start}, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, WellKnownText{GeoFormatTypes.CRS}, Nothing, Y{Colon}}}, Band{DimensionalData.Dimensions.LookupArrays.Categorical{Int64, UnitRange{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}, DiskArrays.SubDiskArray{Int16, 0}, Symbol, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, Int16}})(x::ArchGDAL.RasterDataset{Int16, ArchGDAL.Dataset})
    @ Rasters ~/.julia/dev/Rasters/src/array.jl:203
 [12] call_composed
    @ ./operators.jl:1035 [inlined]
 [13] call_composed
    @ ./operators.jl:1034 [inlined]
 [14] (::ComposedFunction{typeof(Rasters.cleanreturn), Rasters.var"#39#40"{Rasters.var"#31#32"{UnionAll}, Raster{Int16, 0, Tuple{}, Tuple{X{Projected{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Float64}, DimensionalData.Dimensions.LookupArrays.Intervals{DimensionalData.Dimensions.LookupArrays.Start}, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, WellKnownText{GeoFormatTypes.CRS}, Nothing, X{Colon}}}, Y{Projected{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.LookupArrays.ReverseOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Float64}, DimensionalData.Dimensions.LookupArrays.Intervals{DimensionalData.Dimensions.LookupArrays.Start}, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, WellKnownText{GeoFormatTypes.CRS}, Nothing, Y{Colon}}}, Band{DimensionalData.Dimensions.LookupArrays.Categorical{Int64, UnitRange{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}, DiskArrays.SubDiskArray{Int16, 0}, Symbol, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, Int16}}})(x::ArchGDAL.RasterDataset{Int16, ArchGDAL.Dataset}; kw::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Base ./operators.jl:1031
 [15] (::ComposedFunction{typeof(Rasters.cleanreturn), Rasters.var"#39#40"{Rasters.var"#31#32"{UnionAll}, Raster{Int16, 0, Tuple{}, Tuple{X{Projected{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Float64}, DimensionalData.Dimensions.LookupArrays.Intervals{DimensionalData.Dimensions.LookupArrays.Start}, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, WellKnownText{GeoFormatTypes.CRS}, Nothing, X{Colon}}}, Y{Projected{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.LookupArrays.ReverseOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Float64}, DimensionalData.Dimensions.LookupArrays.Intervals{DimensionalData.Dimensions.LookupArrays.Start}, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, WellKnownText{GeoFormatTypes.CRS}, Nothing, Y{Colon}}}, Band{DimensionalData.Dimensions.LookupArrays.Categorical{Int64, UnitRange{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}, DiskArrays.SubDiskArray{Int16, 0}, Symbol, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, Int16}}})(x::ArchGDAL.RasterDataset{Int16, ArchGDAL.Dataset})
    @ Base ./operators.jl:1031
 [16] readraster(f::ComposedFunction{typeof(Rasters.cleanreturn), Rasters.var"#39#40"{Rasters.var"#31#32"{UnionAll}, Raster{Int16, 0, Tuple{}, Tuple{X{Projected{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Float64}, DimensionalData.Dimensions.LookupArrays.Intervals{DimensionalData.Dimensions.LookupArrays.Start}, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, WellKnownText{GeoFormatTypes.CRS}, Nothing, X{Colon}}}, Y{Projected{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.LookupArrays.ReverseOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Float64}, DimensionalData.Dimensions.LookupArrays.Intervals{DimensionalData.Dimensions.LookupArrays.Start}, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, WellKnownText{GeoFormatTypes.CRS}, Nothing, Y{Colon}}}, Band{DimensionalData.Dimensions.LookupArrays.Categorical{Int64, UnitRange{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}, DiskArrays.SubDiskArray{Int16, 0}, Symbol, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, Int16}}}, args::String; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ArchGDAL ~/.julia/packages/ArchGDAL/R3wJR/src/context.jl:267
 [17] readraster(f::Function, args::String)
    @ ArchGDAL ~/.julia/packages/ArchGDAL/R3wJR/src/context.jl:264
 [18] _open(f::Function, ::Type{Rasters.GDALsource}, filename::String; write::Bool, kw::Base.Pairs{Symbol, Nothing, Tuple{Symbol}, NamedTuple{(:key,), Tuple{Nothing}}})
    @ RastersArchGDALExt ~/.julia/dev/Rasters/ext/RastersArchGDALExt/gdal_source.jl:136
 [19] _open
    @ ~/.julia/dev/Rasters/ext/RastersArchGDALExt/gdal_source.jl:122 [inlined]
 [20] #open#21
    @ ~/.julia/dev/Rasters/src/filearray.jl:40 [inlined]
 [21] open
    @ ~/.julia/dev/Rasters/src/filearray.jl:39 [inlined]
 [22] #_open_one#38
    @ ~/.julia/dev/Rasters/src/array.jl:199 [inlined]
 [23] _open_one(f::Function, A::Raster{Int16, 0, Tuple{}, Tuple{X{Projected{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Float64}, DimensionalData.Dimensions.LookupArrays.Intervals{DimensionalData.Dimensions.LookupArrays.Start}, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, WellKnownText{GeoFormatTypes.CRS}, Nothing, X{Colon}}}, Y{Projected{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.LookupArrays.ReverseOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Float64}, DimensionalData.Dimensions.LookupArrays.Intervals{DimensionalData.Dimensions.LookupArrays.Start}, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, WellKnownText{GeoFormatTypes.CRS}, Nothing, Y{Colon}}}, Band{DimensionalData.Dimensions.LookupArrays.Categorical{Int64, UnitRange{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}, DiskArrays.SubDiskArray{Int16, 0}, Symbol, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, Int16}, fa::Rasters.FileArray{Rasters.GDALsource, Int16, 3, Nothing, DiskArrays.GridChunks{3}, DiskArrays.Chunked})
    @ Rasters ~/.julia/dev/Rasters/src/array.jl:198
 [24] open(f::Rasters.var"#31#32"{UnionAll}, A::Raster{Int16, 0, Tuple{}, Tuple{X{Projected{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Float64}, DimensionalData.Dimensions.LookupArrays.Intervals{DimensionalData.Dimensions.LookupArrays.Start}, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, WellKnownText{GeoFormatTypes.CRS}, Nothing, X{Colon}}}, Y{Projected{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.LookupArrays.ReverseOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Float64}, DimensionalData.Dimensions.LookupArrays.Intervals{DimensionalData.Dimensions.LookupArrays.Start}, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, WellKnownText{GeoFormatTypes.CRS}, Nothing, Y{Colon}}}, Band{DimensionalData.Dimensions.LookupArrays.Categorical{Int64, UnitRange{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}, DiskArrays.SubDiskArray{Int16, 0}, Symbol, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, Int16}; kw::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Rasters ~/.julia/dev/Rasters/src/array.jl:191
 [25] open(f::Function, A::Raster{Int16, 0, Tuple{}, Tuple{X{Projected{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Float64}, DimensionalData.Dimensions.LookupArrays.Intervals{DimensionalData.Dimensions.LookupArrays.Start}, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, WellKnownText{GeoFormatTypes.CRS}, Nothing, X{Colon}}}, Y{Projected{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.LookupArrays.ReverseOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Float64}, DimensionalData.Dimensions.LookupArrays.Intervals{DimensionalData.Dimensions.LookupArrays.Start}, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, WellKnownText{GeoFormatTypes.CRS}, Nothing, Y{Colon}}}, Band{DimensionalData.Dimensions.LookupArrays.Categorical{Int64, UnitRange{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}, DiskArrays.SubDiskArray{Int16, 0}, Symbol, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, Int16})
    @ Rasters ~/.julia/dev/Rasters/src/array.jl:184
 [26] modify(f::Type{Array}, A::Raster{Int16, 0, Tuple{}, Tuple{X{Projected{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Float64}, DimensionalData.Dimensions.LookupArrays.Intervals{DimensionalData.Dimensions.LookupArrays.Start}, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, WellKnownText{GeoFormatTypes.CRS}, Nothing, X{Colon}}}, Y{Projected{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.LookupArrays.ReverseOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Float64}, DimensionalData.Dimensions.LookupArrays.Intervals{DimensionalData.Dimensions.LookupArrays.Start}, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, WellKnownText{GeoFormatTypes.CRS}, Nothing, Y{Colon}}}, Band{DimensionalData.Dimensions.LookupArrays.Categorical{Int64, UnitRange{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}, DiskArrays.SubDiskArray{Int16, 0}, Symbol, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, Int16})
    @ Rasters ~/.julia/dev/Rasters/src/array.jl:127
 [27] read(x::Raster{Int16, 0, Tuple{}, Tuple{X{Projected{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Float64}, DimensionalData.Dimensions.LookupArrays.Intervals{DimensionalData.Dimensions.LookupArrays.Start}, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, WellKnownText{GeoFormatTypes.CRS}, Nothing, X{Colon}}}, Y{Projected{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.LookupArrays.ReverseOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Float64}, DimensionalData.Dimensions.LookupArrays.Intervals{DimensionalData.Dimensions.LookupArrays.Start}, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, WellKnownText{GeoFormatTypes.CRS}, Nothing, Y{Colon}}}, Band{DimensionalData.Dimensions.LookupArrays.Categorical{Int64, UnitRange{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}, DiskArrays.SubDiskArray{Int16, 0}, Symbol, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, Int16})
    @ Rasters ~/.julia/dev/Rasters/src/read.jl:9
 [28] top-level scope
    @ ~/projects/force_cutouts/sandbox.jl:69

To me this roughly looks like "help, I can't write an array object into an int". That's why I'm manually rebuilding the raster in the suggested fix.

@maxfreu
Copy link
Contributor Author

maxfreu commented Jul 6, 2023

The only thing I'm unsure about is whether the call to read should return the Raster with all the refdims, or just the value.

Edit: Actually, I think reading a view into a raster should behave like the direct getindex call, returning the value.

@rafaqz
Copy link
Owner

rafaqz commented Jul 6, 2023

It looks like a DiskArrays.jl bug with zero dimensional arrays.

Its calling Array on a SubDiskArray, so if it gets that far its not a Rasters.jl problem.

@maxfreu
Copy link
Contributor Author

maxfreu commented Jul 6, 2023

Ah ok, good to know. I thought it's maybe an error with how DiskArray is called - but I must admit that I'm not familiar with these internals. Maybe @meggart knows what's going wrong? Should I open an issue there?

@meggart
Copy link

meggart commented Jul 7, 2023

Yes, it was a DiskArrays.jl bug. The fix is here: meggart/DiskArrays.jl#100

@rafaqz
Copy link
Owner

rafaqz commented Jul 7, 2023

Thanks for the quick fix @meggart

@rafaqz rafaqz closed this as completed Jul 7, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
3 participants