Add dialect argument for ExecuteSQL

[agila5]

Dear sf developers, I created this PR to propose the inclusion of a dialect argument into st_read(). The new argument can be used to choose the dialect used by ExecuteSQL. The idea for including a new parameter into st_read started here. Still, I think that the possibility of selecting the SQLite dialect can be useful in other contexts. Hence, I created some examples (taken from here) to present the new functionalities:

# packages
library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1

# tests and examples taken from https://gdal.org/user/sql_sqlite_dialect.html
nc_path <- system.file("shape/nc.shp", package = "sf")

# ST_Area
nc <- st_read(nc_path, query = "SELECT GEOMETRY, ST_Area(GEOMETRY) AS AREA FROM nc", dialect = "SQLite")
#> Reading layer `nc' from data source 
#>   `C:\Users\Utente\Documents\R\win-library\4.0\sf\shape\nc.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 100 features and 1 field
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> Geodetic CRS:  NAD27
nc
#> Simple feature collection with 100 features and 1 field
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> Geodetic CRS:  NAD27
#> First 10 features:
#>          AREA                       GEOMETRY
#> 1  0.11428350 MULTIPOLYGON (((-81.47276 3...
#> 2  0.06139976 MULTIPOLYGON (((-81.23989 3...
#> 3  0.14301628 MULTIPOLYGON (((-80.45634 3...
#> 4  0.06977098 MULTIPOLYGON (((-76.00897 3...
#> 5  0.15275930 MULTIPOLYGON (((-77.21767 3...
#> 6  0.09715756 MULTIPOLYGON (((-76.74506 3...
#> 7  0.06188045 MULTIPOLYGON (((-76.00897 3...
#> 8  0.09080191 MULTIPOLYGON (((-76.56251 3...
#> 9  0.11846603 MULTIPOLYGON (((-78.30876 3...
#> 10 0.12385371 MULTIPOLYGON (((-80.02567 3...

# ST_Area + ST_Intersects considering a circle with radius 1m centered in the centroid of North Carolina
st_read(
  dsn = nc_path,
  query = "
  SELECT GEOMETRY, ST_Area(GEOMETRY) AS AREA
  FROM nc
  WHERE ST_Intersects(ST_Transform(GEOMETRY, 32119), BuildCircleMbr(573193.5, 199429.3, 1, 32119))
  ",
  dialect = "SQLite"
)
#> Reading layer `nc' from data source 
#>   `C:\Users\Utente\Documents\R\win-library\4.0\sf\shape\nc.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 1 feature and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -79.55536 ymin: 35.51024 xmax: -78.90572 ymax: 35.86815
#> Geodetic CRS:  NAD27

# Same example with radius 100km
st_read(
  dsn = nc_path,
  query = "
  SELECT GEOMETRY, ST_Area(GEOMETRY) AS AREA
  FROM nc
  WHERE ST_Intersects(ST_Transform(GEOMETRY, 32119), BuildCircleMbr(573193.5, 199429.3, 100000, 32119))
  ",
  dialect = "SQLite"
)
#> Reading layer `nc' from data source 
#>   `C:\Users\Utente\Documents\R\win-library\4.0\sf\shape\nc.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 38 features and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -80.96577 ymin: 34.30457 xmax: -77.80411 ymax: 36.56521
#> Geodetic CRS:  NAD27

# ogr_layer_Extent (st_bbox) + ogr_layer_SRID(CRS) + ogr_layer_GeometryType(st_geometry_type) + ogr_layer_FeatureCount(nrow?)
nc_extent <- st_read(
  dsn = nc_path,
  query = "SELECT 
    ogr_layer_Extent('nc') AS extent, 
    ogr_layer_SRID('nc') AS SRID, 
    ogr_layer_GeometryType('nc') AS GeometryType, 
    ogr_layer_FeatureCount('nc') AS count
  ",
  dialect = "SQLite"
)
#> Reading layer `nc' from data source 
#>   `C:\Users\Utente\Documents\R\win-library\4.0\sf\shape\nc.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 1 feature and 3 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> Geodetic CRS:  NAD27
nc_extent
#> Simple feature collection with 1 feature and 3 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> Geodetic CRS:  NAD27
#>   SRID GeometryType count                         extent
#> 1 4267      POLYGON   100 POLYGON ((-84.32385 33.8819...

st_equals(st_geometry(nc_extent), st_as_sfc(st_bbox(nc)))
#> Sparse geometry binary predicate list of length 1, where the predicate
#> was `equals'
#>  1: 1

# OGR Geocoding
nc <- st_read(
  dsn = nc_path,
  query = "
    SELECT GEOMETRY, ogr_geocode_reverse(
      ST_Centroid(GEOMETRY), 
      'county'
      ) AS x 
    FROM nc 
    LIMIT 5
  ",
  dialect = "SQLite"
)
#> Reading layer `nc' from data source 
#>   `C:\Users\Utente\Documents\R\win-library\4.0\sf\shape\nc.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 5 features and 1 field
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -81.74107 ymin: 36.07282 xmax: -75.77316 ymax: 36.58965
#> Geodetic CRS:  NAD27
nc
#> Simple feature collection with 5 features and 1 field
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -81.74107 ymin: 36.07282 xmax: -75.77316 ymax: 36.58965
#> Geodetic CRS:  NAD27
#>                    x                       GEOMETRY
#> 1        Ashe County MULTIPOLYGON (((-81.47276 3...
#> 2   Alleghany County MULTIPOLYGON (((-81.23989 3...
#> 3       Surry County MULTIPOLYGON (((-80.45634 3...
#> 4   Currituck County MULTIPOLYGON (((-76.00897 3...
#> 5 Northampton County MULTIPOLYGON (((-77.21767 3...

^{Created on 2021-04-15 by the reprex package (v2.0.0)}

I have minimal knowledge of C++/GDAL/GEOS/SQLite, so I'm not sure if I should add more tests, additional checks or implement a completely different approach. Anyway, if you agree, I could one example in ?st_read.

[edzer]

LGTM! Would you mind adding a minimal test for this, somewhere?

[agila5]

Done!

[edzer]

Thanks, Andrea!

[edzer]

See #1663 - I'm reverting this, as it doesn't seem to work on GDAL installs without spatialite.