Explore PostGIS ST_Contour()

defaults

librarian::shelf(
  glue, here, leaflet, lubridate)

dir_apps <- here("../apps_dev-demo")
source(glue("{dir_apps}/libs/db.R"))
source(glue("{dir_apps}/libs/functions.R"))

ctdcast_dates <- dbGetQuery(
  con, "SELECT MIN(date) min, MAX(date) max FROM ctd_casts")

contour_args <- list(
  variable      = "ctd_bottles.t_degc",
  value         = "avg",
  aoi_pattern   = c("standard", "extended"),
  aoi_shore     = c("nearshore", "offshore"),
  date_beg      = as.Date(today() - dyears(5)),
  date_end      = ctdcast_dates$max,
  date_qrtr     = c(1, 2),
  depth_m_min   = 0,            
  depth_m_max   = 515,                 
  n_bins        = 7)

map_r <- function(r){
  rng_r <- global(r, "range", na.rm=T) |> as.numeric()
  rng_p <- range(pts$z)
  rng <- range(c(rng_r, rng_p))
 
  pal <- colorNumeric("Spectral", rng, na.color = NA)
  
  map_base(base_opacity = 0.2) |> 
    addRasterImage(
      r, colors = pal, opacity=0.7, group = "raster") |>
    addLegend(
      pal = pal, values = rng) |> 
    addCircleMarkers(
      data   = pts,
      radius = 5, fillOpacity = 0.2,
      stroke = T, weight = 1, opacity = 0.8,
      color  = ~pal(z), group = "points") |> 
    addLayersControl(
      # baseGroups    = c("OceanBase", "OceanLabels"),
      overlayGroups = c("raster", "points"),
      options       = layersControlOptions(collapsed = T),
      position      = "topleft") |> 
    hideGroup("points")
}

interp_rast <- function(algorithm){
  do.call(
  get_contour,
  c(contour_args, list(
    return_type   = "raster",
    idw_algorithm = algorithm))) |> 
  map_r()
}

inputs <- do.call(
  get_contour,
  c(contour_args, list(
    return_type   = "inputs")))
st_bbox(inputs$ext)

      xmin       ymin       xmax       ymax 
-127.56206   28.88526 -116.61902   38.69428 

inputs |> 
  st_drop_geometry() |> 
  select(-ext)

  pixelsize                       algorithm
1       0.1 invdist:power:2.0:smoothing:2.0

sizes <- do.call(
  get_contour,
  c(contour_args, list(
    return_type   = "sizes")))
sizes

  width height upperleftx upperlefty
1   110     99  -127.5621   38.69428

points

pts <- do.call(
  get_contour,
  c(contour_args, list(
    return_type   = "points")))

pal <- colorNumeric("Spectral", pts$z)

map_base(base_opacity=0.1) |>
  addCircleMarkers(
    data   = pts,
    radius = 5, fillOpacity = 0.2,
    stroke = T, weight = 1, opacity = 0.8,
    color  = ~pal(z))

invdist

• Interpolation algorithms - gdal_grid — GDAL documentation

interp_rast("invdist")

invdistnn

interp_rast("invdistnn")

average

inputs <- do.call(
  get_contour,
  c(contour_args, list(
    return_type   = "inputs",
    idw_algorithm = "average:radius:10.0")))

source(glue("{dir_apps}/libs/functions.R"))
sql <- do.call(
  get_contour,
  c(contour_args, list(
    return_type   = "sql",
    idw_algorithm = "average:radius:10.0")))
cat(sql)

INSERT INTO z_idw (
  args_hash, 
  args_json,
  rast)
WITH 
aoi AS (
  SELECT ST_Union(geom) AS geom
  FROM effort_zones 
  WHERE
  sta_pattern IN ('standard','extended') AND
  sta_shore   IN ('nearshore','offshore') ),
pts AS (
  SELECT 
  ST_SetSRID(
    ST_MakePoint(
      ST_X(c.geom),
      ST_Y(c.geom),
      AVG(t_degc) ),
    4326) AS geom
  FROM ctd_casts AS c 
  JOIN ctd_bottles AS cb USING (cast_count)
  JOIN aoi ON ST_Covers(aoi.geom, c.geom)
  WHERE 
  (depthm BETWEEN 0 AND 515) AND
  (date BETWEEN '2018-03-26' AND '2020-01-26') AND
  -- DATE_PART('year'   , date) = 2020  AND
  DATE_PART('quarter', date) IN (1,2)
  GROUP BY c.geom ),
inputs AS (
  SELECT
  0.1::float8 AS pixelsize,
  -- https://gdal.org/programs/gdal_grid.html#interpolation-algorithms
  -- 'invdist:power:5.5:smoothing:2.0' AS algorithm,
  'average:radius:10.0' AS algorithm, -- default parameters
  ST_Collect(pts.geom) AS geom,
  ST_Expand(ST_Collect(aoi.geom), 0.5) AS ext
  FROM aoi, pts),
-- Calculate output raster geometry
-- Use the expanded extent to take in areas beyond the limit of the
-- temperature stations
sizes AS (
  SELECT
  ceil((ST_XMax(ext) - ST_XMin(ext))/pixelsize)::integer AS width,
  ceil((ST_YMax(ext) - ST_YMin(ext))/pixelsize)::integer AS height,
  ST_XMin(ext) AS upperleftx,
  ST_YMax(ext) AS upperlefty
  FROM inputs)
SELECT 
  '06a34cf6' AS args_hash,
  '{"variable":["ctd_bottles.t_degc"],"value":["avg"],"aoi_pattern":["standard","extended"],"aoi_shore":["nearshore","offshore"],"date_beg":["2018-03-26"],"date_end":["2020-01-26"],"date_qrtr":[1,2],"depth_m_min":[0],"depth_m_max":[515],"n_bins":[7],"idw_algorithm":["average:radius:10.0"]}' AS args_json,
  ST_Clip(
    ST_SetBandNoDataValue(
      ST_InterpolateRaster(
        geom,
        algorithm,
        ST_SetSRID(
          ST_AddBand(
            ST_MakeEmptyRaster(
              width, height, upperleftx, upperlefty, pixelsize), 
            '32BF'), 
          ST_SRID(geom))),
      -9999),
    (SELECT geom FROM aoi)) AS rast
FROM sizes, inputs;

interp_rast("average:radius:0.1")

nearest

interp_rast("nearest")

linear

interp_rast("linear")