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

DUTRIEUX Loic Lo|c@DUTRIEUX @end|ng |rom ec@europ@@eu
Fri Mar 4 13:13:50 CET 2022


Thanks Edzer and Barry, this helps a lot.
For some reason I wrongly assumed that dsn= had to be a database connection for the query= argument to be valid.

The route st_read('/tmp/olinda1.sqlite', query=...) works perfectly fine, also for more complex query involving joins and multiple tables, therefore I will use that.
Happy to help if you consider that these features (reading spatialite via database connection, recognizing/parsing spatialite binary format, ...) deserves some attention.

Kind regards,
Loïc
________________________________________
From: R-sig-Geo <r-sig-geo-bounces using r-project.org> on behalf of Edzer Pebesma <edzer.pebesma using uni-muenster.de>
Sent: 04 March 2022 12:36
To: r-sig-geo using r-project.org
Subject: Re: [R-sig-Geo] 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://urldefense.com/v3/__https://github.com/r-spatial/sf/blob/5dde39c8261dc6b202e6cde9d41ad3bf0f46aa3a/R/db.R*L134__;Iw!!DOxrgLBm!UziVBf_G2dc6TzzBAph_Chp0X8oSDIkLNwS6KCdSNNZOOjj65vqiG_SRaxPztyFqCuKrKfQ$
>
> 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://urldefense.com/v3/__https://gdal.org/drivers/vector/sqlite.html__;!!DOxrgLBm!UziVBf_G2dc6TzzBAph_Chp0X8oSDIkLNwS6KCdSNNZOOjj65vqiG_SRaxPztyFq1ixAQg8$ ), 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://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-geo__;!!DOxrgLBm!UziVBf_G2dc6TzzBAph_Chp0X8oSDIkLNwS6KCdSNNZOOjj65vqiG_SRaxPztyFqp6XITJk$
>>
>
>       [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-geo__;!!DOxrgLBm!UziVBf_G2dc6TzzBAph_Chp0X8oSDIkLNwS6KCdSNNZOOjj65vqiG_SRaxPztyFqp6XITJk$

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

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo using r-project.org
https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-geo__;!!DOxrgLBm!UziVBf_G2dc6TzzBAph_Chp0X8oSDIkLNwS6KCdSNNZOOjj65vqiG_SRaxPztyFqp6XITJk$



More information about the R-sig-Geo mailing list