Home > Blockchain >  Query a geopackage with a simple features polygon using {sf}
Query a geopackage with a simple features polygon using {sf}

Time:11-02

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!

enter image description here

  • Related