Make Grid for Predictions

[AngeliaMiller]

Hello,

I am trying to predict species distributions based on biomass and depth based on a model fit to existing data and across a grid that is constrained to my historical survey area. I have not made a grid using a shapefile in R before, do you have any resources or code that could assist?

Currently, I am using

nefsc_grid <- strata %>%
sf::st_make_grid(cellsize = 0.5, what = "centers") %>%
st_intersection()

But compared to the qcs_grid in the package, my result is a list with points that I need to convert to UTMs and append with my covariates. Additionally, when I plot it, I get the resulting plot.

[seananderson]

Hi @AngeliaMiller, I think one problem above is that your sf::st_intersection() is missing the strata polygon. There are many ways to do this. Here's one example below with the sf package. The area calculation part after intersection is only important if you wanted to calculate something like an area-weighted abundance index.

If you had covariates, you'd still have to find the values at the centroid of each grid cell. There are many ways you could do that and it depends on the format of your environmental/depth data etc. One option would be with the raster package and raster::resample() assuming your environmental raster is at a different resolution.

Creating a grid on a polygon with sf:

library(dplyr)
library(ggplot2)
theme_set(theme_light())

# create some fake data
# you're more likely to start with a shape file of your region
# or maybe create a polygon based on depth contours
fake_polygon_dat <- structure(c(
  3, 4.4, 18.8, 26.8, 33.5, 54.8, 76.5, 83.1, 88.3,
  72.9, 55.4, 48.5, 35.3, 16.6, 7.3, 3, 14.6, 37.7, 52.2, 67.8,
  88, 86.4, 84.2, 71.3, 51.7, 35.5, 27.7, 15.9, 2.2, 4.6, 19.7,
  14.6
), dim = c(16L, 2L))
fake_polygon <- sf::st_polygon(list(fake_polygon_dat)) |>
  sf::st_sfc()

# choose a grid size in units of our polygon shape file
grid_spacing <- 2

# create grid over the bounding box of the polygon
full_grid <- sf::st_make_grid(
  fake_polygon,
  cellsize = c(grid_spacing, grid_spacing)
) %>%
  sf::st_sf()
ggplot(full_grid) + geom_sf()

# subset our grid to cells that intersect our polygon:
intersected <- sf::st_intersects(full_grid, fake_polygon)
selected_grid <- full_grid[lengths(intersected) > 0, ]

# plot it
ggplot() +
  geom_sf(data = fake_polygon) +
  geom_sf(data = selected_grid, fill = NA)

# find how much of each grid cell is within the outer polygon:
overlap <- sf::st_intersection(selected_grid, fake_polygon)
ggplot(overlap) + geom_sf()

calculated_area <- sf::st_area(overlap)

# extract coordinates into a data frame
# append cell area column if desired
coord <- selected_grid |>
  sf::st_centroid() |>
  sf::st_coordinates() |>
  as_tibble() |>
  mutate(area = calculated_area)

# plot the grid (colour is cell area)
ggplot(coord, aes(X, Y, colour = area)) +
  geom_tile(width = 2, height = 2, fill = NA) +
  scale_colour_viridis_c(direction = -1) +
  geom_point(size = 0.5) +
  coord_fixed()

dplyr::glimpse(coord)
#> Rows: 1,213
#> Columns: 3
#> $ X    <dbl> 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 16, 18, 20, 22, 24, 26, 2…
#> $ Y    <dbl> 3.2, 3.2, 3.2, 3.2, 3.2, 3.2, 3.2, 3.2, 3.2, 3.2, 5.2, 5.2, 5.2, …
#> $ area <dbl> 0.10568627, 0.58609626, 1.09946524, 1.61283422, 2.12620321, 2.639…

^{Created on 2022-12-12 with reprex v2.0.2}

[seananderson]

Here's an example with the raster package, which is perhaps simpler and also leaves you with a raster object that may be easier to merge with your environmental data using raster::resample(), raster::aggregate(), or similar.

