I have a large geopackage. I want to read in features from it, that intersect with another geopackage. I am trying this but I'm obviously doing something wrong:
box <- st_read("file1.gpkg", quiet=T) %>% st_bbox()
layer_name <- st_layers("file2.gpkg")$name
my_query <- glue("SELECT * FROM {layer_name} WHERE st_intersects(geom, st_polygon('{box$xmin} {box$ymin}, {box$xmax} {box$ymin}, {box$xmax} {box$ymax}, {box$xmin} {box$ymax}, {box$xmin} {box$ymin}'))")
st_read("file2.gpkg", query = my_query)
If I print the query it looks like this, is it formatted correctly?
SELECT * FROM li_final_town_areas WHERE st_intersects(geom, st_polygon('9.25 47, 9.5 47, 9.5 47.25, 9.25 47.25, 9.25 47'))
The error that I get is:
Error in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
Query execution failed, cannot open layer.
CodePudding user response:
This form of the query, using st_polygonfromtext
with the polygon as WKT, works for me:
q2 = "SELECT * FROM gadm36_BTN_2 WHERE st_intersects(geom, st_polygonfromtext('POLYGON(( 90 27, 91 27, 91 28, 90 28 , 90 27 ))' ))"
BTN2 = st_read(gpkg, query=q2)
Here's the selected level-2 regions of Bhutan inside level-0 Bhutan with the selection polygon lats and longs added. Seems to work!