Using duckdb spatial module from R (with sf)?

Tue, Aug 15, 2023 5:34 PM

Hi list,

One more go at this and I'll stop for other ideas.  Below is a
slightly modified example using dbplyr::sql_render() to construct the
query, which works and to me feels a little cleaner, though maybe is
ill-advised.  Then I try to construct a spatial filter (for points in
a polygon) this way with ST_Contains, but my effort comes back empty.
Not sure where I went wrong.  Any ideas?


## boilerplate setup
library(duckdb)
conn <- DBI::dbConnect(duckdb::duckdb())
status <- DBI::dbExecute(conn, "INSTALL 'spatial';")
status <- DBI::dbExecute(conn, "LOAD 'spatial';")
test <- data.frame(site = letters[1:10], latitude = 1:10, longitude = 1:10)
DBI::dbWriteTable(conn, "test", test)

## Here we go, works!
sql <- tbl(con, "test") |>
  mutate(geom = ST_AsHEXWKB(ST_Point(longitude, latitude))) |>
  dbplyr::sql_render()
ex <- st_read(conn, query=sql, geometry_column = "geom", EWKB=FALSE)
ex

## a trivial search polygon to filter:
p2 <- rbind(c(1,1), c(1,2), c(2,2), c(1,1))
pol <-st_polygon(list(p2))
txt <- st_as_text(pol)

## Can we use this to filter duckdb?
sql <- tbl(conn, "test") |>
  mutate(geometry = ST_Point(longitude, latitude),
         geom = ST_AsHEXWKB(geometry)) |>
  filter(ST_Contains(geometry, ST_GeomFromText({txt}))) |>
  dbplyr::sql_render()
sql
ex <- st_read(conn, query=sql, geometry_column = "geom", EWKB=FALSE)
ex
## oh no! our result comes back empty?  did I need CRS on this?  Or do
I missunderstand "ST_Contains()"?

## Here's what the desired result would be from pure sf:
sf <- st_as_sf(test, coords = c(3,2))
filter <- st_as_sf(tibble(sites = "A", geometry = list(pol)))
sf |> st_filter(filter)




---
Carl Boettiger
http://carlboettiger.info/

On Tue, Aug 15, 2023 at 4:42?PM Carl Boettiger <cboettig at gmail.com> wrote:

Ah ha.  if we ask duckdb to coerce it's geometry format to hex, it
appears sf can understand it just fine!

replacing the above query with the following we are good to go.



query <- paste(
  'SELECT *, ST_AsHEXWKB(ST_Point(longitude, latitude)) AS geom',
  'FROM "test"'
)
ex <- st_read(con, query=query, geometry_column = "geom", EWKB=FALSE)
ex


I'll experiment further.  I'm terrible at SQL, and to my eyes it looks
needlessly verbose.  I'm keen to understand how I can better leverage
sf's notation to compose the sql queries to duckdb....  but it seems
to work!   I'm also still trying to determine if duckdb is using EWKB
or vanilla WKB....

---
Carl Boettiger
http://carlboettiger.info/

On Tue, Aug 15, 2023 at 4:02?PM Carl Boettiger <cboettig at gmail.com> wrote:

Hi folks,


I'm curious if anyone has explored the relatively new spatial
extension in duckdb (https://duckdb.org/docs/extensions/spatial.html)
or has any pointers/tutorials regarding its use from R?

Consider the following minimal example that seeks to use the sf
library to speak to duckdb:

  library(duckdb)
  library(sf)
  conn <- DBI::dbConnect(duckdb::duckdb())
  status <- DBI::dbExecute(conn, "INSTALL 'spatial';")
  status <- DBI::dbExecute(conn, "LOAD 'spatial';")

  test <- data.frame(site = letters[1:10], latitude = 1:10, longitude = 1:10)
  DBI::dbWriteTable(conn, "test", test)

# Ok, let's try and make a geometry column
  query <- paste(
    'SELECT *, ST_Point(longitude, latitude) AS geom',
    'FROM "test"'
  )

  ex <- st_read(con, query=query, geometry_column = "geom")
  ## error: reading wkb type 0 is not supported


  ex <- st_read(con, query=query, geometry_column = "geom", EWKB = FALSE)
  ## prints: wkbType: 1572864
  ## Error in CPL_read_wkb(x, EWKB, spatialite) : unsupported wkbType
dim in switch

 We seem to get closer than I might have hoped (sf doesn't just insist
that conn isn't postgresgis), but is getting stuck on the encoding of
the wkb.  Is this something we can work around?  The queries seem to
run successfully, I just seem to be getting some unsupported ecoding
of the WKB geometry column....

---
Carl Boettiger
http://carlboettiger.info/

Using duckdb spatial module from R (with sf)?

Thread (8 messages)