multilayer vector file data sources with multiple SRS

[brownag]

I noticed this first when other folks were having problems using ESRI File GeoDatabases and vect() that at first I could not reproduce.

Upon closer inspection I found the offending GDBs contained multiple vector layers with different spatial reference systems--one geographic and one projected in case it is relevant.

It appears that there is either an assumption the spatial reference is constant for all layers or only the first SRS is used? I am pretty sure that e.g. GDB and GPKG support multiple SRS. Here is a reproducible example involving a GeoPackage with two layers and two SRS:

library(terra)
#> terra 1.7.19
gdal(lib = "")
#>    gdal    proj    geos 
#> "3.5.2" "8.2.1" "3.9.3"
data(us_states, package = "spData")

x <- vect(us_states)[1,]
y <- project(x, "EPSG:5070")

writeVector(x, "test.gpkg")
writeVector(y, "test.gpkg", layer = "test2", insert = TRUE)

# CRS is as expected
a <- vect("test.gpkg")
#> Warning: [vect] Reading layer: test
#> Other layers: test2
a
#>  class       : SpatVector 
#>  geometry    : polygons 
#>  dimensions  : 1, 6  (geometries, attributes)
#>  extent      : -88.47323, -84.89184, 30.24971, 35.00803  (xmin, xmax, ymin, ymax)
#>  source      : test.gpkg
#>  coord. ref. : lon/lat NAD83 (EPSG:4269) 
#>  names       : GEOID    NAME REGION             AREA total_pop_10 total_pop_15
#>  type        : <chr>   <chr>  <chr>            <chr>        <num>        <num>
#>  values      :    01 Alabama  South 133709.272642372    4.713e+06    4.831e+06
plot(a)

# CRS is wrong (matches first layer)
b <- vect("test.gpkg", layer = "test2")
b
#>  class       : SpatVector 
#>  geometry    : polygons 
#>  dimensions  : 1, 6  (geometries, attributes)
#>  extent      : 704900.9, 1044769, 832261.6, 1376504  (xmin, xmax, ymin, ymax)
#>  source      : test.gpkg (test2)
#>  coord. ref. : lon/lat NAD83 (EPSG:4269) 
#>  names       : GEOID    NAME REGION             AREA total_pop_10 total_pop_15
#>  type        : <chr>   <chr>  <chr>            <chr>        <num>        <num>
#>  values      :    01 Alabama  South 133709.272642372    4.713e+06    4.831e+06
plot(b)
#> Warning: [is.lonlat] coordinates are out of range for lon/lat

# inspect gpkg contents SRS 
con <- DBI::dbConnect(RSQLite::SQLite(), "test.gpkg")
DBI::dbGetQuery(con, "SELECT * FROM gpkg_contents")
#>   table_name data_type identifier description              last_change
#> 1       test  features       test             2023-03-06T21:56:48.897Z
#> 2      test2  features      test2             2023-03-06T21:56:48.942Z
#>          min_x        min_y         max_x        max_y srs_id
#> 1    -88.47323     30.24971     -84.89184 3.500803e+01   4269
#> 2 704900.92430 832261.59123 1044769.37586 1.376504e+06   5070
DBI::dbGetQuery(con, "SELECT * FROM gpkg_spatial_ref_sys")
#>                   srs_name srs_id organization organization_coordsys_id
#> 1  Undefined Cartesian SRS     -1         NONE                       -1
#> 2 Undefined geographic SRS      0         NONE                        0
#> 3                    NAD83   4269         EPSG                     4269
#> 4          WGS 84 geodetic   4326         EPSG                     4326
#> 5     NAD83 / Conus Albers   5070         EPSG                     5070
#>                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         definition
#> 1                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                        undefined
#> 2                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                        undefined
#> 3                                                                                                                                                                                                                                                                                                     GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101],TOWGS84[0,0,0,0,0,0,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4269"]]
#> 4                                                                                                                                                                                                                                                                   GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
#> 5 PROJCS["NAD83 / Conus Albers",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["latitude_of_center",23],PARAMETER["longitude_of_center",-96],PARAMETER["standard_parallel_1",29.5],PARAMETER["standard_parallel_2",45.5],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","5070"]]
#>                                                                description
#> 1                          undefined Cartesian coordinate reference system
#> 2                         undefined geographic coordinate reference system
#> 3                                                                     <NA>
#> 4 longitude/latitude coordinates in decimal degrees on the WGS 84 spheroid
#> 5                                                                     <NA>
DBI::dbDisconnect(con)

[rhijmans]

Thanks! I never considered that. I now get:

vect("test.gpkg", "test") |> crs(proj=TRUE)
#[1] "+proj=longlat +datum=NAD83 +no_defs"
vect("test.gpkg", "test2") |> crs(proj=TRUE)
#[1] "+proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs"

[brownag]

Excellent, thank you Robert!!