library(ggplot2)
library(sf)

# create some fake data
fake_polygon_dat <- structure(c(
  3, 4.4, 18.8, 26.8, 33.5, 54.8, 76.5, 83.1, 88.3,
  72.9, 55.4, 48.5, 35.3, 16.6, 7.3, 3, 14.6, 37.7, 52.2, 67.8,
  88, 86.4, 84.2, 71.3, 51.7, 35.5, 27.7, 15.9, 2.2, 4.6, 19.7,
  14.6
), dim = c(16L, 2L))
fake_polygon <- sf::st_polygon(list(fake_polygon_dat)) |>
  sf::st_sfc()

# create a grid using the raster package
resolution <- 2
r <- raster::raster(as(fake_polygon, "Spatial"), resolution = resolution)
rr <- raster::rasterize(as(fake_polygon, "Spatial"), r, getCover = TRUE)

grid <- as.data.frame(raster::rasterToPoints(rr))
grid$area <- grid$layer * resolution * resolution
grid <- dplyr::filter(grid, area > 0) |> 
  dplyr::select(-layer)

ggplot(grid, aes(x, y, colour = area)) +
  geom_tile(width = 2, height = 2, fill = NA) +
  scale_colour_viridis_c(direction = -1) +
  geom_point(size = 0.5) +
  coord_fixed()

dplyr::glimpse(grid)
#> Rows: 1,201
#> Columns: 3
#> $ x    <dbl> 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 54, 56, 58, 32, 34, 3…
#> $ y    <dbl> 87, 87, 87, 87, 87, 87, 87, 87, 87, 87, 87, 87, 87, 87, 85, 85, 8…
#> $ area <dbl> 0.04, 3.60, 3.60, 3.28, 3.04, 2.76, 2.40, 2.08, 1.84, 1.56, 1.20,…

^{Created on 2022-12-12 with reprex v2.0.2}

[AngeliaMiller]

@seananderson thanks so much for a thorough answer! I will try what you have suggested and let you know if I come across any issues, but this has been super helpful, and I really appreciate it.

[AngeliaMiller]

Hi @seananderson,

I am just getting back to this but was able to successfully create a grid with this process (below overlayed on portions of the NEFSC survey strata).

I followed the rest of your vignettes to use a model fit to year and depth (biomass ~ as.factor(year) + depth) and predict on the grid replicated over the number of years; this resulted in a dataframe of over 275,000 observations (12 years with at most 110 depth data points). My issue now appears when I try to plot my estimates using the plot_map functions from your vignettes:

plot_map <- function(dat, column) {
ggplot(dat, aes(X, Y, fill = {{ column }})) +
geom_raster() +
facet_wrap(~EST_YEAR) +
coord_fixed()
}

saying that the pixels are placed at uneven horizontal intervals, to consider using geom_tile, and that NAs were introduced by coercion to integer range. However, when using geom_tile, I receive a blank plot.

Do you have any other suggestions on how to plot the predictions spatially over the grid?

Thanks!

[seananderson]

On the geom_tile issue, you likely just need to specify a width and height. I think it tries to guess, but if the grid is uneven then it can guess a tiny width/height making the cells invisible.

On the grid, I'm not sure what you mean by "depth data points" and 275,000 sounds like a lot of grid cells, although perhaps. Is it possible you could get away with a coarser grid? It also seems odd you'd end up with an uneven grid if you created it yourself.

[AngeliaMiller]

I'll look into the geom_tile issue some more.

Regarding the grid, each year there were, at most, 110 samples taken at various depths, and these depth values are what I fit to in my model. These values are also what I should have in my grid dataframe to predict on correct? I do think 275,000 is a lot, but I wasn't sure how else to get my depth attached to my grid.

grid_yrs <- sdmTMB::replicate_df(coord, "EST_YEAR", unique(sf_fall$EST_YEAR))
grid_yrs2 <- grid_yrs |>
left_join(depths, by = "EST_YEAR")

