N-dimensional interpolation with missing and nothing values

`missing` values propagate in the interpolation and the 1D trick suggestion in the documentation of filtering `missing` values out does not generalize to multiple dimensions. Unless I missed something? This begs the question: Is it possible to interpolate data with `missing` values in *n* dimensions with *n*>1? Here is a minimal example to illustrate the kind of thing I am thinking about:
```julia
julia> nx, ny = 3, 4

julia> xs, ys = 1:nx, 1:ny

julia> data = [cos(x) + sin(y) + y for x in xs, y in ys] # some data
3×4 Array{Float64,2}:
2.38177 3.4496 3.68142 3.7835
1.42532 2.49315 2.72497 2.82705
0.851478 1.9193 2.15113 2.25321

julia> data = convert(Array{Union{Float64, Missing}}, data) ;

julia> data[2:end-1, 2:end-1] .= missing ; data # with a hole of missing values
3×4 Array{Union{Missing, Float64},2}:
2.38177 3.4496 3.68142 3.7835
1.42532 missing missing 2.82705
0.851478 1.9193 2.15113 2.25321
```
With the `data` above, I'm not even sure how one should interpolate values *near* the `missing` entries. Currently, as mentioned at the top of this post, the `missing` values propagate:
```julia
julia> knots = (xs, ys) ;

julia> itp = interpolate(knots, data, Gridded(Linear()))
3×4 interpolate((1:3,1:4), ::Array{Union{Missing, Float64},2}, Gridded(Linear())) with element type Float64:
missing missing missing missing
missing missing missing missing
missing missing missing missing
```
Is there a good way to interpolate data close to the `missing` entries? Maybe by applying the filtering trick in each direction as a default? (Of course, one could use the `Constant()` interpolation instead of the `Linear()` one, but the problem of interpolation *near* the original `missing` values remains.)

I would also be interested in having interpolation *boundaries* in the sense that the interpolation should not go across those boundaries. For that, I thought using `nothing` values made sense? Here is an example of what I am thinking about, adding a `nothing` in place of one of the missing entries:
```julia
julia> data = convert(Array{Union{Float64, Missing, Nothing}}, data) ;

julia> data[2,2] = nothing ; data # Now with a nothing in there
3×4 Array{Union{Missing, Nothing, Float64},2}:
2.38177 3.4496 3.68142 3.7835
1.42532 nothing missing 2.82705
0.851478 1.9193 2.15113 2.25321
```
In this case, I feel like querying for values near the remaining `missing` should not use the `[2,1]` entry (the `1.42532`) because it is separated from the query by a `nothing` value.

Maybe these features are dumb to request, I'm not sure 😅, but I think they could be useful for some applications, like geographic data interpolation. For example, consider a marine tracer for which a `missing` value at location `(x,y)` could mean "no data was collected at `(x,y)`" while `nothing` could represent "`(x,y)` location is an island". In this case it would make sense to be able to interpolate the data in 2D across holes of `missing` data, while *blocking* the interpolation through islands where, in terms of the marine tracer data, `nothing` is observed.