Nothing
## pgGetRast
##' Load raster from PostGIS database into R.
##'
##' Retrieve rasters from a PostGIS table into a \code{terra SpatRaster} object
##'
##'
##' Since version 1.5.0, this function retrieve SpatRaster objects from
##' \code{terra} package by default. The argument \code{returnclass} can be
##' used to return \code{raster} objects instead.
##'
##' The argument \code{bands} can take as argument:
##'
##' * The index of the desirable band (e.g. bands = 2 will fetch the second band
##' of the raster).
##'
##' * More than one index for several bands (e.g. bands = c(2,4) will return a
##' \code{SpatRaster} with two bands).
##'
##' * All bands in the raster (bands = TRUE).
##'
##'
##' @param conn A connection object to a PostgreSQL database
##' @param name A character string specifying a PostgreSQL schema and
##' table/view name holding the geometry (e.g., \code{name =
##' c("schema","table")})
##' @param rast Name of the column in \code{name} holding the raster object.
##' Defaults to "rast".
##' @param clauses character, optional SQL to append to modify select
##' query from table. Must begin with 'WHERE'.
##' @param bands Index number(s) for the band(s) to retrieve (defaults to 1).
##' The special case (\code{bands = TRUE}) returns all bands in the raster. See
##' also 'Details'
##' @param boundary \code{sf} object, \code{SpatVector} object, or numeric. If a
##' spatial object is provided, its bounding box will be used to select
##' the part of the raster to import. Alternatively, a numeric vector
##' (\code{c([top], [bottom], [right], [left])}) indicating the
##' projection-specific limits with which to clip the raster. If not value
##' is provided, the default \code{boundary = NULL} will return the
##' full raster.
##' @param returnclass 'terra' by default; or 'raster' for \code{raster} objects.
##' @param progress whether to show a progress bar (TRUE by default). The progress
##' bar mark the progress of reading bands from the database.
##'
##' @author David Bucklin \email{david.bucklin@@gmail.com} and Adrián Cidre
##' González \email{adrian.cidre@@gmail.com}
##' @importFrom terra rast ext crs crop
##' @importFrom sf st_crs st_bbox
##' @importFrom raster raster stack
##' @export
##' @return \code{SpatRaster}; \code{raster}; or \code{RasterStack} object
##' @examples
##' \dontrun{
##' pgGetRast(conn, c("schema", "tablename"))
##' pgGetRast(conn, c("schema", "DEM"), boundary = c(55,
##' 50, 17, 12))
##' }
pgGetRast <- function(conn, name, rast = "rast", bands = 1,
boundary = NULL, clauses = NULL,
returnclass = "terra", progress = TRUE) {
## Message
message("Since version 1.5 this function outputs SpatRaster objects by default. Use returnclass = 'raster' to return raster objects.")
## Check connection and PostGIS extension
dbConnCheck(conn)
if (!suppressMessages(pgPostGIS(conn))) {
stop("PostGIS is not enabled on this database.")
}
## Check and prepare the schema.name
name1 <- dbTableNameFix(conn, name)
nameque <- paste(name1, collapse = ".")
namechar <- gsub("'", "''", paste(gsub('^"|"$', '', name1), collapse = "."))
## Raster query name
rastque <- dbQuoteIdentifier(conn, rast)
## Fix user clauses
clauses2 <- sub("^where", "AND", clauses, ignore.case = TRUE)
## Check table exists and return error if it does not exist
tmp.query <- paste0("SELECT r_raster_column AS geo FROM raster_columns\n WHERE (r_table_schema||'.'||r_table_name) = '",
namechar, "';")
tab.list <- dbGetQuery(conn, tmp.query)$geo
if (is.null(tab.list)) {
stop(paste0("Table '", namechar, "' is not listed in raster_columns."))
} else if (!rast %in% tab.list) {
stop(paste0("Table '", namechar, "' raster column '", rast,
"' not found. Available raster columns: ", paste(tab.list,
collapse = ", ")))
}
## Check bands
tmp.query <- paste0("SELECT st_numbands(", rastque, ") FROM ", nameque, " WHERE ", rastque, " IS NOT NULL LIMIT 1;")
nbs <- 1:dbGetQuery(conn, tmp.query)[1,1]
if (isTRUE(bands)) {
bands <- nbs
} else if (!all(bands %in% nbs)) {
stop(paste0("Selected band(s) do not exist in PostGIS raster: choose band numbers between ", min(nbs), " and ", max(nbs), "."))
}
## Retrieve the SRID
tmp.query <- paste0("SELECT DISTINCT(ST_SRID(", rastque, ")) FROM ", nameque, " WHERE ", rastque, " IS NOT NULL;")
srid <- dbGetQuery(conn, tmp.query)
## Check if the SRID is unique, otherwise throw an error
if (nrow(srid) > 1) {
stop("Multiple SRIDs in the raster")
} else if (nrow(srid) < 1) {
stop("Database table is empty.")
}
## Retrieve the proj4string
p4s <- NA
# tmp.query <- paste0("SELECT proj4text AS p4s FROM spatial_ref_sys WHERE srid = ",
# srid$st_srid, ";")
tmp.query.sr <- paste0("SELECT r_proj4 AS p4s FROM ", nameque, ";")
try(db.proj4 <- dbGetQuery(conn, tmp.query.sr)$p4s, silent = TRUE)
# if db.proj4 doesnt exist (error because raster was not loaded using rpostgis)
if (!exists("db.proj4")) {
tmp.query.sr <- paste0("SELECT proj4text AS p4s FROM spatial_ref_sys WHERE srid = ",
srid$st_srid, ";")
db.proj4 <- dbGetQuery(conn, tmp.query.sr)$p4s
}
if (!is.null(db.proj4)) {
try(p4s <- terra::crs(db.proj4[1]), silent = TRUE)
}
if (is.na(p4s)) {
warning("Table SRID not found. Projection will be undefined (NA)")
}
## Check alignment of raster
tmp.query <- paste0("SELECT ST_SameAlignment(", rastque, ") FROM ", nameque, ";")
# needs postgis version 2.1+, so just try
al <- FALSE
try(al <- dbGetQuery(conn, tmp.query)[1,1])
if (!al) {
# get alignment from upper left pixel of all raster tiles
tmp.query <- paste0("SELECT min(ST_UpperLeftX(", rastque, ")) ux, max(ST_UpperLeftY(", rastque, ")) uy FROM ", nameque, ";")
aligner <- dbGetQuery(conn, tmp.query)
aq <- c("ST_SnapToGrid(", paste0(aligner[1,1],","), paste0(aligner[1,2],"),"))
} else {
aq <- NULL
}
## Get band function
get_band <- function(band) {
## Get raster information (bbox, rows, cols)
info <- dbGetQuery(conn, paste0("select
st_xmax(st_envelope(rast)) as xmax,
st_xmin(st_envelope(rast)) as xmin,
st_ymax(st_envelope(rast)) as ymax,
st_ymin(st_envelope(rast)) as ymin,
st_width(rast) as cols,
st_height(rast) as rows
from
(select st_union(",aq[1],rastque,",",aq[2],aq[3],band,") rast from ",nameque," ", clauses,") as a;"))
## Retrieve values of the cells
vals <- dbGetQuery(conn, paste0("select
unnest(st_dumpvalues(rast, 1)) as vals
from
(select st_union(",aq[1],rastque,",",aq[2],aq[3],band,") rast from ",nameque," ", clauses,") as a;"))$vals
rout <- terra::rast(nrows = info$rows, ncols = info$cols, xmin = info$xmin,
xmax = info$xmax, ymin = info$ymin, ymax = info$ymax,
crs = p4s, vals = vals)
return(rout)
}
## Get band with boundary function
get_band_boundary <- function(band) {
## Get info
info <- dbGetQuery(conn, paste0("select
st_xmax(st_envelope(rast)) as xmx,
st_xmin(st_envelope(rast)) as xmn,
st_ymax(st_envelope(rast)) as ymx,
st_ymin(st_envelope(rast)) as ymn,
st_width(rast) as cols,
st_height(rast) as rows
from
(select st_union(",aq[1],rastque,",",aq[2],aq[3],band,") rast from ",nameque, "\n
WHERE ST_Intersects(",
rastque, ",ST_SetSRID(ST_GeomFromText('POLYGON((", boundary[4],
" ", boundary[1], ",", boundary[4], " ", boundary[2],
",\n ", boundary[3], " ", boundary[2], ",", boundary[3],
" ", boundary[1], ",", boundary[4], " ", boundary[1],
"))'),", srid, "))", clauses2,") as a;"))
if (is.na(info$cols) & is.na(info$rows)) {
stop("No data found within geographic subset defined by 'boundary'.")
}
vals <- dbGetQuery(conn,paste0("select
unnest(st_dumpvalues(rast, 1)) as vals
from
(select st_union(",aq[1],rastque,",",aq[2],aq[3],band,") rast from ",nameque, "\n
WHERE ST_Intersects(",
rastque, ",ST_SetSRID(ST_GeomFromText('POLYGON((", boundary[4],
" ", boundary[1], ",", boundary[4], " ", boundary[2],
",\n ", boundary[3], " ", boundary[2], ",", boundary[3],
" ", boundary[1], ",", boundary[4], " ", boundary[1],
"))'),", srid, "))", clauses2,") as a;"))$vals
rout <- terra::rast(nrows = info$rows, ncols = info$cols,
xmin = info$xmn, xmax = info$xmx, ymin = info$ymn, ymax = info$ymx,
crs = p4s, vals = vals)
return(rout)
}
## Get raster
if (is.null(boundary)) {
## Get bands
if (progress) {
rout <- purrr::map(bands, get_band, .progress = "Reading bands")
} else {
rout <- purrr::map(bands, get_band)
}
rb <- terra::rast(rout)
## Else: when boundary is provided
} else {
## Bbox of terra and sf objects
if (inherits(boundary, "sf")) {
boundary <- sf::st_bbox(boundary)
boundary <- c(boundary[4], boundary[2], boundary[3], boundary[1])
} else if (inherits(boundary, "SpatVector")) {
boundary <- c(terra::ext(boundary)[4], terra::ext(boundary)[3],
terra::ext(boundary)[2], terra::ext(boundary)[1])
}
## Extent to clip the Rast
extclip <- terra::ext(boundary[4], boundary[3], boundary[2], boundary[1])
## Get bands
rout <- purrr::map(bands, get_band_boundary, .progress = "Reading bands")
rb <- terra::rast(rout)
}
## Set layer names
if ("band_names" %in% dbTableInfo(conn,name)$column_name) {
try({
ct <- 1
for (b in bands) {
lnm <- dbGetQuery(conn, paste0("SELECT DISTINCT band_names[",b,
"][1] as nm FROM ",nameque," ", clauses,";"))
names(rb)[ct] <- lnm$nm
ct <- ct + 1
}
})
}
# precise cropping
if (!is.null(boundary)) {
rb_final <- terra::crop(rb, extclip)
} else {
rb_final <- rb
}
# Output terra or raster
if (returnclass == "terra") {
return(rb_final)
} else if (returnclass == "raster") {
if (terra::nlyr(rb_final) == 1) {
return(raster::raster(rb_final))
} else {
return(raster::stack(rb_final))
}
} else {
stop("returnclass must be one of 'terra' or 'raster'")
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.