Setting a fine-scale and coarse-scale modeling area

[echelleburns]

Hi Dan,

We're hoping to use marlin to model a tuna fishery. I know we'll need to model the entire stock region in order to get results that make sense, but that requires some balancing of spatial resolutions and run times.

I was wondering if it was possible to set up marlin to model at two spatial resolutions:

• a coarse resolution for the species range
• a finer resolution for the "fishing grounds"

The study region would then look something like this:

The hope would be to still model the tuna species at the broader, stock level, but then to also be able to model fleet outputs (catch, effort, etc) at the finer resolution for the area of interest with/without an MPA.

Wondering about your thoughts on this, whether it's something that marlin might be able to do in the nearish future, and/or if you have thought of any other workarounds that might be worth considering.

Thanks in advance!

[DanOvando]

Hmmm interesting question! In theory it might be possible. As part of the model you need to set a patch area, which at the moment is a constant (e.g. 100km2). In theory, patch area could be made patch specific instead, since all that happens before the matrix multiplication stuff.

In practice would take some work though. The most obvious thing is adding in that functionality to directly affected lines of code see instances in lines below

marlin/R/fish_class.R

Line 553 in 771fb16

taxis_matrix[[i]] <- pmin(exp((time_step * outer(taxis_matrix[[i]], taxis_matrix[[i]], "-")) / sqrt(patch_area)), max_hab_mult) # convert habitat gradient into diffusion multiplier

marlin/R/simmar.R

Line 47 in 771fb16

patch_area <- unique(purrr::map_dbl(fauna, "patch_area"))

After that though, would have to think through downstream effects. A few that come to mine

• Would need to adjust carrying capacity in each patch to reflect the different areas
• Would need to adjust fleet model calculations based on distance from port

It's a cool idea, but realistically not something that I'd have time to get to anytime soon. But, you or anyone else is welcome to make a branch and try to implement it, and I could advise.

Workarounds

• Set steepness to 1 and scale the size of the population to the area in question. The assumption here would be that there is no immigration or emigration of post-recruits, but that there is a subsidy of recruits that arrive from somewhere independent of the spawning biomass in the study area, which sort of gets at the idea of this being part of a bigger population, but depending on the size of the area in question raises other problems.

• Consider whether the resolution of interest is really meaningful given the population dynamics and the data. e.g. a map might look nicer if it has pixels separating out protection in the east side versus west side of a 1KM bay, but we clearly don't have the data or models to actually make reliable statements about the impact of that choice on something like a tuna population. So I guess recommendation is consider what is actually the finest resolution the data and model can provide useful insights on, and if the computational bottlenecks still crop up at that resolution?

[DanOvando]

I like the idea though, and might be able to get to it in a few months, but can't promise anything if a project is riding on this now!

[echelleburns]

Thanks, this is super useful information! This project is still in the scoping/prep stage, so if you do find time to implement this in the next couple of months, please let me know :)

If we end up trying to incorporate this before you get around to it, I'll let you know and maybe we can put some time on our schedules to figure out all the kinks.

[echelleburns]

Hi Dan,

Just following up on this thread to see if this is still something you'd be willing to advise on. I'm not sure of your capacity, but I'd love to work with you to get this feature added to marlin.

[DanOvando]

Hi Echelle,

I don't have much time to develop on this right now but would be happy to advise some on how to tackle it. I could basically provide a list of the things that would need to be done if you'd be willing to create a branch and take a stab at it

[echelleburns]

Hi Dan - that would be great if you have the capacity to make a to-do list for me. I'd love to create a branch and put some time in to see if we can get something working!

[DanOvando]

@echelleburns:

• Create a new branch to work on this
• Convert patch_area into flexible parameter that can either be a scalar or a array of some type (matrix or flat array) with dimensions matching resolution of the system
• for now, going to be on user to think about what to set as habitat in the different sized patches
• Propagate that patch area correctly through the model, so that the right patch area is affecting the right patch

  marlin/R/fish_class.R

  Line 553 in 771fb16

  taxis_matrix[[i]] <- pmin(exp((time_step * outer(taxis_matrix[[i]], taxis_matrix[[i]], "-")) / sqrt(patch_area)), max_hab_mult) # convert habitat gradient into diffusion multiplier

  marlin/R/simmar.R

  Line 47 in 771fb16

  patch_area <- unique(purrr::map_dbl(fauna, "patch_area"))

  
  

  marlin/R/sim_habitat.R

  Line 50 in 025e94c

  patch_width <- sqrt(patch_area) # the width of a patch

