# Interpolate over large HDF5 arrays

**URL:** <https://discourse.julialang.org/t/interpolate-over-large-hdf5-arrays/127079>\
**Category:** General Usage\
**Tags:** hdf5, interpolations\
**Created:** [March 18, 2025, 4:17am UTC](https://discourse.julialang.org/t/interpolate-over-large-hdf5-arrays/127079 "2025-03-18T04:17:06Z")\
**Posts on this page:** 19\
**Page:** 1

<div class="post-metadata">

**Author:** ![henry2004y](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/henry2004y/32/9284_2.png) [@henry2004y](https://discourse.julialang.org/u/henry2004y)\
**Post date:** [March 18, 2025, 4:17am UTC](https://discourse.julialang.org/t/interpolate-over-large-hdf5-arrays/127079/1 "2025-03-18T04:17:06Z")

</div>

Hi,

I have a extremely large standard HDF5 file (2 TB) which saves a 3D array, let’s say `B`. With HDF5.jl, I can access individual elements of the array in the same way as this thread shows: [load-hdf5-file-larger-than-memory](https://discourse.julialang.org/t/load-hdf5-file-larger-than-memory/107412/1).

```julia
using HDF5

fid = h5open(fname, "r")
B = fid["B"] # size(B) = (10000, 10000, 4992), eltype(B) = Float32
B[10,10, 1] # this works

```

However, when I tried to perform interpolations over `B` using Interpolations.jl, I got

```julia
using Interpolations

itp = extrapolate(interpolate(B, BSpline(Linear())), NaN)

```

```julia
ERROR: MethodError: no method matching interpolate(::HDF5.Dataset, ::BSpline{Linear{Throw{OnGrid}}})
The function `interpolate` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  interpolate(::AbstractArray, ::IT, ::GT) where {IT<:Union{NoInterp, Tuple{Vararg{Union{NoInterp, BSpline}}}, BSpline}, GT<:Union{NoInterp, Tuple{Vararg{Union{NoInterp, Interpolations.GridType}}}, Interpolations.GridType}}
   @ Interpolations ~/.julia/packages/Interpolations/91PhN/src/deprecations.jl:19
  interpolate(::AbstractArray, ::IT) where IT<:Union{NoInterp, Tuple{Vararg{Union{NoInterp, BSpline}}}, BSpline}
   @ Interpolations ~/.julia/packages/Interpolations/91PhN/src/b-splines/b-splines.jl:190
  interpolate(::AbstractArray, ::IT, ::Real, ::Int64) where IT<:Union{NoInterp, Tuple{Vararg{Union{NoInterp, BSpline}}}, BSpline}
   @ Interpolations ~/.julia/packages/Interpolations/91PhN/src/b-splines/b-splines.jl:217

```

This means that Interpolations.jl expects an AbstractArray as input, but HDF5.Dataset is not an AbstractArray. This array is so large that I cannot possibly read it into memory as a whole. Is it possible to interpolate over such large array?

---

<div class="post-metadata">

**Author:** ![Salmon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/salmon/32/22968_2.png) [@Salmon](https://discourse.julialang.org/u/Salmon)\
**Post date:** [March 18, 2025, 5:54am UTC](https://discourse.julialang.org/t/interpolate-over-large-hdf5-arrays/127079/2 "2025-03-18T05:54:24Z")

</div>

Sounds like [HDF5 memory mapping](https://juliaio.github.io/HDF5.jl/stable/#Memory-mapping) would be good here.

---

<div class="post-metadata">

**Author:** ![henry2004y](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/henry2004y/32/9284_2.png) [@henry2004y](https://discourse.julialang.org/u/henry2004y)\
**Post date:** [March 18, 2025, 12:41pm UTC](https://discourse.julialang.org/t/interpolate-over-large-hdf5-arrays/127079/3 "2025-03-18T12:41:40Z")

</div>

Yes, but unfortunately when I checked with `HDF5.ismmappable`, it returned false. I don’t know why, but probably it’s compressed.

---

<div class="post-metadata">

**Author:** ![Salmon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/salmon/32/22968_2.png) [@Salmon](https://discourse.julialang.org/u/Salmon)\
**Post date:** [March 18, 2025, 1:04pm UTC](https://discourse.julialang.org/t/interpolate-over-large-hdf5-arrays/127079/4 "2025-03-18T13:04:09Z")

</div>

I see, it would make sense to compress such a big data set of course.

What should probably still work is to simply create a wrapper that subtypes `AbstractArray{T,N}`.  
The following seems to work on first glance, though I am not sure if Interpolations does not allocate the full array under the hood somehow. The reason I am suspecting something like that is because I can still access the interpolation after deleting the file.

```julia
using HDF5, Interpolations

struct HDF5MemoryArray{T,N} <: AbstractArray{T,N}
    data::HDF5.Dataset
    function HDF5MemoryArray(data::HDF5.Dataset)
        T = eltype(data)
        N = ndims(data)
        new{T,N}(data)
    end
end
Base.size(a::HDF5MemoryArray) = size(a.data)
Base.getindex(a::HDF5MemoryArray, i::Vararg{Int,N}) where N = a.data[i...]
Base.getindex(a::HDF5MemoryArray, i::Int) = a.data[i]

itp = h5open("a.h5","r") do file
        
    B = HDF5MemoryArray(file["data"])
        
    itp = extrapolate(interpolate(B, BSpline(Linear())), NaN)
end

```

But maybe this is a starting point

Edit:looks like the information is stored in itp.itp.coefs. I think this is more likely an issue with interpolations not being able to use a lazy interpolation here…

---

<div class="post-metadata">

**Author:** ![henry2004y](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/henry2004y/32/9284_2.png) [@henry2004y](https://discourse.julialang.org/u/henry2004y)\
**Post date:** [March 18, 2025, 1:09pm UTC](https://discourse.julialang.org/t/interpolate-over-large-hdf5-arrays/127079/5 "2025-03-18T13:09:16Z")

</div>

I’ll try! Meanwhile, I found this issue in the HDF5.jl repository: [Dataset as an AbstractArray](https://github.com/JuliaIO/HDF5.jl/issues/930#issue-1225987340). What is stopping this from being implemented?

---

<div class="post-metadata">

**Author:** ![Salmon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/salmon/32/22968_2.png) [@Salmon](https://discourse.julialang.org/u/Salmon)\
**Post date:** [March 18, 2025, 1:14pm UTC](https://discourse.julialang.org/t/interpolate-over-large-hdf5-arrays/127079/6 "2025-03-18T13:14:24Z")

</div>

At first glance at the github issue it looks like its more about technicalities, i.e. avoiding breaking changes. As far as i can tell there should be no issue in creating your own wrapper until `HDF5.Dataset` is an `AbstractArray`

---

<div class="post-metadata">

**Author:** ![henry2004y](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/henry2004y/32/9284_2.png) [@henry2004y](https://discourse.julialang.org/u/henry2004y)\
**Post date:** [March 18, 2025, 2:04pm UTC](https://discourse.julialang.org/t/interpolate-over-large-hdf5-arrays/127079/7 "2025-03-18T14:04:19Z")

</div>

I quickly tried this on the large HDF5 array. As you mentioned, this might be a problem now for Interpolations.jl:

```julia
julia> itp = extrapolate(interpolate(B, BSpline(Linear())), NaN)
ERROR: OutOfMemoryError()
Stacktrace:
  [1] GenericMemory
    @ ./boot.jl:516 [inlined]
  [2] new_as_memoryref
    @ ./boot.jl:535 [inlined]
  [3] Array
    @ ./boot.jl:585 [inlined]
  [4] Array
    @ ./boot.jl:593 [inlined]
  [5] Array
    @ ./boot.jl:599 [inlined]
  [6] padded_similar(::Type{Float32}, inds::Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}, Base.OneTo{Int64}})
    @ Interpolations ~/.julia/packages/Interpolations/91PhN/src/b-splines/prefiltering.jl:5
  [7] copy_with_padding(::Type{Float32}, A::HDF5MemoryArray{Float32, 3}, it::BSpline{Linear{Throw{OnGrid}}})
    @ Interpolations ~/.julia/packages/Interpolations/91PhN/src/b-splines/prefiltering.jl:23
  [8] prefilter(::Type{Float32}, ::Type{Float32}, A::HDF5MemoryArray{Float32, 3}, it::BSpline{Linear{Throw{OnGrid}}})
    @ Interpolations ~/.julia/packages/Interpolations/91PhN/src/b-splines/prefiltering.jl:39
  [9] interpolate(::Type{Float32}, ::Type{Float32}, A::HDF5MemoryArray{Float32, 3}, it::BSpline{Linear{Throw{OnGrid}}})
    @ Interpolations ~/.julia/packages/Interpolations/91PhN/src/b-splines/b-splines.jl:167
 [10] interpolate(A::HDF5MemoryArray{Float32, 3}, it::BSpline{Linear{Throw{OnGrid}}})
    @ Interpolations ~/.julia/packages/Interpolations/91PhN/src/b-splines/b-splines.jl:191
 [11] top-level scope
    @ REPL[13]:1

```

Maybe I should try some other interpolation libraries?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 18, 2025, 4:50pm UTC](https://discourse.julialang.org/t/interpolate-over-large-hdf5-arrays/127079/8 "2025-03-18T16:50:28Z")

</div>

> [@henry2004y](#):
>
> Maybe I should try some other interpolation libraries?

Or you could implement it yourself. Linear interpolation on a regular grid is just a few lines of code.

---

<div class="post-metadata">

**Author:** ![henry2004y](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/henry2004y/32/9284_2.png) [@henry2004y](https://discourse.julialang.org/u/henry2004y)\
**Post date:** [March 18, 2025, 6:49pm UTC](https://discourse.julialang.org/t/interpolate-over-large-hdf5-arrays/127079/9 "2025-03-18T18:49:14Z")

</div>

Following the suggestions, I quickly came up the following code for handling linear interpolations:

```julia
using HDF5

struct HDF5MemoryArray{T,N} <: AbstractArray{T,N}
   data::HDF5.Dataset
   function HDF5MemoryArray(data::HDF5.Dataset)
      T = eltype(data)
      N = ndims(data)
      new{T,N}(data)
   end
end
Base.size(a::HDF5MemoryArray) = size(a.data)
Base.getindex(a::HDF5MemoryArray{T,N}, i::Vararg{Int,N}) where {T, N} = a.data[i...]::T
Base.getindex(a::HDF5MemoryArray{T,N}, i::Int) where {T, N} = a.data[i]::T

"""
    trilin_reg(x, y, z, Q)

Trilinear interpolation for x1,y1,z1=(0,0,0) and x2,y2,z2=(1,1,1)
Q's are surrounding points such that Q000 = F[0,0,0], Q100 = F[1,0,0], etc.
"""
function trilin_reg(x::T, y::T, z::T, Q000, Q100, Q010, Q110, Q001, Q101, Q011, Q111) where {T<:Real}
   oneT = one(T)
   mx = oneT - x
   my = oneT - y
   mz = oneT - z
   fout =
      Q000 * mx * my * mz +
      Q100 * x * my * mz +
      Q010 * y * mx * mz +
      Q110 * x * y * mz +
      Q001 * mx * my * z +
      Q101 * x * my * z +
      Q011 * y * mx * z +
      Q111 * x * y * z
end

"""
    grid_interp(x, y, z, ix, iy, iz, field)

Interpolate a value at (x,y,z) in a field. `ix`,`iy` and `iz` are indexes for x, y and z
locations (0-based).
"""
grid_interp(x::T, y::T, z::T, ix::Int, iy::Int, iz::Int, field::AbstractArray{U,3}) where
   {T<:Real, U<:Number} =
   trilin_reg(x-ix, y-iy, z-iz,
      field[ix+1, iy+1, iz+1],
      field[ix+2, iy+1, iz+1],
      field[ix+1, iy+2, iz+1],
      field[ix+2, iy+2, iz+1],
      field[ix+1, iy+1, iz+2],
      field[ix+2, iy+1, iz+2],
      field[ix+1, iy+2, iz+2],
      field[ix+2, iy+2, iz+2]
   )

fname = "data.h5"

fid = h5open(fname, "r")

B = HDF5MemoryArray(fid["B"])

xGrid = range(-0.5, 0.5, length=size(B,1))
yGrid = range(-0.5, 0.5, length=size(B,2))
zGrid = range(-1, 1, length=size(B,3))

dx = xGrid[2] - xGrid[1]
dy = yGrid[2] - yGrid[1]
dz = zGrid[2] - zGrid[1]

locx = 0.0
locy = 0.0
locz = 0.0

x = (locx - xGrid[1]) / dx
y = (locy - yGrid[1]) / dy
z = (locz - zGrid[1]) / dz

# Find surrounding points (0-based indices)
ix = floor(Int, x)
iy = floor(Int, y)
iz = floor(Int, z)

grid_interp(x, y, z, ix, iy, iz, B)

```

Additional type declarations for `Base.getindex` are required for type stability.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [March 19, 2025, 6:12am UTC](https://discourse.julialang.org/t/interpolate-over-large-hdf5-arrays/127079/10 "2025-03-19T06:12:21Z")

</div>

@henry2004y do you have an example HD5F file that you could share?

---

<div class="post-metadata">

**Author:** ![henry2004y](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/henry2004y/32/9284_2.png) [@henry2004y](https://discourse.julialang.org/u/henry2004y)\
**Post date:** [March 19, 2025, 2:10pm UTC](https://discourse.julialang.org/t/interpolate-over-large-hdf5-arrays/127079/11 "2025-03-19T14:10:21Z")

</div>

Maybe a newly generated small HDF5 file can be used as an example?

```julia
using HDF5

h5open("example.h5", "w") do fid
   g = create_group(fid, "mygroup")
   dset = create_dataset(g, "myvector", Float32, (10,10,10))
   values = [Float32(i) for i in 1:10, j in 1:10, k in 1:10]
   write(dset, values)
end

fid = h5open("example.h5", "r")
g = fid["mygroup"]
# Assume the new struct and linear interpolation codes are loaded
B = g["myvector"] |> HDF5MemoryArray

xGrid = range(-0.5, 0.5, length=10)
yGrid = range(-0.5, 0.5, length=10)
zGrid = range(-1, 1, length=10)

dx = xGrid[2] - xGrid[1]
dy = yGrid[2] - yGrid[1]
dz = zGrid[2] - zGrid[1]

startx = 0.0
starty = 0.0
startz = 0.0

x = (startx - xGrid[1]) / dx
y = (starty - yGrid[1]) / dy
z = (startz - zGrid[1]) / dz

# Find surrounding points (0-based indices)
ix = floor(Int, x)
iy = floor(Int, y)
iz = floor(Int, z)

grid_interp(x, y, z, ix, iy, iz, bx) # returns 5.5

```

---

<div class="post-metadata">

**Author:** ![fabiangans](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fabiangans/32/2624_2.png) [@fabiangans](https://discourse.julialang.org/u/fabiangans)\
**Post date:** [March 19, 2025, 4:08pm UTC](https://discourse.julialang.org/t/interpolate-over-large-hdf5-arrays/127079/12 "2025-03-19T16:08:14Z")

</div>

Just for reference, depending on the use case, this might not be the most performant solution. An alternative AbstractArray implementation could use the [DiskArrays.jl](https://github.com/JuliaIO/DiskArrays.jl) interface which will always try to access the array in large blocks instead of relying on random access to single points as the pure AbstractArray interface. The following example would create an interpolated HDF5-based array.

### Create a small wrapper around HDF5 Dataset that implements the DiskArrays interface

```julia
using HDF5, DiskArrays
import DiskArrays: eachchunk, haschunks, readblock!, writeblock!, GridChunks, Chunked, Unchunked

#Implement the DiskArray interface for HDF5 as proposed in https://github.com/JuliaIO/HDF5.jl/issues/615
struct HDF5DiskArray{T,N,CS} <: AbstractDiskArray{T,N}
  ds::HDF5.Dataset
  cs::CS
end
Base.size(x::HDF5DiskArray) = size(x.ds)
haschunks(x::HDF5DiskArray{<:Any,<:Any,Nothing}) = Unchunked()
haschunks(x::HDF5DiskArray) = Chunked()
eachchunk(x::HDF5DiskArray{<:Any,<:Any,<:GridChunks}) = x.cs
readblock!(x::HDF5DiskArray, aout, r::AbstractUnitRange...) = aout .= x.ds[r...]
writeblock!(x::HDF5DiskArray, v, r::AbstractUnitRange...) = x.ds[r...] = v
function HDF5DiskArray(ds::HDF5.Dataset)
    cs = try
        GridChunks(ds, HDF5.get_chunk(ds))
    catch
        nothing
    end
    HDF5DiskArray{eltype(ds),ndims(ds),typeof(cs)}(ds,cs)
end

```

### Create a Test dataset and wrap it

```julia
#Create example file
h5open("example.h5", "w") do fid
   g = create_group(fid, "mygroup")
   dset = create_dataset(g, "myvector", Float32, (10,10,10),chunk=(5,5,5))
   values = [Float32(i) for i in 1:10, j in 1:10, k in 1:10]
   write(dset, values)
end

fid = h5open("example.h5", "r")
da = HDF5DiskArray(fid["mygroup/myvector"])

```

### Use InterpolatedDiskArray from DiskArrayTools to create an interpolated view into the array

```julia
using DiskArrayTools: InterpolatedDiskArray
using Interpolations
target_indices = (range(1,10,91),range(1,10,91),1:10)
newsize = length.(target_indices)
newchunks = GridChunks(newsize,(50,50,5))

di = InterpolatedDiskArray(da, newchunks, target_indices...,order=Quadratic())
#Compute sum over the last dimension, every chunk will only be accessed once, so this is quite efficient and works on arbitrary array sizes without running OOM
sum(di,dims=3)

```

If your workflow depends on random access to single locations into your dataset, this solution will bring no benefit. However, if you rather do some larger batch-processing processing chunks of data piece-by-piece it might be worth to look into the DiskArrays machinery.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [March 19, 2025, 4:17pm UTC](https://discourse.julialang.org/t/interpolate-over-large-hdf5-arrays/127079/13 "2025-03-19T16:17:01Z")

</div>

Given any `AbstractArray` from the above comments, you can lazily `georef` it for interpolation with GeoStats.jl:

```julia
using GeoStats

data = georef((; B)) # B from the example above

grid = CartesianGrid(10, 10, 10) # domain of interpolation

data |> InterpolateNeighbors(grid, model = NN())

```

That will give you access to various interpolation models and various types of interpolation domain (point sets, grids, meshes). The result of `InterpolateNeighbors` is currently composed of non-lazy columns, but it should be easy to work on a hotfix to reproduce the exact input array type.

---

<div class="post-metadata">

**Author:** ![henry2004y](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/henry2004y/32/9284_2.png) [@henry2004y](https://discourse.julialang.org/u/henry2004y)\
**Post date:** [March 23, 2025, 1:24pm UTC](https://discourse.julialang.org/t/interpolate-over-large-hdf5-arrays/127079/14 "2025-03-23T13:24:43Z")

</div>

Thanks for the detailed example! This indeed requires some understanding of the internals of DiskArrays.jl. I will compare the performance with the first naive example next.

When I searched for related packages, I noticed a higher level YAXArrays.jl that’s built upon DimensionalData.jl and DiskArrays.jl. Maybe this task can be directly accomplished by taking advantage of YAXArrays.jl?

---

<div class="post-metadata">

**Author:** ![henry2004y](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/henry2004y/32/9284_2.png) [@henry2004y](https://discourse.julialang.org/u/henry2004y)\
**Post date:** [November 29, 2025, 10:52pm UTC](https://discourse.julialang.org/t/interpolate-over-large-hdf5-arrays/127079/15 "2025-11-29T22:52:51Z")

</div>

Now I get another question about the combination of DiskArrays.jl with Interpolations.jl. This example shows how to get an interpolated array directly. For my use case, I need an efficient interpolation function. How would I instead get an interpolation function that accepts a location and output an interpolated value?

I guess the wrapping with Geostats as shown above is one possible solution.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [November 29, 2025, 11:04pm UTC](https://discourse.julialang.org/t/interpolate-over-large-hdf5-arrays/127079/16 "2025-11-29T23:04:06Z")

</div>

> [@henry2004y](#):
>
> I guess the wrapping with Geostats as shown above is one possible solution.

Happy to help if you need assistance to understand the code base.

---

<div class="post-metadata">

**Author:** ![henry2004y](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/henry2004y/32/9284_2.png) [@henry2004y](https://discourse.julialang.org/u/henry2004y)\
**Post date:** [November 29, 2025, 11:12pm UTC](https://discourse.julialang.org/t/interpolate-over-large-hdf5-arrays/127079/17 "2025-11-29T23:12:17Z")

</div>

Thanks @juliohm! You are so helpful as always~ I will later post an example with an ugly wrapper I have.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [November 29, 2025, 11:16pm UTC](https://discourse.julialang.org/t/interpolate-over-large-hdf5-arrays/127079/18 "2025-11-29T23:16:55Z")

</div>

No problem 🙂

We could evolve the MWE offline, and come back with a final solution to avoid polluting the discussion. What do you think? Please feel free to reach me out in private here or in our Zulip channel.

---

<div class="post-metadata">

**Author:** ![henry2004y](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/henry2004y/32/9284_2.png) [@henry2004y](https://discourse.julialang.org/u/henry2004y)\
**Post date:** [December 3, 2025, 5:52pm UTC](https://discourse.julialang.org/t/interpolate-over-large-hdf5-arrays/127079/19 "2025-12-03T17:52:48Z")

</div>

Getting DiskArrays.jl and Interpolations.jl to work together is actually pretty simple:

```julia-auto
using HDF5, DiskArrays
import DiskArrays: eachchunk, haschunks, readblock!, writeblock!, GridChunks, Chunked,
                   Unchunked
using Interpolations

# ## Implement HDF5DiskArray
#
# First, we define a wrapper around HDF5 Dataset that implements the `DiskArrays` interface.

struct HDF5DiskArray{T, N, CS} <: AbstractDiskArray{T, N}
   ds::HDF5.Dataset
   cs::CS
end

Base.size(x::HDF5DiskArray) = size(x.ds)
haschunks(x::HDF5DiskArray{<:Any, <:Any, Nothing}) = Unchunked()
haschunks(x::HDF5DiskArray) = Chunked()
eachchunk(x::HDF5DiskArray{<:Any, <:Any, <:GridChunks}) = x.cs
readblock!(x::HDF5DiskArray, aout, r::AbstractUnitRange...) = aout .= x.ds[r...]
writeblock!(x::HDF5DiskArray, v, r::AbstractUnitRange...) = x.ds[r...] = v

function HDF5DiskArray(ds::HDF5.Dataset)
   chunks = try
      HDF5.get_chunk(ds)
   catch
      nothing
   end
   cs = isnothing(chunks) ? nothing : GridChunks(size(ds), chunks)
   HDF5DiskArray{eltype(ds), ndims(ds), typeof(cs)}(ds, cs)
end

# ## Create Example Data
#
# We create a dummy HDF5 file with some field data.

filename = "example_field.h5"
h5open(filename, "w") do fid
   g = create_group(fid, "mygroup")
   ## Create a dataset with chunking enabled
   dset = create_dataset(g, "myvector", Float32, (10, 10, 10), chunk = (5, 5, 5))
   values = [Float32(i + j + k) for i in 1:10, j in 1:10, k in 1:10]
   write(dset, values)
end

# ## Interpolation Function
#
# We define a function `itp` to query the field at a single location `x`.

# Open the file and wrap the dataset
h5open(filename, "r") do fid
   da = HDF5DiskArray(fid["mygroup/myvector"])

   cached = DiskArrays.cache(da)
   itp = extrapolate(
      interpolate(cached, BSpline(Linear(Periodic(OnCell())))), Periodic(OnCell()))

   # Evaluate at a point
   loc_int = (5.0, 5.0, 5.0)
   println("Value at $loc_int: ", itp(loc_int...))

   loc_float = (5.5, 5.5, 5.5)
   println("Value at $loc_float: ", itp(loc_float...))

   loc_out = (-0.5, 1.0, 1.0)
   val_periodic = itp(loc_out...)
   println("Value at $loc_out (Periodic): ", val_periodic)
   ## Check correctness of Periodic
   ## Note that we assume cell center values. -0.5 wraps to 9.5 (since period is 10).
   ## Value at 9.5, 1.0, 1.0. Data is i+j+k.
   println("Expected Periodic: ", 9.5 + 1 + 1)
end

# ## Cleanup
# Remove the temporary file.

rm(filename)

```