Currently my grid spacing is 0.25, besides that and the calculated area piece, I copied your code above almost verbatim. Should my spacing be higher or lower than that?

[AngeliaMiller]

Hi Sean,

I think I was able to figure it out, with a few changes:

1. Used a coarser grid based on your suggestion at 0.5 grid spacing to have a total of 90k grid cells;
   -- it does look pretty coarse (plot below), so I'm hoping to find a happy medium between 0.25 and 0.5 grid spacing. Currently 0.25 spacing and 275,000 cells overloads R for it to map them all.
2. While predictions were made on the UTM columns, I added my decimal degree coordinates back to my dataframe for mapping to geom_raster;
3. Also found that my coordinates were too precise at 4 decimal places for geom_raster to handle, and had to round them to 2 decimal places, applied it both to my UTM and long/lats (ras_df <- preds |> mutate(across(c(X, Y, longitude, latitude), round, digits = 2)))
   -- these two resources helped with the two errors/warnings I was getting:
   -- Stack Overflow: Problem with geom_raster size vector cannot be allocated
   -- geom_raster.R in ggplot2 repo

[seananderson]

Typically your grid would be made in UTMs so the area per grid cell is constant. There's some help in the 'pretty maps' vignette https://pbs-assess.github.io/sdmTMB/articles/ on then plotting it back in lat/long with a projection. Lat/Long cells are fine, but you'll have to make sure your area is calculated appropriately for every cell first.

I'm a bit confused by:

grid_yrs <- sdmTMB::replicate_df(coord, "EST_YEAR", unique(sf_fall$EST_YEAR))
grid_yrs2 <- grid_yrs |> 
left_join(depths, by = "EST_YEAR")

It looks like you're replicating your grid at all depths. Every grid cell should have a single depth associated with it. Often this is bottom depth (unless this is some 3D projection across a bunch of depths at once?).

The estimated spatial range can help inform a reasonable grid resolution. Often the survey design already has some sampling resolution built in to work with. I'd tend to start coarser (for speed) and lower it to make sure you've hit a point of diminishing returns where you see little impact on the index.

We discuss some of this here https://peerj.com/articles/12783/

[AngeliaMiller]

Okay, I did make my grid in UTMs and will look at the vignette you linked. I am now having an issue with calculating the area though; I followed the code from above, also copied here but with my respective data:

# subset our grid to cells that intersect our polygon:
intersected <- sf::st_intersects(full_grid, fall_strata) # results in list of 741 
selected_grid <- full_grid[lengths(intersected) > 0, ] # results in 222 rows

# find how much of each grid cell is within the outer polygon:
overlap <- sf::st_intersection(selected_grid, fake_polygon) # results in 357 rows

calculated_area <- sf::st_area(overlap) # results in 357 rows

But I cannot add the area back to my grid because of the mismatch in data - 222 grid cells, and 357 area values.

Re: replicating your grid at all depths - I guess this is what I am doing, unfortunately, and definitely not what I want. I'm attaching the summer flounder data I'm using to potentially help with this discussion. I fit my model using the EXPCATCHWT values to the AVGDEPTH values (EXPCATCHWT ~ AVGDEPTH _ as.factor(EST_YEAR)). I'm not sure how to attach these depth values back to my grid cells, again if the number of values differ (222 grid cells and 110 depths for each year).

I'm sorry for all the questions and am really appreciating all your help!

sf_fall_data.csv

[seananderson]

@AngeliaMiller provided me with code and data by email. I may provide a full worked example later, maybe as a vignette, but in the meantime, the two missing pieces were:

1. Converting to UTMs first with sf::st_transform() to get a regular grid in UTMs.
2. Combining the stratum polygons into a single outer polygon first with sf::st_union(). The inner stratum boundaries were contributing to extra cells from sf::st_intersection (which does an intersection with every polygon) that didn't match sf::st_intersects (which just checks if a grid cell touches any of the polygons).