explore_pg_contour2.qmd

---
title: "Explore PostGIS ST_Contour()"
format:
  html:
    toc: true
    toc-depth: 3
    code-fold: true
    self-contained: true
editor_options: 
  chunk_output_type: console
execute:
  warning: false
---


## defaults

```{r}
librarian::shelf(
  glue, here, leaflet, lubridate, mapview, purrr)

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)

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

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

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

## points

```{r}
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))
```

* [Interpolation algorithms - gdal_grid — GDAL documentation](https://gdal.org/programs/gdal_grid.html#interpolation-algorithms)

## `stars::st_contour()`

```{r}
librarian::shelf(
  gstat, stars)

# template raster from buffered area of interest
buffer_dd    <- 0.5
aoi_b <- st_buffer(aoi, buffer_dd)
b <- st_bbox(aoi_b)
dx <- b$xmax - b$xmin
dy <- b$ymax - b$ymin
# calculate pixel size of at least 100 pixels for max dim (x or y of aoi)
d100 <- max(c(dx, dy)) / 100
brks <- c(-Inf, 0.01, 0.1, 1, 10, Inf)
d <- cut(d100, brks, include.lowest=T)
pxl_size <- brks[which(a==levels(d))]
pxl_size <- case_match(
    pxl_size,
    -Inf ~  0.01,  # min
     Inf ~ 10.00,  # max
    .default = pxl_size)
nx <- ceiling( dx / pxl_size )
ny <- ceiling( dy / pxl_size )
# create template raster of extent
r_ext <- rast(
  nrows = nx, 
  ncols = ny,
  xmin  = b$xmin,
  xmax  = b$xmax,
  ymin  = b$ymin,
  ymax  = b$ymax,
  crs = "EPSG:4326")
suppressWarnings({
  r_ext <- st_as_stars(r_ext) })

# alg <- "invdistnn:power=1.3:smoothing=0.2:radius=2.0:max_points=12:nodata=-9999.0"
r_i <- gstat::idw(
  z ~ 1, locations = pts, newdata = r_ext)

g <- gstat(z ~ x + y, idp = 0.5, data = pts)
r_i <- predict(g, r_ext)
r_i <- predict(g, r_ext, BLUE = T)

plot(r_i)
names(r_i)
r_i |> 
  slice("band", "var1.pred")
r_i[[1]]

r_ic <- r_i |> 
  rast() |> 
  mask(aoi_b) |> 
  trim() |> 
  st_as_stars()
ply <- st_contour(r_ic) |> 
  st_intersection(aoi)

mapView(r_ic) + 
  mapView(ply, zcol = "Max")
```

## `terra::interpIDW` -> `stars::st_contour()`

```{r}
librarian::shelf(
  gstat, stars, terra)
select <- dplyr::select

contour_args <- list(
  variable      = "ctd_bottles.t_degc",
  value         = "avg",
  # aoi_pattern   = c("standard", "extended"),
  # aoi_shore     = c("nearshore", "offshore"),
  aoi_pattern   = c("standard"),
  aoi_shore     = c("nearshore"),
  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)

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

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

# template raster from buffered area of interest
b <- st_bbox(aoi)
dx <- b$xmax - b$xmin
dy <- b$ymax - b$ymin
# calculate pixel size approaching 100 pixels for max dim (x or y of aoi)
d100 <- max(c(dx, dy)) / 100
brks <- c(-Inf, 0.01, 0.1, 1, 10, Inf)
d <- cut(d100, brks, include.lowest=T)
pxl_size <- brks[which(a==levels(d))]
pxl_size <- case_match(
    pxl_size,
    -Inf ~  0.01,  # min
     Inf ~ 10.00,  # max
    .default = pxl_size)
nx <- ceiling( dx / pxl_size )
ny <- ceiling( dy / pxl_size )
# create template raster of extent
r_ext <- rast(
  nrows = nx, 
  ncols = ny,
  xmin  = b$xmin,
  xmax  = b$xmax,
  ymin  = b$ymin,
  ymax  = b$ymax,
  crs = "EPSG:4326") |> 
  extend(c(2, 2)) # extend by 2 pixels in x & y

m_pts <- pts |> 
  mutate(
    x = st_coordinates(geom)[,"X"],
    y = st_coordinates(geom)[,"Y"]) |> 
  st_drop_geometry() |>
  select(x, y, z) |> 
  as.matrix()

# alg <- "invdistnn:power=1.3:smoothing=0.2:radius=2.0:max_points=12:nodata=-9999.0"
r_i <- interpIDW(
  x = r_ext, y = m_pts, 
  radius=2.0, power=1.3, smooth=0.2, 
  maxPoints=12, minPoints=1, near=FALSE, fill=NA)
# plot(r_i)

