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 <cboet...@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 <cboet...@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 <cboet...@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@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo