Optimize initial hex grid generation

[MathewBiddle]

When trying to build the hex grid for higher resolutions (e.g. 6) the grid generation process takes a long time.

In some cases, a person doesn't necessarily need to create the hexagons for the entire globe. Would it be possible to provide obisindicators::make_hex_res() your own bounding box to narrow down the region. For example,

obisindicators::make_hex_res(RES, xmin = -180, xmax = -40, ymin = 0, ymax = 75)

See this code:

obisindicators/R/h3.R

Lines 35 to 38 in 9b1ef2f

	hex_ids <- c(
	h3::polyfill(east, res=hex_res),
	h3::polyfill(west, res=hex_res)
	)

[mimidiorio]

Here is the H3 resolution 6 grid for US waters only.

[MathewBiddle]

I adjusted make_hex_res to use a smaller bounding box:

function (hex_res = 2) 
{
  CRS <- sf::st_crs(4326)
  east <- sf::st_sf(geom = sf::st_as_sfc(sf::st_bbox(c(xmin = 160, 
    xmax = 180, ymin = 0, ymax = 74.7), crs = CRS)))
  west <- sf::st_sf(geom = sf::st_as_sfc(sf::st_bbox(c(xmin = -180, 
    xmax = -40, ymin = 0, ymax = 74.7), crs = CRS)))
  hex_ids <- c(h3::polyfill(east, res = hex_res), h3::polyfill(west, 
    res = hex_res))
  dl_offset <- 60
  hex_sf <- purrr::map_df(hex_ids, h3::h3_to_geo_boundary_sf) %>% 
    sf::st_wrap_dateline(c("WRAPDATELINE=YES", glue::glue("DATELINEOFFSET={dl_offset}"))) %>% 
    dplyr::mutate(hexid = hex_ids)
  return(hex_sf)
}

Still took a long time, but it finished.

[MathewBiddle]

I was under the assumption that you had to generate the full grid before you could do any indicator calculations. That was an incorrect assumption.

I re-wrote a little function to generate h3 gridded indicators at a resolution of your choosing for only the data you provide. This drastically reduces the amount of memory and time to create the spatial feature as it only generates geometries for applicable hexagon identifiers.

h3_indicators <- function(occ, resolution = 9) {
  ## Compute the indicator on the grid resolution of interest for the points in the identified polygons
  # return h3 cell index for occurrences in polygon
  occ_h3 <- occ %>%
    mutate(cell = h3::geo_to_h3(data.frame(decimalLatitude, decimalLongitude), res = RES))
  
  # group by cell index and compute indicators
  idx <- obisindicators::calc_indicators(occ_h3)
  
  # convert hexagon ids to spatial features
  # NOTE: DATELINEOFFSET is inv proportional to hex_res b/c we need to look
  #       further from the dateline as hex sizes get bigger.
  dl_offset <- 60  # 60 is enough for hex_res >= 1. res 0 is weird; don't use it.
  hex_sf <- purrr::map_df(idx$cell, h3::h3_to_geo_boundary_sf) %>%
    sf::st_wrap_dateline(c(
      "WRAPDATELINE=YES",
      glue::glue("DATELINEOFFSET={dl_offset}")
    )) %>%
    dplyr::mutate(hexid = idx$cell)
  
  # merge geometry into indicator table
  grid <- hex_sf %>%
    inner_join(idx, by = c("hexid" = "cell"))
  
  return(grid)

If this seems like a decent solution that should be added to the package, I can make the change to

obisindicators/R/h3.R

Lines 23 to 38 in 9b1ef2f

	east <- sf::st_sf(
	geom = sf::st_as_sfc(
	sf::st_bbox(
	c(xmin = 0, xmax = 180, ymin = -90, ymax = 90),
	crs = CRS)))
	west <- sf::st_sf(
	geom = sf::st_as_sfc(
	sf::st_bbox(
	c(xmin = -180, xmax = 0, ymin = -90, ymax = 90),
	crs = CRS)))
	# mapview(world)
	hex_ids <- c(
	h3::polyfill(east, res=hex_res),
	h3::polyfill(west, res=hex_res)
	)

.

It would probably have to be some if statement that says, if your dataframe has a column of hexids, use those to generate the spatial feature. Otherwise, build the global grid and go from there.

[MathewBiddle]

See https://github.com/NOAA-GIS4Ocean/H3_indicators/blob/main/GIS4Ocean_h3_indicators.Rmd