[R-sig-Geo] Error reading spatialite table with sf::st_read

Sat Mar 5 16:52:45 CET 2022


On 05/03/2022 16:29, Micha Silver wrote:
> 
> On 04/03/2022 13:36, Edzer Pebesma wrote:
>> 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:
> 
> 
> I'm certainly no expert in WKB, but I did notice on the spatialite group 
> list the following thread:
> 
> https://groups.google.com/g/spatialite-users/c/fJSlg6qGUfc
> 
> where Sandro explains that the spatialite "geom" column is NOT in 
> standard WKB...

thanks - that's what sf takes care of when setting spatialite = TRUE in 
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 using 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
>>>>
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> R-sig-Geo using r-project.org
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>
>>>     [[alternative HTML version deleted]]
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo using r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>

-- 
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081