aoi_b <- aoi |> 
  st_buffer(pxl_size * 2)

r_ic <- r_i |> 
  mask(aoi_b) |> 
  trim()

brks = classInt::classIntervals(
  values(r_ic, na.rm=T), n=contour_args$n_bins+1, 
  style = "pretty")$brks

ply <- st_contour(st_as_stars(r_ic), breaks = brks) |> 
  st_intersection(aoi) |> 
  mutate(
    avg = (Min + Max) / 2 )

mapView(r_ic) + 
mapView(ply, zcol = "avg")
# 
```

TODO:
- show legend with discrete breaks

## `automap::autofitVariogram()`

```{r}
librarian::shelf(automap)

xy <- st_coordinates(pts$geom)

xy
pts <- pts |> 
  mutate(
    x = st_coordinates(pts$geom)[,"X"],
    y = st_coordinates(pts$geom)[,"Y"])
v_mod_ok = autofitVariogram(z ~ x + y, pts)
plot(v_mod_ok)
g = gstat(formula = z ~ x + y, model = v_mod_ok$var_model, data = pts)

terra::interpIDW()
terra::colFromX()
r_ext
r_k = predict(g, r_ext)
plot(r_k)

mapView(r_k)
r_i <- gstat::idw(
  z ~ 1, locations = pts, newdata = r_ext, idp = 0.5)
plot(r_i)
names(r_i)
r_i |> 
  slice("band", "var1.pred")
r_i[[1]]

r_ic <- r_i |> 
  rast() |> 
  mask(aoi_b) |> 
  trim() |> 
  st_as_stars()

v_ic <- values(rast(r_ic)[[1]], na.rm=T)
brks = classInt::classIntervals(
  v_ic, n=contour_args$n_bins+1, style = "pretty")$brks

ply <- st_contour(r_ic, breaks = brks) |> 
  st_intersection(aoi)

mapView(r_ic) + 
  mapView(ply, zcol = "Max")







librarian::shelf(
  raster, rgdal, sp, tmap)

# Load precipitation data
z <- gzcon(url("https://github.com/mgimond/Spatial/raw/main/Data/precip.rds"))
P <- readRDS(z)

# Load Texas boudary map
z <- gzcon(url("https://github.com/mgimond/Spatial/raw/main/Data/texas.rds"))
W <- readRDS(z)

# Replace point boundary extent with that of Texas
P@bbox <- W@bbox

# library(gstat) # Use gstat's idw routine
# library(sp)    # Used for the spsample function

# Create an empty grid where n is the total number of cells
grd              <- as.data.frame(sp::spsample(P, "regular", n=50000))
names(grd)       <- c("X", "Y")
coordinates(grd) <- c("X", "Y")
gridded(grd)     <- TRUE  # Create SpatialPixel object
fullgrid(grd)    <- TRUE  # Create SpatialGrid object

# Add P's projection information to the empty grid
proj4string(P) <- proj4string(P) # Temp fix until new proj env is adopted
proj4string(grd) <- proj4string(P)

# Interpolate the grid cells using a power value of 2 (idp=2.0)
P.idw <- gstat::idw(Precip_in ~ 1, P, newdata=grd, idp=2.0)

# Convert to raster object then clip to Texas
r       <- raster(P.idw)
r.m     <- mask(r, W)

# Plot
tm_shape(r.m) + 
  tm_raster(n=10,palette = "RdBu", auto.palette.mapping = FALSE,
            title="Predicted precipitation \n(in inches)") + 
  tm_shape(P) + tm_dots(size=0.2) +
  tm_legend(legend.outside=TRUE)


?st_warp

gstat::idw()

stars::st_as_stars(pts)

brks = classInt::classIntervals(
  v, n=contour_args$n_bins+1, style = "pretty")$brks



stars::st_contour()
pts

```


## `invdist`

```{r}
interp_rast("invdist")
```

## `invdist` tweaked

```{r}
interp_rast("invdist:power=1.5:smoothing=0.2")
```

## `invdistnn`

```{r}
interp_rast("invdistnn")
```

## `invdistnn` tweaked

```{r}
alg <- "invdistnn:power=1.3:smoothing=0.2:radius=2.0:max_points=12:nodata=-9999.0"

r <- do.call(
  get_contour,
  c(contour_args, list(
    idw_algorithm = alg,
    return_type   = "raster")))

p_eqbrks <- do.call(
  get_contour,
  c(contour_args, list(
    idw_algorithm = alg,
    return_type   = "polygon")))
attr(p_eqbrks, "breaks")
# [1]  9 10 11 12 13 14 15 16
mapView(p_eqbrks, zcol = "k_avg")
mapView(p_eqbrks, zcol = "val_ctr")

clns_eqbrks <- do.call(
  get_contour,
  c(contour_args, list(
    idw_algorithm = alg,
    return_type   = "closed_lines")))
