Load vector dataset: Essential Fish Habitat

Read dataset

library(here)
library(glue)
library(readr)
library(fs)
library(sf)
library(raster)
library(mapview)
library(leaflet)
library(DT)
library(dplyr)

dir_mc       <- '/share/data/marinecadastre.gov'
csv_mc       <- file.path(dir_mc, '_datasets.csv')
csv_mc_paths <- file.path(dir_mc, '_datasets_paths.csv')

d <- read_csv(csv_mc_paths) %>% 
  filter(title == "Essential Fish Habitat (EFH)")

# get shapefile(s)
shp <- list.files(d$dir_path, ".*\\.shp$", recursive = T, full.names = T)

# read in shapefile
d_ply <- read_sf(shp)

Connect to database

# connect to database
library(DBI)
library(RPostgres)

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

tbls <- dbListTables(con)
tbls

## [1] "geography_columns" "geometry_columns"  "spatial_ref_sys"  
## [4] "raster_columns"    "raster_overviews"  "d_cetmap_bia"     
## [7] "topology"          "layer"             "d_efh"

# helper functions
drop_d <- function(d_tbl){
  dbSendQuery(con, glue("SELECT DropGeometryTable ('public','{d_tbl}');"))
}

dbSafeNames = function(names) {
  # make names db safe: no '.' or other illegal characters,
  # all lower case and unique
  names = gsub('[^a-z0-9]+','_',tolower(names))
  names = make.names(names, unique=TRUE, allow_=TRUE)
  names = gsub('.','_',names, fixed=TRUE)
  names
}

dbRenameTable <- function(con, old, new){
  if (old %in% dbListTables(con))
    dbExecute(con, glue("ALTER TABLE {old} RENAME TO {new}"))
}

# rename old table
dbRenameTable(con, "cetmap_bia", "d_cetmap_bia")

Load dataset

# set field names
names(d_ply) <- dbSafeNames(names(d_ply))

d_tbl  <- "d_efh" # prefix: mc_ for MarineCadastre, ds_ for dataset?
d_redo <- FALSE

# project to geographic coordinate reference system if need be
if (st_crs(d_ply) != st_crs(4326)){
  d_ply <- st_transform(d_ply, crs = 4326)
}

if (!d_tbl %in% tbls | d_redo){
  st_write(d_ply, con, d_tbl)
}
dbListTables(con)

## [1] "geography_columns" "geometry_columns"  "spatial_ref_sys"  
## [4] "raster_columns"    "raster_overviews"  "d_cetmap_bia"     
## [7] "topology"          "layer"             "d_efh"

Spatial query

library(leaflet)
library(lwgeom)
library(dplyr)

aoi_geo <- here("report/data/geom.geojson")

aoi_ply <- read_sf(aoi_geo)
aoi_wkt <- st_as_text(aoi_ply$geometry, EWKT=T)

d_cols <- tbl(con, d_tbl) %>% head(0) %>% collect() %>% names()
#d_cols_notgeom <- setdiff(d_cols, c("geometry", "geom"))

#dist_m  <- 10*1000
dist_m  <- 0
sql <- glue("
  SELECT {paste(d_cols, collapse=',')}
  FROM public.{d_tbl} d
  WHERE ST_DWithin(Geography(d.geometry), '{aoi_wkt}', {dist_m});")
d_win <- st_read(con, query = sql)

bb <- st_bbox(aoi_ply) %>% as.vector()

leaflet() %>% 
  addProviderTiles(providers$Esri.OceanBasemap, group = "Esri Ocean") %>% 
  addProviderTiles( providers$Stamen.TonerLite, group = "Toner Lite") %>%
  addPolygons(data = aoi_ply, color =    "red", group = "AOI") %>% 
  addPolygons(data =   d_win, color =  "green", group = "EFH w/in 0km") %>% 
  #addPolygons(data =   d_ply, color =   "#03F", group = "EFHs") %>% 
  addScaleBar() %>% 
  fitBounds(bb[1], bb[2], bb[3], bb[4]) %>% 
  addLayersControl(
    baseGroups = c("Esri Ocean", "Toner Lite"),
    #overlayGroups = c("AOI", "EFH w/in 0km", "EFHs"),
    overlayGroups = c("AOI", "EFH w/in 0km"),
    options = layersControlOptions(collapsed = FALSE)) # %>% 

  #hideGroup("EFHs")

d_win %>% 
  st_drop_geometry() %>% 
  datatable()