• adjusting r0s (the number of recruits each patch can support under unfished conditions) such that they reflect patch area in some way

  marlin/R/fish_class.R

  Line 669 in 025e94c

  r0s <- recruit_habitat %>%

• as a start, ensure that rec_form == "local_habitat" or rec_form == "global_habitat" if using this hack
• Create a simple reprex showing how this all works, focusing on once you set patch area how do the resulting distributions of recruits, numbers, and biomass play out in space

[DanOvando]

Here's a minimal reprex you can play with

library(tidyverse)

library(marlin)

resolution <- c(10,6) # resolution is in squared patches, so 20 implies a 20X20 system, i.e. 400 patches

years <- 25

patch_area <- 1000 # in square km

patches <- expand_grid(
  x = seq(1, resolution[1]),
  y = seq(1, resolution[2])
) |> 
  mutate(
    patch_id = row_number(),
    area = patch_area
  )

recruit_habitat <- sim_habitat("bigeye",kp = .001, resolution = resolution, patch_area = patch_area, output = "list")

adult_habitat <- sim_habitat("bigeye",kp = .001, resolution = resolution, patch_area = patch_area, output = "list")


fauna <-
  list(
    "bigeye" = create_critter(
      common_name = "bigeye tuna",
      adult_home_range = 10,
      recruit_home_range = 4,
      habitat = adult_habitat$critter_distributions$bigeye,
      recruit_habitat = recruit_habitat$critter_distributions$bigeye,
      density_dependence = "local_habitat",
      resolution = resolution,
      steepness = 0.6,
      ssb0 = 1000
    )
  )


fauna$bigeye$plot_movement()

patch_things <- patches

patches$r0 <- fauna$bigeye$r0s

patches$ssb0 <- fauna$bigeye$ssb0_p

patches |> 
  pivot_longer(cols = c(r0, ssb0), names_to = "variable", values_to = "value") |> 
  group_by(variable) |>
  mutate(svalue = value / max(value)) |>
  ggplot(aes(x = x, y = y, fill = svalue)) +
  geom_tile() +
  facet_wrap(~variable, scales = "free") +
  scale_fill_viridis_c() +
  theme_minimal()

[echelleburns]

Thanks, Dan! I took a first rough crack at this. Here's what I've done so far in new branch flexible-patch-area

• added conditional here to grab the value at the same position as the taxis_matrix if the patch_area is also a matrix or vector
• I was getting an error here after allowing a matrix for patch_area - I documented the error in the notes. I think the way around this is to also vectorize patch_area if it's a matrix, so I added that piece here
• Edited here so that we just get the values of patch_area for each species as a list (if patch_area is not a matrix or vector, it just returns one value per list)
• Edited here so that we pull from the patch_area list described above for the correct species
• Scaled the r0s here where if patch_area is just a single value, no scaling happens, but if patch_area is a vector or matrix, we create a scalar equal to patch_area[[i]]/max(patch_area)) and then apply that to rec_habitat. The idea here was that cells with higher areas would have higher capacity for recruits (scalar = 1) compared to cells with smaller areas (scalar < 1).
• Added a stop sequence here if people try to pass a matrix or vector to patch_area if density_dependence is something other than global_habitat or local_habitat
• Started a vignette here using the reprex you provided - assuming we'll turn this into a true vignette once things are working properly

Here are the primary outputs using a modified version of the reprex you supplied:

Using a single value for patch_area (the norm)

library(tidyverse)

library(marlin)

# Set marlin parameters
resolution <- c(20,20) # resolution is in squared patches, so 20 implies a 20X20 system, i.e. 400 patches
years <- 25
patch_area <- 1000 # in square km
 
# Produce patch dataframe
patches <- expand_grid(
  x = seq(1, resolution[1]),
  y = seq(1, resolution[2])
) |> 
  mutate(
    patch_id = row_number(),
    area = patch_area
  )

