analyze_spatial.Rmd

---
title: "Spatial Query Optimizing"
output: 
  html_document:
    toc: true
    number_sections: true
editor_options: 
  chunk_output_type: console
---

```{r setup, include=F}
library(DBI)
db = dbConnect(
  RPostgres::Postgres(),
  dbname   = "gis",
  host     = "marineenergy.app",
  port     = 5432,
  user     = "admin",
  password = pass
)

knitr::opts_chunk$set(echo = T, warning = F, message = F, connection = "db")
```


# CTE

```{sql CTE, class.source='fold-show'}
explain analyze

with
tt_aoi as (
	select st_buffer(st_collect (
		'SRID=4326;
		POLYGON ((-67.06819 44.99416, -67.1857 44.94707, 
				  -67.21651 44.88058, -67.15834 44.78871, 
				  -67.04385 44.81789, -66.91015 44.86279, 
				  -67.06819 44.99416))'::geometry), 
		5) as aoi),
				  
tt_seldata as (
	SELECT 
	  sitename_l, lifestage, geometry 
	FROM 
	  shp_nationwide_efh as ds
)

SELECT 
	sitename_l AS Species, string_agg(lifestage, ', ') AS Lifestage 
FROM 
	tt_seldata as ds
	inner join tt_aoi on st_intersects(ds.geometry, tt_aoi.aoi)
GROUP BY 
     sitename_l;
```

---

# Temporary Tables


## Option 1

### Drop Temp Tables
```{sql temp-table-drop1, class.source='fold-show'}
drop table if exists tt_aoi cascade;
```

```{sql temp-table-drop2, class.source='fold-show'}
drop table if exists tt_seldata cascade;
```

### Select AOI
```{sql, temp-table-aoi, class.source='fold-show'}
--drop table if exists tt_aoi cascade;
explain analyze
--drop table if exists tt_aoi cascade;
create temporary table tt_aoi as
select st_buffer(st_collect (
	'SRID=4326;
	POLYGON ((-67.06819 44.99416, -67.1857 44.94707, 
			  -67.21651 44.88058, -67.15834 44.78871, 
			  -67.04385 44.81789, -66.91015 44.86279, 
			  -67.06819 44.99416))'::geometry), 5) as aoi;
```

### Select Data
```{sql temp-table-seldata, class.source='fold-show'}
explain analyze
create temporary table tt_seldata as
SELECT 
  sitename_l, lifestage, geometry 
FROM 
  shp_nationwide_efh as ds
```

### Summarize
```{sql temp-table-summarize, class.source='fold-show'}
explain analyze
select 
  ds.sitename_l AS Species, string_agg(ds.lifestage, ', ') AS Lifestage 
from
  tt_seldata as ds
  inner join tt_aoi on st_intersects(ds.geometry, tt_aoi.aoi)
group by
  sitename_l
```

---

## Option 2 (indexing)

> Standard database indexes create a hierarchical tree based on the values of the column being indexed. Spatial indexes are a little different – they are unable to index the geometric features themselves and instead index the bounding boxes of the features.

