Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

speed up oceano app #22

Open
bbest opened this issue Oct 30, 2023 · 1 comment
Open

speed up oceano app #22

bbest opened this issue Oct 30, 2023 · 1 comment

Comments

@bbest
Copy link
Contributor

bbest commented Oct 30, 2023

  • bogs down with increasing time
@bbest
Copy link
Contributor Author

bbest commented Nov 1, 2023

Current approach

Pull all point data from database into R to interpolate to raster before creating contour polygons:

  1. apps:oceano/server.R get_map_data
  2. r_idw <- pts_to_rast_idw(pts, value, aoi, out_tif = idw_tif)
  3. calcofi4r::pts_to_rast_idw()

New approach

Generate interpolated raster in the database and the contour polygons in R using previous attempt (that was successful and fast with interpolation to raster, but wonky with contour creation):

args_in <- list(
variable = variable,
value = value,
aoi_pattern = aoi_pattern,
aoi_shore = aoi_shore,
date_beg = date_beg,
date_end = date_end,
date_qrtr = date_qrtr,
depth_m_min = depth_m_min,
depth_m_max = depth_m_max,
n_bins = n_bins,
idw_algorithm = idw_algorithm)
args_json <- jsonlite::toJSON(args_in)
hash <- digest(args_in, algo="crc32")
message(glue("hash: {hash}"))
# TODO: variable, value; use MATERIALIZED VIEW to combine all vars into single table in advance
aoi_sql <- glue("
aoi AS (
SELECT ST_Union(geom) AS geom
FROM effort_zones
WHERE
sta_pattern IN ('{paste(aoi_pattern, collapse = '\\',\\'')}') AND
sta_shore IN ('{paste(aoi_shore , collapse = '\\',\\'')}') )")
if (return_type == "aoi"){
sql <- glue("
WITH
{aoi_sql}
SELECT geom FROM aoi
")
aoi <- st_read(con, query = sql)
return(aoi)
}
aoi_pts_sql <- glue("
{aoi_sql},
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
(depth_m BETWEEN {depth_m_min} AND {depth_m_max}) AND
(date BETWEEN '{date_beg}' AND '{date_end}') AND
-- DATE_PART('year' , date) = 2020 AND
DATE_PART('quarter', date) IN ({paste(date_qrtr, collapse=',')})
GROUP BY c.geom )")
if (return_type == "points"){
sql <- glue("
WITH
{aoi_pts_sql}
SELECT * FROM pts
")
pts <- st_read(con, query = sql) |>
mutate(
z = st_coordinates(geom)[,"Z"]) |>
st_zm()
return(pts)
}
inputs_sql <- glue("
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,
'{idw_algorithm}' AS algorithm, -- default parameters
ST_Collect(pts.geom) AS geom,
ST_Expand(ST_Collect(aoi.geom), 0.5) AS ext
FROM aoi, pts)")
if (return_type == "inputs"){
sql <- glue("
WITH
{aoi_pts_sql},
{inputs_sql}
SELECT * FROM inputs
")
inputs <- st_read(con, query = sql)
return(inputs)
}
sizes_sql <- glue("
-- 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)")
if (return_type == "sizes"){
sql <- glue("
WITH
{aoi_pts_sql},
{inputs_sql},
{sizes_sql}
SELECT * FROM sizes
")
sizes <- dbGetQuery(con, sql)
return(sizes)
}
# do once
# q <- dbSendStatement(con, "DROP TABLE z_idw"); dbClearResult(q)
# q <- dbSendStatement(
# con,
# "CREATE TABLE z_idw (
# rid SERIAL PRIMARY KEY,
# args_hash TEXT,
# args_json JSON,
# rast RASTER)")
# dbSendStatement(
# con,
# "DROP TABLE idw_plys")
# dbSendStatement(
# con,
# "CREATE TABLE idw_plys (
# oid SERIAL PRIMARY KEY,
# hash_id TEXT,
# poly_id NUMERIC,
# k_min NUMERIC,
# k_max NUMERIC,
# k_avg NUMERIC,
# val_ctr NUMERIC,
# val_rndm NUMERIC,
# geom GEOMETRY(POLYGON, 4326) )")
# dbClearResult(q)
# q <- dbSendStatement(con, "SET postgis.gdal_enabled_drivers = 'ENABLE_ALL'")
# dbClearResult(q)
rast_idw_exists <- tbl(con, "z_idw") |>
filter(args_hash == hash) |>
collect() |>
nrow() > 0
if (rast_idw_exists & rast_redo){
q <- dbSendStatement(con, glue("DELETE FROM z_idw WHERE args_hash = '{hash}'"))
dbClearResult(q)
rast_idw_exists <- F
}
sql <- glue("
INSERT INTO z_idw (
args_hash,
args_json,
rast)
WITH
{aoi_pts_sql},
{inputs_sql},
{sizes_sql}
SELECT
'{hash}' AS args_hash,
'{args_json}' 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;")
if (return_type == "sql"){
return(sql)
}
if (!rast_idw_exists){
# cat(sql)
q <- dbSendStatement(con, sql)
dbClearResult(q)
}
# pgListRast(con)
if (return_type == 'raster'){
# r <- pgGetRast(con, name = rast_idw)
r <- pgGetRast(
con, name = "z_idw", clauses = glue("WHERE args_hash = '{hash}'")) |>
rast()
return(r)
}

@bbest bbest added this to Management Feb 26, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
Status: No status
Development

No branches or pull requests

1 participant