# Generate habitat matrices
recruit_habitat <- sim_habitat("bigeye",kp = .001, resolution = resolution, patch_area = patch_area, output = "list")
adult_habitat <- sim_habitat("bigeye",kp = .001, resolution = resolution, patch_area = patch_area, output = "list")

# Build fauna object
fauna <-
  list(
    "bigeye" = create_critter(
      common_name = "bigeye tuna",
      adult_home_range = 10,
      recruit_home_range = 4,
      habitat = adult_habitat$critter_distributions$bigeye,
      recruit_habitat = recruit_habitat$critter_distributions$bigeye,
      density_dependence = "local_habitat",
      resolution = resolution,
      patch_area = patch_area,
      steepness = 0.6,
      ssb0 = 1000
    )
  )

# Plot movement for fauna 
fauna$bigeye$plot_movement()

patch_things <- patches

patches$r0 <- fauna$bigeye$r0s

patches$ssb0 <- fauna$bigeye$ssb0_p

patches |> 
  pivot_longer(cols = c(r0, ssb0), names_to = "variable", values_to = "value") |> 
  group_by(variable) |>
  mutate(svalue = value / max(value)) |>
  ggplot(aes(x = x, y = y, fill = svalue)) +
  geom_tile() +
  facet_wrap(~variable, scales = "free") +
  scale_fill_viridis_c() +
  theme_minimal()

Using a matrix for patch_area

Assuming the new patch areas look like this:

# Set marlin parameters
resolution <- c(20,20) # resolution is in squared patches, so 20 implies a 20X20 system, i.e. 400 patches
years <- 25
patch_area <- matrix(c(rep(50*50, 20*5), rep(c(rep(50*50, 5), rep(10*10, 10), rep(50*50, 5)), 10),
                  rep(50*50, 20*5)), ncol = resolution[1]) # in square km

# Produce patch dataframe
patches <- expand_grid(
  x = seq(1, resolution[1]),
  y = seq(1, resolution[2])
) |>
  mutate(
    patch_id = row_number()) |> 
  left_join(patch_area %>% 
      as.data.frame() %>%
      dplyr::mutate(y = 1:nrow(.)) %>%
      tidyr::pivot_longer(
        -y,
        names_to = "x",
        values_to = "rec_habitat",
        names_prefix = "V",
        names_transform = list(x = as.integer)
        ) 
  )

# Generate habitat matrices
recruit_habitat <- sim_habitat("bigeye",kp = .001, resolution = resolution, patch_area = patch_area, output = "list")
adult_habitat <- sim_habitat("bigeye",kp = .001, resolution = resolution, patch_area = patch_area, output = "list")

# Build fauna object
fauna <-
  list(
    "bigeye" = create_critter(
      common_name = "bigeye tuna",
      adult_home_range = 10,
      recruit_home_range = 4,
      habitat = adult_habitat$critter_distributions$bigeye,
      recruit_habitat = recruit_habitat$critter_distributions$bigeye,
      density_dependence = "local_habitat",
      resolution = resolution,
      patch_area = patch_area,
      steepness = 0.6,
      ssb0 = 1000
    )
  )

# Plot movement for fauna
fauna$bigeye$plot_movement()

# Calculate other patch metrics
patch_things <- patches # save the original output
patches$r0 <- fauna$bigeye$r0s # add the number of recruits each patch can support under unfished conditions
patches$ssb0 <- fauna$bigeye$ssb0_p # add spawning stock biomass by patch

# Plot the recruits/ssb0 starting conditions over space
patches |>
  pivot_longer(cols = c(r0, ssb0), names_to = "variable", values_to = "value") |>
  group_by(variable) |>
  mutate(svalue = value / max(value)) |>
  ggplot(aes(x = x, y = y, fill = svalue)) +
  geom_tile() +
  facet_wrap(~variable, scales = "free") +
  scale_fill_viridis_c() +
  theme_minimal()

The movement plot is obviously wonky, so something needs to be fixed there, but I'm thinking that maybe the recruits and ssb0 looks promising - lower values in cells with lower areas? I'm also wondering if I'm thinking about the patch_area matrix input the right way... it feels a little backwards to me, but that might just be some brain fatigue happening.

emLab is going offline next week, so I just wanted to give you this update before the break. I'll be back the week of new years so happy to chat then if that seems worthwhile (or we can iterate here). Thanks!