Read raster into memory automatically for copy and deepcopy?

[vlandau]

This code doesn't work

using Rasters
url_base = "https://raw.githubusercontent.com/Circuitscape/datasets/main/"

download(string(url_base, "data/nlcd_2016_frederick_md.tif"),
         "local/nlcd_2016_frederick_md.tif")

landscape = Raster("local/nlcd_2016_frederick_md.tif", missingval = -9999)

the_copy = deepcopy(landscape)
the_copy[1,1,1] = 0x05

Ouput:

julia> the_copy[1,1,1] = 0x05
ERROR: LoadError: StackOverflowError:
Stacktrace:
  [1] Array
    @ ./boot.jl:452 [inlined]
  [2] Array
    @ ./boot.jl:459 [inlined]
  [3] similar
    @ ./array.jl:359 [inlined]
  [4] writeblock!(::DiskArrays.BroadcastDiskArray{Float64, 3, Base.Broadcast.Broadcasted{DiskArrays.ChunkStyle{3}, Nothing, typeof(convert), Tuple{Base.RefValue{Type{Float64}}, Rasters.FileArray{Rasters.GDALfile, UInt8, 3, Nothing, DiskArrays.GridChunks{3}, DiskArrays.Unchunked}}}}, ::Array{Int64, 3}, ::UnitRange{Int64}, ::UnitRange{Int64}, ::UnitRange{Int64})
    @ DiskArrays ~/.julia/packages/DiskArrays/NTqYs/src/DiskArrays.jl:46
  [5] writeblock!(::DiskArrays.BroadcastDiskArray{Float64, 3, Base.Broadcast.Broadcasted{DiskArrays.ChunkStyle{3}, Nothing, typeof(convert), Tuple{Base.RefValue{Type{Float64}}, Rasters.FileArray{Rasters.GDALfile, UInt8, 3, Nothing, DiskArrays.GridChunks{3}, DiskArrays.Unchunked}}}}, ::Array{Int64, 3}, ::UnitRange{Int64}, ::UnitRange{Int64}, ::UnitRange{Int64}) (repeats 8097 times)
    @ DiskArrays ~/.julia/packages/DiskArrays/NTqYs/src/DiskArrays.jl:48
  [6] setindex_disk!(::DiskArrays.BroadcastDiskArray{Float64, 3, Base.Broadcast.Broadcasted{DiskArrays.ChunkStyle{3}, Nothing, typeof(convert), Tuple{Base.RefValue{Type{Float64}}, Rasters.FileArray{Rasters.GDALfile, UInt8, 3, Nothing, DiskArrays.GridChunks{3}, DiskArrays.Unchunked}}}}, ::Array{Int64, 3}, ::Int64, ::Vararg{Int64, N} where N)
    @ DiskArrays ~/.julia/packages/DiskArrays/NTqYs/src/DiskArrays.jl:70
  [7] setindex!
    @ ~/.julia/packages/DiskArrays/NTqYs/src/DiskArrays.jl:202 [inlined]
  [8] setindex!(::DiskArrays.BroadcastDiskArray{Float64, 3, Base.Broadcast.Broadcasted{DiskArrays.ChunkStyle{3}, Nothing, typeof(convert), Tuple{Base.RefValue{Type{Float64}}, Rasters.FileArray{Rasters.GDALfile, UInt8, 3, Nothing, DiskArrays.GridChunks{3}, DiskArrays.Unchunked}}}}, ::Int64, ::Int64, ::Int64, ::Int64)
    @ DiskArrays ~/.julia/packages/DiskArrays/NTqYs/src/DiskArrays.jl:206
  [9] setindex!(::Raster{Float64, 3, Tuple{X{Projected{Float64, LinRange{Float64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Float64}, DimensionalData.Dimensions.LookupArrays.Intervals{DimensionalData.Dimensions.LookupArrays.Start}, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALfile, Dict{Any, Any}}, WellKnownText{GeoFormatTypes.CRS, String}, Nothing, X{Colon}}}, Y{Projected{Float64, LinRange{Float64}, DimensionalData.Dimensions.LookupArrays.ReverseOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Float64}, DimensionalData.Dimensions.LookupArrays.Intervals{DimensionalData.Dimensions.LookupArrays.Start}, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALfile, Dict{Any, Any}}, WellKnownText{GeoFormatTypes.CRS, String}, Nothing, Y{Colon}}}, Band{DimensionalData.Dimensions.LookupArrays.Categorical{Int64, UnitRange{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}, Tuple{}, DiskArrays.BroadcastDiskArray{Float64, 3, Base.Broadcast.Broadcasted{DiskArrays.ChunkStyle{3}, Nothing, typeof(convert), Tuple{Base.RefValue{Type{Float64}}, Rasters.FileArray{Rasters.GDALfile, UInt8, 3, Nothing, DiskArrays.GridChunks{3}, DiskArrays.Unchunked}}}}, Symbol, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALfile, Dict{Symbol, Any}}, Int64}, ::Int64, ::Int64, ::Int64, ::Int64)
    @ Rasters ~/.julia/packages/Rasters/W22nT/src/array.jl:250

It looks like the problem is due to the Raster being entirely disk-backed as of the call of copy.

The code below does work

using Rasters

url_base = "https://raw.githubusercontent.com/Circuitscape/datasets/main/"

download(string(url_base, "data/nlcd_2016_frederick_md.tif"),
         "local/nlcd_2016_frederick_md.tif")

landscape = Raster("local/nlcd_2016_frederick_md.tif")[:,:,:] # <- Force raster to be read into memory

the_copy = deepcopy(landscape)
the_copy[1,1,1] = 0x05

To be clear, I recognize at this time this may be expected behavior, but it may be user-friendly to add new methods to copy and deepcopy that first read the raster into memory (if it isn't already), and then return a fully memory-backed copy, so the user can skip that step?

Of course, this is a major edge case, and therefore low priority. I can use the workaround above in the mean time!

Thanks for all your work to provide GIS functionality in Julia @rafaqz!

EDIT: Simplified code for clarity

[rafaqz]

Thanks, that makes sense.

So it's making a perfectly fine copy of the object, but the copy actually points at the same file. And that's kind of weird for how people think about cpoy/deepcopy. I guess copy and deepcopy should be the same here?

So maybe copy/deepcopy:

• work as usual if the file is memory backed
• if the file is disk backed (Rasters.isdisk(A) == true) then follow copy/deepcopy wtih read it to move to memory?

(btw you can use read to move to memory as you do with indexing)

We will have to go through and manually copy/deepcopy the dims/metadata etc fields if we override copy/deepcopy, and rebuild the object.

[vlandau]

Ah yes, looks like copy works just fine. I'm still slightly unsure about when to use copy vs deepcopy in general...

I think your proposal above makes sense.

[rafaqz]

I'm wondering if this needs a warning if the file is disk-based?

"Warning: this raster is disk based. If you mean to copy the file, use cp(src, dst). For this operation, raster read to memory with read.

Also, we could actually define cp:

cp(src::Raster, dst::AbstractString; kw...)` = `Raster(copy(filename(src), dst; kw...)`

And have a similar warning if it isn't disk-based.

[vlandau]

I think that makes sense -- could do a warning or even an error