[R-sig-Geo] Using duckdb spatial module from R (with sf)?

[Josiah Parry]

Quick note on the contains: a polygon has 5 points! The first has to be the
same as the last. And they shuold be going in Counter clock wise order if
memory serves ! :)

On Tue, Aug 15, 2023 at 8:34 PM Carl Boettiger <cboettig using gmail.com> wrote:

> 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 using 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 using 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/
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

	[[alternative HTML version deleted]]