2

I have a GeoPackage (.gpkg) database with only one layer:

ogrListLayers("data/LM_data/zipped/ay_riks.gpkg")
[1] "ay_riks"
attr(,"driver")
[1] "GPKG"
attr(,"nlayers")
[1] 1

There are too many objects and I'd like to read only a subset of them based on an attribute:

library(sf)
library(rgdal)
library(tidyverse)    
dta<-st_read("data/LM_data/zipped/ay_riks.gpkg",
                 layer = "ay_riks",
                 EWKB = T,
                 geom_column="geom",
                 query= "SELECT * FROM ay_riks WHERE KOMMUNKOD='2480';")

I get an error:

Reading layer `ay_riks' from data source `/Users/chihayakenji/Documents/spatial_r/data/LM_data/zipped/ay_riks.gpkg' using driver `GPKG'
Error in st_sf(x, ..., agr = agr, sf_column_name = sf_column_name) : 
  no simple features geometry column present

What is the correct way of doing this?

  • I would try to read the GeoPackage first with ogrinfo and verify if "no simple features geometry column present" is misleading error. – user30184 Nov 13 '19 at 16:24
  • There are geometry columns present. If I read the package without a query statement they are there. – Kenji Nov 13 '19 at 16:31
  • 1
    Try geometry_column instead of geom_column: docs: st_read(dsn, layer, ..., query = NA, options = NULL, quiet = FALSE, geometry_column = 1L, type = 0, promote_to_multi = TRUE, although you dont need to specify it unless your table has more than one geom column or isn't a proper geometry table... – Spacedman Nov 13 '19 at 17:13
  • @Spacedman you were spot on. The following worked: dta<-st_read("data/LM_data/zipped/ay_riks.gpkg", layer = "ay_riks", geometry_column="geom", query= "SELECT * FROM ay_riks WHERE KOMMUNKOD = '2480';") – Kenji Nov 13 '19 at 17:29
  • Since this question was essentially due to my lack of attention to the documentation, should I delete it? The answer boils down to "read the documentation carefully." – Kenji Nov 13 '19 at 17:32
2

Although this does not explicitly resolve your issue, the following example should demonstrate that the st_read function can be used successfully in conjunction with an SQL query. As mentioned above, the first step is to confirm that the layer within the GeoPackage contains simple feature geometry column:

ogrinfo lut_data.gpkg 
...
1: lut_data (Multi Polygon)

We can the read in the gpkg into R initially without any SQL query to confirm that all features are read:

library(sf)
library(rgdal)
library(tidyverse)

data <- st_read('lut_data.gpkg', layer='lut_data')

We can confirm that all features have been read in:

dim(data)
[1] 60  3

and plot the output

plot(data)

enter image description here

We can now repeat the above steps, but only import the features where the lut == 'CHF', as follows:

 data <- st_read('lut_data.gpkg', layer='lut_data', \ 
                query='SELECT * FROM lut_data WHERE LUT ="CHF"')

We should now have imported a subset of the features:

dim(data)
[1] 42  3

and we can plot the data to confim the output (and check for missing features from the above figure). enter image description here

  • +1 for the nice guide for doing it. It turns out I was specifying the geometry column argument wrongly and that led to my error. – Kenji Nov 13 '19 at 17:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for?Browse other questions tagged or ask your own question.