mapView(clns_eqbrks)

p0 <- p

source(glue("{dir_apps}/libs/functions.R"))
p <- do.call(
  get_contour,
  c(contour_args, list(
    idw_algorithm = alg,
    # return_type   = "polygons")))
    return_type   = "polygons3")))

mapView(p, zcol = "val_ctr") + 
  mapView(r)

st_drop_geometry(p) |> 
  View()

source(glue("{dir_apps}/libs/functions.R"))
p <- do.call(
  get_contour,
  c(contour_args, list(
    idw_algorithm = alg,
    return_type   = "insert_polygons")))

p |> 
  
p |> 
  group_by()
hash = '0b99ea0b'

sql = glue("
  WITH
    plys_ctr_rast AS (
      SELECT * FROM
      idw_plys
      WHERE hash_id = '{hash}'),
    pts_rndm AS (
      SELECT 
        poly_id,
        ST_GeneratePoints(geom, 10, 42) AS geom_multipt -- MULTIPOINT
      FROM plys_ctr_rast
      WHERE oid IN (SELECT generate_series(1, 7)) ),
    multipts_rndm AS (
      SELECT
        poly_id,
        (ST_Dump(geom_multipt)).geom AS geom_pt
      FROM pts_rndm)
    SELECT
      poly_id,
      ST_Value(rast, 1, geom_pt) AS val_rndm,
      geom_pt
    FROM multipts_rndm p, z_idw z
    WHERE 
      z.args_hash = '{hash}' AND
      ST_Intersects(z.rast, p.geom_pt)
  ")
# message(sql)
pts <- st_read(con, query = sql)

sql = glue("
  DO $$
    DECLARE
    row record;
  BEGIN
    FOR row IN (SELECT DISTINCT oid FROM idw_plys WHERE hash_id = '{hash}') LOOP
    
      WITH
        pts_rndm AS (
          SELECT 
            (ST_Dump(ST_GeneratePoints(geom, 10, 42)).geom AS geom
          FROM idw_plys
          WHERE oid = row.oid),
        pts_rast AS (
          SELECT
            ST_Value(rast, 1, geom) AS val_rast
          FROM pts_rndm p, z_idw z
          WHERE 
            z.args_hash = '{hash}' AND
            ST_Intersects(z.rast, p.geom)
      UPDATE idw_plys
        SET val_rndm = (
          SELECT AVG(val_rast) AS val_rndm FROM pts_rast)
        WHERE oid = row.oid

    END LOOP;
  END $$ ;
  ")
message(sql)
dbSendStatement(con, sql)

# message(sql)
pts <- st_read(con, query = sql)

mapView(pts, zcol = "val_rndm")

p <- p |> 
  mutate(
    k_avg   = (k_min + k_max) / 2,
    k_dif    = k_max - k_min)
p


p |> 
  st_drop_geometry() |> 
  View()
mapView(
  p |> 
    filter(k_dif == 1), 
  zcol = "val_ctr",
  alpha.regions = 0.3) + 
mapView(
  p |> 
    filter(k_dif != 1), 
  zcol = "val_ctr",
  alpha.regions = 0.9) + 
  mapView(r)

# TODO: consider dissolving same breaks

map_rp(r, p)
```


## new reclassify approach in R

```{r}
# TODO: polygon dump after reclassifying

source(glue("{dir_apps}/libs/functions.R"))

aoi <- do.call(
  get_contour,
  c(contour_args, list(
    idw_algorithm = alg,
    return_type   = "aoi")))

mapView(r) + 
  mapView(aoi)

r_a <- crop(r, a) |> 
  mask(a)
# TODO: expand by a percent buffer
mapView(r_a)

v <- values(r_a, na.rm=T)
# q <- quantile(v, seq(0, 1, length.out = contour_args$n_bins+1)) |> 
#   discr.round()
q <- pretty(v, n = contour_args$n_bins + 1, bounds = T)

d_q <- tibble(
  from = q[1:(length(q)-1)],
  to   = q[-1]) |> 
  tibble::rowid_to_column("new") |> 
  mutate(
    mean = map2_dbl(from, to, ~mean(c(.x, .y))))
# TODO: show color ramp with `new` values at ticks and full range of min(`from`), max(`to`)

source(glue("{dir_apps}/libs/functions.R"))

r_c <- terra::classify(
  r_a,
  d_q |> 
    select(from, to, new))
mapView(r_a) + 
  mapView(r_c)

librarian::shelf(stars)

x <- stars::st_as_stars(r_a)
v <- na.omit(as.vector(x[[1]]))
range(v)

breaks = classInt::classIntervals(
  v, n=7, style = "pretty")$brks
# "fixed", "sd", "equal", "pretty", "quantile", "kmeans", "hclust", "bclust", "fisher", "jenks", "dpih", "headtails", or "maximum"
breaks
diff(breaks)
r_k <- stars::st_contour(x, breaks=breaks)
#r_k$Min[1] = min(na.omit(as.vector(x[[1]])))
r_k <- r_k |> 
  mutate(
    avg = (Min + Max) / 2)
mapView(r_k, zcol = "avg")

range(na.omit(as.vector(x[[1]]))) # 9.61488723755 15.26715946198




p_r <- terra::as.polygons(r_c) |> st_as_sf(crs=4326) |> 
  left_join(
    d_q, 
    by = c("layer"="new"))
mapView(r_a) + 
  mapView(r_c) + 
  mapView(p_r, zcol="mean") + 
  mapView(a)


write_sf(p_r, "data/test_p_r.geojson", delete_dsn=T)
p_r <- read_ "data/test_p_r.geojson"
p_s <- st_simplify(p_r, preserveTopology = T)
mapView(p_s)
  
plot(r_a)
interp_rast(alg)
```

TODO:

- bootstrap sample and minimize error

## `polygons`

```{r}
l <- do.call(
  get_contour,
  c(contour_args, list(
    idw_algorithm = alg,
    return_type   = "lines")))

mapView(l)
```

```{r}
source(glue("{dir_apps}/libs/functions.R"))
args <- c(contour_args, list(
  idw_algorithm = alg,
  return_type   = "polygons"))
args$n_bins = 7
k <- do.call(
  get_contour, args)


a <- do.call(
  get_contour,
  c(contour_args, list(
    idw_algorithm = alg,
    return_type   = "aoi")))

cl <- do.call(
  get_contour,
  c(contour_args, list(
    idw_algorithm = alg,
    return_type   = "closed_lines")))

source(glue("{dir_apps}/libs/functions.R"))
args <- c(contour_args, list(
  idw_algorithm = alg,
  return_type   = "polygons"))
args$n_bins    = 5
args$rast_redo = T
k <- do.call(get_contour, args)
# mapView(k, zcol = "val_ctr")
# mapView(r)

qbrks <- attr(k, "breaks") |> pull(value) |> sort() |> unique()
qpal <- colorBin("Spectral", bins = qbrks, reverse = T)

map_base() |> 
  clearGroup("contours") |> 
  clearControls() |> 
  addPolygons(
    data = k, fillOpacity = 0.5, 
    weight = 0.5, opacity = 0.7,
    color = ~qpal(val_ctr), group = "contours") |> 
  addLegend(
    pal = qpal, values = k$val_ctr, opacity = 1, 
    title = 'contour', group = "contours") |> 
  addLayersControl(
    baseGroups = c("OceanBase", "OceanLabels"),
    overlayGroups = c("aoi", "contours"),
    options = layersControlOptions(collapsed = T))

table(cut(k$val_ctr, qbrks))

# TODO next:
# - store vectors too
# - breaks option
# - return trimmed raster

sf::write_sf(k, "_k.geojson")

k |> 
  st_drop_geometry() |> 
  View()

k_na <- k |> 
  # filter(is.na(val_avg)) |> 
  filter(poly_id == 72) |> 
  slice(1)
mapView(k, zcol = "val_avg") + 
  mapView(k_na, col.regions="red")

st_centroid(k_na$geom) |> st_coordinates() |> as.numeric()
   

b698cc34

SELECT ST_Value(rast, 1, foo.pt_geom) as val
FROM z_idw z CROSS JOIN (SELECT ST_SetSRID(ST_Point(-119.92052, 34.25617), 4326) As pt_geom) As foo
WHERE z.args_hash = 'b698cc34';

ST_Value()

k <- k |> 
  mutate(
    area_km2 = st_area(geom) |> units::set_units("km^2") |> as.numeric()) |> 
  arrange(desc(area_km2), val_avg)

k |> 
  st_drop_geometry() |> 
  View()

librarian::shelf(
  units)

k$area_km2[1:3] 



k |> 
  filter()

k |> 
  st_drop_geometry() |> 
  View()

r0 <- r
r <- do.call(
  get_contour,
  c(contour_args, list(
    idw_algorithm = alg,
    return_type   = "raster",
    rast_redo     = T)))
mapView(r)

k

l_max <- l |> 
  filter(value == max(l$value))
mapView(a, color="black", col.regions="yellow") + 
  mapView(k) +
  mapView(r) + 
  mapView(l) + 
  mapView(cl) + 
  mapView(l_max, color="red")

k

map_rpk(r, p, k)
```


## `average`

```{r}
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)

interp_rast("average:radius:0.1")
```

## `nearest`

```{r}
interp_rast("nearest")
```

## `linear`

```{r}
interp_rast("linear")
```