Error reading spatialite table with sf::st_read

Loic,

I agree with Barry that there is a (logic) bug there, leading to a 
confusing error message. The problem seems the DBI driver failing to 
interpret a binary column ("blob") as wkb, which is needed when taking 
this path:

tb = dbReadTable(con, "olinda1")
head(tb)

? ogc_fid??? id????? cd_geocodi?? tipo?? cd_geocodb??? nm_bair v014 
GEOMETRY
1?????? 1 28801 260960005000001 URBANO 260960005020 Ouro Preto 1119 
blob[484 B]
2?????? 2 28802 260960005000002 URBANO 260960005020 Ouro Preto 1267 
blob[532 B]
3?????? 3 28803 260960005000003 URBANO 260960005020 Ouro Preto 557 
blob[276 B]
4?????? 4 28804 260960005000004 URBANO 260960005020 Ouro Preto 709 
blob[468 B]
5?????? 5 28805 260960005000005 URBANO 260960005020 Ouro Preto 1045 
blob[452 B]
6?????? 6 28806 260960005000006 URBANO 260960005020 Ouro Preto 727 
blob[388 B]

st_as_sfc(tb$GEOMETRY)

wkbType: 30078980
Error in CPL_read_wkb(x, EWKB, spatialite) :
? unsupported wkbType dim in switch

so there seems to be something special with spatialite's WKB, which is 
not recognized by sf's WKB reader. Reading ?st_as_sfc:

st_as_sfc(tb$GEOMETRY, spatialite = TRUE)

Geometry set for 470 features
Geometry type: POLYGON
Dimension:???? XY
Bounding box:? xmin: -34.91692 ymin: -8.044467 xmax: -34.82779 ymax: 
-7.954672
CRS:?????????? NA
First 5 geometries:
POLYGON ((-34.86406 -7.992488, -34.86397 -7.992...
POLYGON ((-34.86129 -7.989947, -34.86087 -7.989...
POLYGON ((-34.85856 -7.992407, -34.85899 -7.992...
POLYGON ((-34.85752 -7.994827, -34.85722 -7.994...
POLYGON ((-34.8582 -7.997543, -34.858 -7.996937...
Warning message:
In CPL_crs_from_input(x) :
? GDAL Error 1: PROJ: proj_create_from_database: crs not found

This can no doubt be made more automatic (e.g. by adding a spatialite 
argument to st_read.DBIObject), but the the question remains why, when

x = st_read('/tmp/olinda1.sqlite')

works OK, that couldn't be used instead?

Many regards,


On 03/03/2022 19:19, Barry Rowlingson wrote:

Looks like that error might be generated about here:
https://github.com/r-spatial/sf/blob/5dde39c8261dc6b202e6cde9d41ad3bf0f46aa3a/R/db.R#L134 


where the code does:

try(sfc <- st_as_sfc(tbl[[i]]), silent = TRUE)
if (!inherits(sfc, "try-error")) {
tbl[[i]] = sfc
success = TRUE
}

If the first "try" fails then `sfc` doesn't get created and so the error
happens on `if(!inherits(sfc,....)`.

The code is looping over columns trying to find a valid geometry column,
but not trapping this error properly. A bug, I think. Feel free to 
report
as an issue.

However I don't think that gets us out of the woods - that code only 
talks
about PostGIS databases, but it gets called for any object of class
DBIObject, which your SQLite connection object is:

? > inherits(con,"DBIObject")
? [1] TRUE

Perhaps that method should be something like st_read.PGConnection if it
only works on PostGIS databases? Experimenting with bits of that 
function
don't get me the geometry from my test spatialite DB anyway, and I think
its due to it being PostGIS-specific. I kept getting empty geometries.

But why not read the Spatialite directly without going through DBI?? 
Here's
a test I made of some points:

d = st_read("/tmp/olinda1.sqlite","pts")

Reading layer `pts' from data source `/tmp/olinda1.sqlite' using driver
`SQLite'
Simple feature collection with 1000 features and 1 field
Geometry type: POINT
Dimension:???? XY
Bounding box:? xmin: 234250 ymin: -198750 xmax: 909750 ymax: 502750
CRS:?????????? unknown

THe CRS is printed as unknown but the WKT is there. Maybe this is a 
funny
CRS I was using for testing or something:

st_crs(d)

Coordinate Reference System:
?? User input: unknown
?? wkt:
PROJCRS["unknown",
???? BASEGEOGCRS["GCS_unknown",
DATUM["D_Unknown_based_on_Australian_Natl_S_Amer_1969_ellipsoid",
ELLIPSOID["Australian_Natl_S_Amer_1969",6378160,298.25,
???????????????? LENGTHUNIT["metre",1,

and you can add SQL queries to `st_read` if that;s part of your 
motivation
for using a spatial file database:

d = st_read("/tmp/olinda1_no.sqlite", query="select * from pts where

pampa_srtm>300")
Reading query `select * from pts where pampa_srtm>300' from data source
`/tmp/olinda1_no.sqlite' using driver `SQLite'
Simple feature collection with 189 features and 1 field

Barry

On Thu, Mar 3, 2022 at 4:48 PM DUTRIEUX Loic 
<Loic.DUTRIEUX at ec.europa.eu>
wrote:

Hi everyone,

When trying to read data contained in a sqlite table with the full
spatialite options (see what FORMAT=SPATIALITE implies here -->
https://gdal.org/drivers/vector/sqlite.html), I get a
wkbType: 30078980 Error in inherits(sfc, "try-error") : object 'sfc' 
not
found.

See the example below:

# In bash
ogr2ogr -F SQLITE -dsco SPATIALITE=YES /tmp/olinda1.sqlite
~/R/x86_64-pc-linux-gnu-library/4.0/sf/shape/olinda1.shp

# In R
library(sf)
library(DBI)

con <- dbConnect(RSQLite::SQLite(), '/tmp/olinda1.sqlite')
st_read(con, query='select * from olinda1;')


If I create the database without the SPATIALITE=YES creation option
(geometries encoded as WKB), reading with sf partly works (no CRS) 
but then
I lose the spatialite functionalities. Any suggestion?

Kind regards,
Lo?c

????[[alternative HTML version deleted]]