0

I have a 3D array that represents a lake basin, some of the 3D array is land (-999 values) and is not of interest. Within the lake, cells can contain either pelagic resource or littoral resources. These two resources grow at different rates and can be treated as independent, non-interacting parts of the lake.

Within the littoral and within the pelagic cells, I want to average cell values with neighboring cells (Chebyshev distance = 1), but only if the neighboring cells are of the same type, e.g. pelagic cells should only average with neighboring pelagic cells, not littoral or land cells and littoral cells should only average with neighboring littoral cells, not pelagic or land cells.

I have created a small hypothetical example.

using Distributions

# create hypothetical lake bathymetry
foo(N) = [min(i, j, N+1-i, N+1-j) for i in 1:N, j in 1:N]

# 10 x 10 surface area, with max depth of 5
bathy = foo(10)

# delineate land (-999) - area not interested in
bathy[bathy .< 2] .= -999

# delineate either littoral or pelagic part of lake
lake_type = ones(Int, size(bathy))
lake_type .= bathy

# if depth is between 0 and 4 cell is littoral (1)
lake_type[lake_type .> -1 .&& lake_type .< 4] .= 1

# if depth is > than 3 cell is pelagic (0)
lake_type[lake_type .> 3] .= 0

# create 3D lake resource array --------

# max depth 
mx_dpth = maximum(bathy)

# dims of basin
dims = (10,10, mx_dpth)

# 3d array of basin - to store resources
# [:, :, 1] = deepest layer of the lake
# [:, :, 5] = surface layer of the lake
basal_resource = zeros(dims)

# populate with resource - varies depending if littoral (0) or pelagic (1)
for i in 1:dims[1], j in 1:dims[2]
    if lake_type[i,j] == 0    # initial pelagic resource amount 
        basal_resource[i, j, (mx_dpth - bathy[i,j] + 1):mx_dpth] .= round.(rand(Uniform(100, 150), 1), sigdigits=3)
    end
    if lake_type[i,j] == 1 # initial littoral  resource amount
        basal_resource[i, j, (mx_dpth - bathy[i,j] + 1):mx_dpth] .= round.(rand(Uniform(1, 5), 1), sigdigits=3)
    end
end

# convert "land" cells to -999
basal_resource[basal_resource .== 0.0] .= -999

I'm now stuck trying to average only neighboring cells (in 3 dimensions) in basal_resource of the same type. Any advice would be much appreciated, hopefully my question makes sense.

EDIT

The 2D lake_type array is only relevant at the surface layer, as cells may be water at the surface, but then become land deeper down. Here is a 3D array indicating lake type through the whole water column.

# create 3d "map" of the lake indicating if cell is pelagic (0) or littoral (1)
lake_type_3d = copy(basal_resource)
replace!(x -> x>100.0 ? 0.0 : x, lake_type_3d)
replace!(x -> x>1.0 ? 1.0 : x, lake_type_3d)
4
  • If you are okay with using a package for this, check out github.com/JuliaDynamics/Agents.jl.
    – xlxs4
    Commented Jun 12 at 19:23
  • I am familiar with the agents package, but I'm not sure how it would be used in this situation?
    – flee
    Commented Jun 12 at 21:06
  • I assume you could define a GridSpace and query for neighboring cells.
    – xlxs4
    Commented Jun 12 at 21:14
  • Thanks, you might be correct. nearby_ids iterates over agents rather than the GridSpace, but nearby_positions could work, but it excludes the current position, and I'd still need to separate littoral and pelagic cells. I suspect this could end up being quite a slow approach as nearby_positions would be gathering lots of cells that aren't needed.
    – flee
    Commented Jun 12 at 21:44

1 Answer 1

2

Assuming the borders of the lake_type variable are always land (so the average doesn't really matter), you could do this with a loop like this:

avg_same = zeros(size(basal_resource))

for i in 2:size(basal_resource, 1)-1, j in 2:size(basal_resource, 2)-1
    neighbors2 = repeat(lake_type[i-1:i+1, j-1:j+1] .== lake_type[i, j], outer = [1, 1, 2])
    neighbors3 = repeat(lake_type[i-1:i+1, j-1:j+1] .== lake_type[i, j], outer = [1, 1, 3])
    avg_same[i, j, 1] = sum(basal_resource[i-1:i+1, j-1:j+1, 1:2] .* neighbors2) / sum(neighbors2)
    avg_same[i, j, end] = sum(basal_resource[i-1:i+1, j-1:j+1, end-1:end] .* neighbors2) / sum(neighbors2)
    for k in 2:size(basal_resource, 3)-1
        avg_same[i, j, k] = sum(basal_resource[i-1:i+1, j-1:j+1, k-1:k+1] .* neighbors3) / sum(neighbors3)
    end
end

This might not be the most performant way of accomplishing the task, but for decently sized arrays a nested loop shouldn't be too taxing.

EDIT

If you use a 3D version of lake_type, you can do the same kind of matrix multiplication logic:

avg_same = zeros(size(basal_resource))

for i in 2:size(basal_resource, 1)-1, j in 2:size(basal_resource, 2)-1, k in 1:size(basal_resource, 3)
    kmin = max(k-1, 1)
    kmax = min(k+1, size(basal_resource, 3))
    neighbors = lake_type_3d[i-1:i+1, j-1:j+1, kmin:kmax] .== lake_type_3d[i, j, k]
    avg_same[i, j, k] = sum(basal_resource[i-1:i+1, j-1:j+1, kmin:kmax] .* neighbors) / sum(neighbors)
end
5
  • Thanks for answering, although it does not appear to work as I imagined. if you look at basal_resource[:, :, 1] (lowest layer in the "lake") there are 4 non-land cells in the middle., which should average with each other and the adjoining cells in the layer above. But if we look at avg_same[:, :, 1] the 4 lake cells have averaged with the land (-999) cells, resulting in negative values.
    – flee
    Commented Jun 12 at 1:18
  • I think this is because lake_type has a 4x4 grid of pelagic cells in the center, but basal_resource has only a 2x2 grid in the center. Is this by design? I assumed lake_type contained the correct definition. Commented Jun 12 at 14:15
  • Ah you are correct, that is confusing from my part. I have added an edit, which gives a 3D version of lake_type which matches basal_resource, if you are willing to have another look.
    – flee
    Commented Jun 12 at 21:37
  • 1
    Edited the answer to have the version for 3D lake_type. Note this still assumes that all the outer borders in i and j dimensions are land, but you can swap out the same logic for kmin and kmax if that assumption is wrong. Commented Jun 13 at 15:09
  • Thank you! Indeed in my real data not all outer cells are land, but applying the kmin and kmax logic to i and j does the trick. Much appreciated.
    – flee
    Commented Jun 13 at 22:13

Not the answer you're looking for? Browse other questions tagged or ask your own question.