[Reference](https://postgis.net/workshops/postgis-intro/indexing.html)

### Drop Temp Tables
```{sql temp-table-drop1-i, class.source='fold-show'}
drop table if exists tt_aoi2 cascade;
```

```{sql temp-table-drop2-i, class.source='fold-show'}
drop table if exists tt_seldata2 cascade;
```

### Select AOI
```{sql, temp-table-aoi-i, class.source='fold-show'}
--drop table if exists tt_aoi cascade;
explain analyze
--drop table if exists tt_aoi cascade;
create temporary table tt_aoi2 as
select st_buffer(st_collect (
	'SRID=4326;
	POLYGON ((-67.06819 44.99416, -67.1857 44.94707, 
			  -67.21651 44.88058, -67.15834 44.78871, 
			  -67.04385 44.81789, -66.91015 44.86279, 
			  -67.06819 44.99416))'::geometry), 5) as aoi;
```

### Select Data
```{sql temp-table-seldata-i, class.source='fold-show'}
explain analyze
create temporary table tt_seldata2 as
SELECT 
  sitename_l, lifestage, geometry 
FROM 
  shp_nationwide_efh as ds;
```

### Create index
```{sql temp-table-seldata-index, class.source='fold-show'}
create index tt_seldata_geom_idx
  on tt_seldata2
  using gist (geometry)
```


### Summarize
```{sql temp-table-summarize-i, class.source='fold-show'}
explain analyze
select 
  ds.sitename_l AS Species, string_agg(ds.lifestage, ', ') AS Lifestage 
from
  tt_seldata2 as ds
  inner join tt_aoi2 on st_intersects(ds.geometry, tt_aoi2.aoi)
group by
  sitename_l
```

---

# Subquery

```{sql subquery-combined-i, class.source='fold-show'}
explain analyze
SELECT 
        sitename_l AS Species, string_agg(lifestage, ', ') AS Lifestage 
FROM 
	(SELECT 
			sitename_l, lifestage, geometry 
	FROM 
			shp_nationwide_efh as ds
	where
		st_dwithin(ds.geometry,
		(select st_collect (
		'SRID=4326;
		POLYGON ((-67.06819 44.99416, -67.1857 44.94707, 
				  -67.21651 44.88058, -67.15834 44.78871, 
				  -67.04385 44.81789, -66.91015 44.86279, 
				  -67.06819 44.99416))'::geometry)),
		5)) as tmp_aoi
GROUP BY 
        sitename_l
```

---

# CURRENT tabulate_dataset_shp_within_aoi

```{r ds-func, eval=F, class.source='fold-show'}
tabulate_dataset_shp_within_aoi <- function(dataset_code, aoi_wkt, output = "kable"){
  # summarize shapefile dataset from area of interest
  
  # TODO: pull from db: https://docs.google.com/spreadsheets/d/1MMVqPr39R5gAyZdY2iJIkkIdYqgEBJYQeGqDk1z-RKQ/edit#gid=936111013
  # datasets_gsheet2db()
  
  # dataset_code = "cetacean-bia"; aoi_wkt = params$aoi_wkt
  # dataset_code = "efh"; aoi_wkt = "POLYGON ((-67.06819 44.99416, -67.1857 44.94707, -67.21651 44.88058, -67.15834 44.78871, -67.04385 44.81789, -66.91015 44.86279, -67.06819 44.99416))"
  message(glue("tab..._shp_within_aoi(dataset_code='{dataset_code}', aoi_wkt='{paste(aoi_wkt, collapse=';')}')"))

  # dataset_code = "cetacean-bia";
  # params <- yaml::yaml.load("
  # title: Testing
  # aoi_wkt:
  # - POLYGON ((-115.9245 32.26236, -115.9245 32.26565, -115.9206 32.26565, -115.9206
  #               32.26236, -115.9245 32.26236))
  # - POLYGON ((-121.9503 33.01519, -121.9503 35.51658, -115.8711 35.51658, -115.8711
  #             33.01519, -121.9503 33.01519))")
  # aoi_wkt <- params$aoi_wkt
  # dataset_code = "oil-gas-wells"
  # aoi_wkt      = "POLYGON ((-157.4273 55.22198, -157.4273 61.76097, -143.1428 61.76097, -143.1428 55.22198, -157.4273 55.22198))"
  # dataset_code='fed-sand-gravel-lease'
  # aoi_wkt='POLYGON ((-175.4932 15.34568, -175.4932 27.93566, -151.813 27.93566, -151.813 15.34568, -175.4932 15.34568))'
  
  # dataset_code='monuments'
  # aoi_wkt='POLYGON ((-180.0668 16.98081, -180.0668 29.87807, -153.4797 29.87807, -153.4797 16.98081, -180.0668 16.98081))'
  
  if (is.null(aoi_wkt))
    return("Please draw a Location to get a summary of the intersecting features for this dataset.")

  # if (length(aoi_wkt) > 1)
  #   return("Please draw only ONE polygon to get a summary of the intersecting features for this dataset.")
    
  ds <- tbl(con, "datasets") %>% 
    filter(code == !!dataset_code) %>% 
    replace_na(list(buffer_nm = 0)) %>% 
    collect()
  
  if (length(aoi_wkt) > 1){
    aoi_wkts <- glue("'SRID=4326;{aoi_wkt}'::geometry")
    #aoi_sql  <- glue("ST_COLLECT(\n{paste(aoi_wkts, collapse=',\n')})")
    aoi_sql  <- glue("ST_COLLECT(\n{paste(aoi_wkts, collapse=',\n')})")
    # cat(aoi_sql)
  } else {
    aoi_sql <- glue("'SRID=4326;{aoi_wkt}'")
    # aoi_sql <- glue("'SRID=4326;{aoi_wkt[2]}'")
  }

  
  dbSendQuery(con, glue("DROP TABLE IF EXISTS tmp_aoi CASCADE;"))
  if (!is.na(ds$summarize_sql)){
    sql_intersection <- glue("
      {ds$select_sql} AS ds
      WHERE ST_DWithin(Geography(ds.geometry), {aoi_sql}, {ds$buffer_nm} * 1852);")
    dbExecute(con, glue("CREATE TEMPORARY TABLE tmp_aoi AS {sql_intersection};"))
    # sql_summarize <- "SELECT sitename_l AS Species, string_agg(lifestage, ', ') AS Lifestage FROM tmp_aoi GROUP BY sitename_l"
    # sql_summarize <- "SELECT * FROM tmp_aoi"
    #cat(ds$summarize_sql)
    x_df <- dbGetQuery(con, ds$summarize_sql)
  } else {
    x_sql <- glue("
      {ds$select_sql} AS ds
      WHERE ST_DWithin(Geography(ds.geometry), {aoi_sql}, {ds$buffer_nm} * 1852);")
    x_sf <- st_read(con, query = x_sql)
    x_df <- st_drop_geometry(x_sf)
    
    if (!is.na(ds$summarize_r))
      eval(parse(text=ds$summarize_r))
  }
  
  if (output == "tibble"){
    return(x_df)
  }
  
  x_spatial <- ifelse(
    ds$buffer_nm == 0,
    glue("\n\nSpatial: within site", .trim = F),
    glue("\n\nSpatial: within {ds$buffer_nm} nautical miles of site", .trim = F))

  if (knitr::is_html_output()){
    x_caption <- HTML(markdownToHTML(
      text = glue("Source: [{ds$src_name}]({ds$src_url}){x_spatial}"),
      fragment.only = T))
    
    tbl <- x_df %>% 
      kbl(caption = x_caption) %>%
      kable_styling(
        # full_width = F, position = "left", # position = "float_right"
        bootstrap_options = c("striped", "hover", "condensed", "responsive"))
    
  } else {
    x_caption <- glue("Source: [{ds$src_name}]({ds$src_url}){x_spatial}")
    
    tbl <- x_df %>% 
      kable(caption = x_caption, format = "pipe")
  }
  
  tbl
}
```


# UPDATED tabulate_dataset_shp_within_aoi

```{r new-func, eval=F, class.source='fold-show'}
tabulate_dataset_shp_within_aoi <- function(dataset_code, aoi_wkt, output = "kable"){
  # summarize shapefile dataset from area of interest

  message(glue("tab..._shp_within_aoi(dataset_code='{dataset_code}', aoi_wkt='{paste(aoi_wkt, collapse=';')}')"))

  if (is.null(aoi_wkt))
    return("Please draw a Location to get a summary of the intersecting features for this dataset.")

  ds <- tbl(con, "datasets") %>% 
    filter(code == !!dataset_code) %>% 
    replace_na(list(buffer_nm = 0)) %>% 
    collect()
  
  if (length(aoi_wkt) > 1){
    aoi_wkts <- glue("'SRID=4326;{aoi_wkt}'::geometry")
    aoi_sql  <- glue("ST_COLLECT(\n{paste(aoi_wkts, collapse=',\n')})") # Is this recreating the ST_COLLECT statement
                                                                        # for every item in <aoi_wkt> array?
  } else {
    # aoi_sql <- glue("'SRID=4326;{aoi_wkt}'")
    aoi_sql <- glue("'SRID=4326;{aoi_wkt}'::geometry")
  }

  # Use CTE instead of temporary tables
  # TODO
  #    Add conditional to check if ds$summarize_r
  #    Drop geometry column in x_df?
  if (!is.na(ds$sumarize_sql)){
    x_df <- dbGetQuery(con,
              glue("with 
                    tmp_selarea as (
                      select ST_BUFFER({aoi_sql}, {ds$buffer_nm}) as geom
                    )
                    {ds$select_sql} as ds
                    inner join tmp_selarea on ST_INTERSECTS(ds.geometry, tmp_selarea.geom)
                    "))
  } else {
    x_df <- dbGetQuery(con,
              glue("with 
                    tmp_selarea as (
                      select ST_BUFFER({aoi_sql}, {ds$buffer_nm}) as geom
                    ),
                    tmp_aoi as (
                      {ds$select_sql} as ds
                      inner join tmp_selarea on ST_INTERSECTS(ds.geometry, tmp_selarea.geom)
                    )
                   {ds$summarize_sql}
                   "))
  }
  if (output == "tibble"){
    return(x_df)
  }

  x_spatial <- ifelse(
    ds$buffer_nm == 0,
    glue("\n\nSpatial: within site", .trim = F),
    glue("\n\nSpatial: within {ds$buffer_nm} nautical miles of site", .trim = F))

  if (knitr::is_html_output()){
    x_caption <- HTML(markdownToHTML(
      text = glue("Source: [{ds$src_name}]({ds$src_url}){x_spatial}"),
      fragment.only = T))
    
    tbl <- x_df %>% 
      kbl(caption = x_caption) %>%
      kable_styling(
        # full_width = F, position = "left", # position = "float_right"
        bootstrap_options = c("striped", "hover", "condensed", "responsive"))
    
  } else {
    x_caption <- glue("Source: [{ds$src_name}]({ds$src_url}){x_spatial}")
    
    tbl <- x_df %>% 
      kable(caption = x_caption, format = "pipe")
  }
  
  tbl
}
```


# Speed Testing

Gathered from `/share/github/apps/logging_tmp.txt`:

```yaml
title: Spatial Speed Test
aoi_wkt: POLYGON ((-122.0821 33.74988, -118.3637 32.54795, -117.1548 34.66526, -120.5
  35.57111, -122.0821 33.74988))
technology: Marine Energy
lit_queries:
- Fish AND Marine Energy
- Marine Mammals AND Marine Energy
spatial_receptors:
- Fish
- Marine Mammals
lit_tethys: yes
lit_ferc: ~
spatial: yes
mgt_tethys: ~
```

```{r}
source("functions.R")

params <- list(
  spatial_receptors = c("Fish", "Marine Mammals"),
  aoi_wkt           = "POLYGON ((-122.0821 33.74988, -118.3637 32.54795, -117.1548 34.66526, -120.5
  35.57111, -122.0821 33.74988))")

datasets <- tbl(con, "datasets") %>% 
  collect() %>%
  filter(ready) %>% 
  replace_na(list(buffer_km = 0)) %>% 
  select(-notes, -issues) %>% 
  separate_rows(tags, sep = ";") %>% 
  rename(tag = tags) %>% 
  mutate(
    tag = str_trim(tag)) %>% 
  filter(
    tag %in% params$spatial_receptors) %>% 
  arrange(tag, title)

datasets

source("functions.R")
system.time({
  for (i in 1:nrow(datasets)){
    message(glue("dataset: {datasets$code[i]}"))
    res <- tabulate_dataset_shp_within_aoi2(datasets$code[i], params$aoi_wkt)
  }
})
# before adding datasets: hapc, cetacean-pacific-summer
#   tabulate_dataset_shp_within_aoi: 10.498 12.029 10.846 
#   tabulate_dataset_shp_within_aoi2: 0.551 0.555
# after adding datasets: hapc, cetacean-pacific-summer
#   tabulate_dataset_shp_within_aoi: 36.069 35.531
#   tabulate_dataset_shp_within_aoi2: 2.849 2.236
```

20x improvement!



# Intersects vs Intersection

This was a pain to figure out.

ST_Area(geometry) for SRID:4326 returns area in square meters.

```{r intersection_intersects}
if (!require(librarian)){
  remotes::install_github("DesiQuintans/librarian")
  library(librarian)
}
shelf(
  # database
  DBI, RPostgres,
  # spatial
  ggmap, leaflet,
  r-spatial/mapview, # https://github.com/r-spatial/mapview/issues/324
  sf, sp,
  # scrape
  rvest,
  # tidyverse
  dplyr, purrr, readr, tibble, tidyr,
  # someday
  # googledrive, zeallot,
  # report
  DT, gt, htmltools, htmlwidgets, kableExtra, knitr, markdown, rmarkdown, shiny, webshot,
  # utility
  fs, glue, here, png, scales, stringr, urltools)
here <- here::here


pass <- readLines("/share/.password_mhk-env.us")
con = dbConnect(
  RPostgres::Postgres(),
  dbname   = "gis",
  host     = "marineenergy.app",
  port     = 5432,
  user     = "admin",
  password = pass
)


x_df <- dbGetQuery(con,
          glue("with 
  tmp_selarea as (
  select ST_BUFFER(
    st_collect(
  	'SRID=4326;
  	 POLYGON ((-123.99169921875 43.89789239125797, 
  		-124.91455078125 43.18114705939968, 
  		-125.52978515625001 43.46886761482925, 
  		-125.44189453125 42.84375132629021, 
  		-125.96923828125 42.5530802889558, 
  		-125.88134765625 42.09822241118974, 
  		-125.0244140625 42.24478535602799, 
  		-124.91455078125 42.65012181368022, 
  		-124.47509765625 42.5530802889558, 
  		-124.07958984375001 42.84375132629021, 
  		-123.99169921875 43.89789239125797))'::geometry),
  	0) as geom ),
  
  
  tmp_seldata as (
    select
      unnest(
        array['Blue whale', 'Fin whale', 'Berardius', 'Short-beaked common dolphin',
                  'Rissos dolphin', 'Pacific white-sidded dolphin', 'Northern right whale dolphin', 'Humpback whale',
                  'Dalls porpoise', 'Sperm whale', 'Striped dolphin', 'Small beaked whale']
        ) as common_name,
      unnest(
        array['Balaenoptera musculus', 'Balaenoptera physalus', 'Berardius', 'Delphinus delphis',
                  'Grampus griseus', 'Lagenorhynchus obliquidens', 'Lissodelphis borealis', 'Megaptera novaeangliae',
                  'Phocoenoides dalli', 'Physeter macrocephalus', 'Stenella coeruleoalba', 'Ziphius and Mesoplodon']
        ) as taxon,
      unnest(
        array[bmu_i_u_d, bph_r_u_d, bba_r_u_d, dde_i_u_d, 
                  ggr_r_u_d, lob_i_u_d, lbo_i_u_d, mno_r_u_d, 
                  pda_i_u_d, pma_i_u_d, sco_r_u_d, zsm_r_u_d]
        ) as avg_density,
    	geometry,
    	st_intersection(ds.geometry, tmp_selarea.geom) as ixn_msq
    from 
      shp_serdp_swfsc_cce as ds, tmp_selarea
    where 
      st_intersects(ds.geometry, tmp_selarea.geom)
  ),
  
  
  tmp_ixn as (
    select
    	common_name,
    	taxon,
    	ds.geometry,
    	ds.area_int,
    	st_area(ds.geometry::geography) as ixs_msq,   -- whole features
    	st_area(ds.area_int::geography) as ixn_msq, -- clipped features
    	st_area(tmp_selarea.geom::geography) as aoi_msq,     -- area of interest
    	100 * (st_area(ds.area_int::geography) / st_area(ds.geometry::geography)) as pct_overlap, -- percent of feature that remains after intersection
    	avg_density, -- feature density value
    	avg_density * (st_area(ds.area_int::geography) / st_area(ds.geometry::geography)) as new_avg_density -- area weighted density
    from 
      tmp_seldata as ds, tmp_selarea
  )
  
  
  select 
    common_name, taxon, sum(avg_density) as avg_density, sum(new_avg_density) as area_weighted_avg_density
  from
    tmp_ixn
  group by 
    common_name, taxon
  order by
    common_name, taxon
  "))

x_df
```


---


# Generalize CTE
```{r generalize_cte, eval=F, class.source='fold-show'}

  if (ds$st_intersection = T){
    # Add st_intersection() to data selection query by finding the "geometry"
    #  column and appending st_intersection()
    selsql <- str_replace({ds$select_sql}, 'geometry', 'geometry, st_intersection(ds.geometry, tmp_selarea.geom) as ixn ')
    
    x_df <- dbGetQuery(con,
              glue("with
                      tmp_selarea as (
                        select ST_BUFFER({aoi_sql}, {ds$buffer_nm}) as geom
                      ),
                      tmp_aoi as (
                        {selsql} as ds, tmp_selarea
                        where st_intersects(ds.geometry, tmp_selarea.geom)
                      )
                    {ds$summarize_sql}
                    "))
  } else {
    if (!is.na(ds$sumarize_sql)){
      x_df <- dbGetQuery(con,
                glue("with 
                      tmp_selarea as (
                        select ST_BUFFER({aoi_sql}, {ds$buffer_nm}) as geom
                      )
                      {ds$select_sql} as ds
                      inner join tmp_selarea on ST_INTERSECTS(ds.geometry, tmp_selarea.geom)
                      "))
    } else {
      x_df <- dbGetQuery(con,
                glue("with 
                      tmp_selarea as (
                        select ST_BUFFER({aoi_sql}, {ds$buffer_nm}) as geom
                      ),
                      tmp_aoi as (
                        {ds$select_sql} as ds
                        inner join tmp_selarea on ST_INTERSECTS(ds.geometry, tmp_selarea.geom)
                      )
                     {ds$summarize_sql}
                     "))
    }
  }

```


# tabulate_dataset_shp_within_aoi3
```{r tabulate_dataset_shp_within_aoi3, eval=F, class.source='fold-show'}
tabulate_dataset_shp_within_aoi3 <- function(dataset_code, aoi_wkt, output = "kable"){
  # summarize shapefile dataset from area of interest
  
  # dataset_code = "cetacean-bia"; aoi_wkt = params$aoi_wkt; output = "kable"
  # dataset_code = "cetacean-pacific-summer"; aoi_wkt = params$aoi_wkt; output = "kable"

  message(glue("tab..._shp_within_aoi(dataset_code='{dataset_code}', aoi_wkt='{paste(aoi_wkt, collapse=';')}')"))

  if (is.null(aoi_wkt))
    return("Please draw a Location to get a summary of the intersecting features for this dataset.")

  ds <- tbl(con, "datasets") %>% 
    filter(code == !!dataset_code) %>% 
    replace_na(list(buffer_nm = 0)) %>% 
    collect()
  
  if (length(aoi_wkt) > 1){
    aoi_wkts <- glue("'SRID=4326;{aoi_wkt}'::geometry")
    aoi_sql  <- glue("ST_COLLECT(\n{paste(aoi_wkts, collapse=',\n')})") # Is this recreating the ST_COLLECT statement
                                                                        # for every item in <aoi_wkt> array?
  } else {
    # aoi_sql <- glue("'SRID=4326;{aoi_wkt}'")
    aoi_sql <- glue("'SRID=4326;{aoi_wkt}'::geometry")
  }

  # Different set of queries required for data sets that do or
  #   do not need area weighted statistics 
  if (ds$st_intersection){
    # Area weighted statistics ARE required
    ixn_sql <- str_replace({ds$select_sql}, 'geometry', 'geometry, st_intersection(ds.geometry, buf_aoi.geom) as ixn ')

    if (!is.na(ds$summarize_sql)){
      x_df <- dbGetQuery(
        con,
        glue("
          with
            buf_aoi as (
              select ST_BUFFER({aoi_sql}, {ds$buffer_nm}) as geom),
            tmp_aoi as (
              {ixn_sql} as ds, buf_aoi
              where st_intersects(ds.geometry, buf_aoi.geom))
            {ds$summarize_sql}
          "))
    } else {
      x_sf <- st_read(
        con, 
        glue("
          with
            buf_aoi as (
              select ST_BUFFER({aoi_sql}, {ds$buffer_nm}) as geom)
            {ixn_sql} as ds, buf_aoi
            where st_intersects(ds.geometry, buf_aoi.geom)
          "))
      x_df <- st_drop_geometry(x_sf)
      
      if (!is.na(ds$summarize_r))
        eval(parse(text=ds$summarize_r))
    }
    
  } else {
    # Area weighted statistics NOT required
    if (!is.na(ds$summarize_sql)){
      x_df <- dbGetQuery(
        con, glue("
          with 
            buf_aoi as (
              select ST_BUFFER({aoi_sql}, {ds$buffer_nm}) as geom ),
            tmp_aoi as (
              {ds$select_sql} as ds
              inner join buf_aoi on st_intersects(ds.geometry, buf_aoi.geom) )
           {ds$summarize_sql}
           "))
    } else {
      x_sf <- st_read(
        con, query = glue("
          with 
            buf_aoi as (
              select ST_BUFFER({aoi_sql}, {ds$buffer_nm} * 1852) as geom)
            {ds$select_sql} as ds
            inner join buf_aoi on st_intersects(ds.geometry, buf_aoi.geom )
            "))
      x_df <- st_drop_geometry(x_sf)
      
      if (!is.na(ds$summarize_r))
        eval(parse(text=ds$summarize_r))
    }
  }
  if (output == "tibble"){
    return(x_df)
  }

  x_spatial <- ifelse(
    ds$buffer_nm == 0,
    glue("\n\nSpatial: within site", .trim = F),
    glue("\n\nSpatial: within {ds$buffer_nm} nautical miles of site", .trim = F))

  if (knitr::is_html_output()){
    x_caption <- HTML(markdownToHTML(
      text = glue("Source: [{ds$src_name}]({ds$src_url}){x_spatial}"),
      fragment.only = T))
    
    tbl <- x_df %>% 
      kbl(caption = x_caption) %>%
      kable_styling(
        # full_width = F, position = "left", # position = "float_right"
        bootstrap_options = c("striped", "hover", "condensed", "responsive"))
    
  } else {
    x_caption <- glue("Source: [{ds$src_name}]({ds$src_url}){x_spatial}")
    
    tbl <- x_df %>% 
      kable(caption = x_caption, format = "pipe")
  }
  
  tbl
}
```

```{r testing_func}
# source("functions.R")
# datasets_gsheet2db() # update datasets table from spatial Google Sheet

params <- list(
  spatial_receptors = c("Cetaceans"),
  aoi_wkt           = "POLYGON ((-123.99169921875 43.89789239125797, 
  		-124.91455078125 43.18114705939968, 
  		-125.52978515625001 43.46886761482925, 
  		-125.44189453125 42.84375132629021, 
  		-125.96923828125 42.5530802889558, 
  		-125.88134765625 42.09822241118974, 
  		-125.0244140625 42.24478535602799, 
  		-124.91455078125 42.65012181368022, 
  		-124.47509765625 42.5530802889558, 
  		-124.07958984375001 42.84375132629021, 
  		-123.99169921875 43.89789239125797))")

datasets <- tbl(con, "datasets") %>% 
  collect() %>%
  filter(ready) %>% 
  replace_na(list(buffer_km = 0)) %>% 
  select(-notes, -issues) %>% 
  separate_rows(tags, sep = ";") %>% 
  rename(tag = tags) %>% 
  mutate(
    tag = str_trim(tag)) %>% 
  filter(
    tag %in% params$spatial_receptors) %>% 
  arrange(tag, title)


source("functions.R")

# All datasets matching tags
for (i in 1:nrow(datasets)){
  message(glue("dataset: {datasets$code[i]}"))
  res <- tabulate_dataset_shp_within_aoi3(datasets$code[i], params$aoi_wkt)
  res
}

# First dataset matching tags
message(glue("dataset: {datasets$code[i]}"))
res <- tabulate_dataset_shp_within_aoi3(datasets$code[i], params$aoi_wkt)
res
```