st_read and requiring layer name

[rundel]

This is a small pet peeve of mine with regard to rgdal - is there a compelling reason to continue to force the user to provide the layer name when reading in data?

There isn't a technical reason within gdal to require this and in my experience most geometry files I've worked with only ever have one layer. It seems to me a more sensible approach is to make the layer name an optional argument and in the case of a single layer to just automatically return that object.

If more than one layer are present then the current error could be thrown, or alternatively throw a warning and return a list containing all layers.

[edzer]

Nice idea, would improve usability! What would be a good approach to try to derive layer name from a single string: remove the file name extension and possibly path from the file name?

[rundel]

I'm not sure, my use cases have almost always been geojson and shapefiles so I don't have direct experience with the more exotic vector types and their usage of layer names.

GDAL peculiarities mean that both of the following work without issue

nc <- st_read(system.file("shape/nc.shp", package="sf"), "nc", crs = 4267)
nc <- st_read(system.file("shape/", package="sf"), "nc", crs = 4267)

so I'm not sure if just looking at the path will be helpful.

The example I have the most issue with has always been with reading geojson where the user is expected to know that the layer is called OGRGeoJSON - that has always been a pain point when my students are first learning this stuff.

[Robinlovelace]

I agree with the changes suggested by @rundel. That will increase usability lots. It was a godsend when I found out about raster::shapefile() and I think st_read() would benefit from working in the same or a similar way.

I think it's easiest to just point it at the file automatically deciding the driver/layer based on the extension. This is what rio does and I think sf should too: https://cran.r-project.org/web/packages/rio/vignettes/rio.html

[edzer]

So these might be some heuristics then: if st_read is called with a single argument, it must be an existing file, and

1. if the file name ends in .shp, we assume it is a shapefile (with full path), and call the layer the basename (i.e. extension and path stripped) of the file name
2. Geopackages end in .gpkg (see req 3, so use the same trick as for shapefiles;
3. Other file types with fixed extension insert here
4. other file names, we assume will be GeoJSON, and layer will be called OGRGeoJSON (see here)

[Robinlovelace]

Some other common extensions that I think would be good for sf to recognise:

.geojson
.osm
.gml
.kml
.kmz

I think raster already does a good job of handling the multiple raster formats.

[edzer]

Examples of using st_read on .geojson .osm, .gml, .kml and .kmz would be helpful!

[rsbivand]

Please be careful to separate out st_read for dummies and the real st_read, which shouldn't make assumptions. A typical problem is with the shapefile driver which falls over in a dsn=directory with a given layer=, but where the directory contains a *.dbf without being part of a shapefile. Robert's use of shapefile as a general term for vector data is very unhelpful. The border between sensible simplification and going too far is very real. Something not (yet) in the frame is encoding, which is often very unpleasant (even CP1252 is a problem), and far from everyone is UTF-8 yet. So either do a wrapper for st_read which handles frequent cases, and explain that more often than we like, real data is messy and needs the underlying complications in the bare interface, or just provide the bare interface and enough guidance. Fixing GeoJSON for the current odd OGR naming may change as soon as the driver gets updated - there are reasons for exposing users to the actual complications, because it makes derailments their responsibility.

[Robinlovelace]

Example of a 5 Mb file I created - cycling potential on major roads in London:

u = "https://github.com/npct/pct-data/raw/master/london/rnet.geojson"
download.file(url = u, destfile = "/tmp/rnet.geojson")
rnet = sf::st_read("/tmp/rnet.geojson") # fails
rnet = sf::st_read(dsn = "/tmp/rnet.geojson", layer = "OGRGeoJSON") # works

I take Roger's point about precision and robustness vs ease of use and like the idea of having a separate wrapper function for magically guessing, most of the time correctly, how to read in a file. For that I would propose something with a name like st_import(), akin to rio::import().

Would that satisfy all parties?

rnet <- sf::st_import("/tmp/rnet.geojson") 

Would work for me and think it's a function I can even help writing!

[edzer]

I will not provide two functions. We can give layer a default name where this makes sense, and an error message when it is missing otherwise (e.g. where dsn is a directory, or an extension is missing). Like rgdal, the function responds with printing which layer is being read from which data source, on success, so that is there. Even ogr2ogr doesn't force you to specify layers if they're obvious.

[Robinlovelace]

Sounds reasonable and can understand desire to minimise number of functions. So, to clarify, that would mean st_read("file.geojson") and st_read("file.shp") would work? I think that would make life easier for lots of people, myself included, especially if I end up teaching sf which I hope will happen based on the impressive rate of development.

[edzer]

no, because that doesn't point to a file with an extension that tells us what the layer name should be. I just pushed a commit that accepts

 st_read(system.file("shape/nc.shp", package="sf"))

though -- please test (now supports gpkg, geojson, shp).

To be discussed: what will it return as default?

[Robinlovelace]

• here's an example xml (.osm) file. Not sure if this can be read with rgdal or sf at present. Btw - fast work Edzer! Will aim to test during the week, maybe Weds.
  doc_txt.osm.zip

[Robinlovelace]

Could you clarify what you mean by "that doesn't point to a file with an extension that tells us what the layer name should be" - not understanding 100%.

[edzer]

I meant: literally the command you typed. It has no file extension since the "/" comes after the ".". If you mean, either

st_read(file.geojson")

or

st_read("file.shp")

then yes, but please don't imply meaning in github issues -- I can only read what you write.

Which layer in your .osm? ogrinfo reports 5:

edzer@gin-edzer:/tmp$ ogrinfo doc_txt.osm
Had to open data source read-only.
INFO: Open of `doc_txt.osm'
      using driver `OSM' successful.
1: points (Point)
2: lines (Line String)
3: multilinestrings (Multi Line String)
4: multipolygons (Multi Polygon)
5: other_relations (Geometry Collection)

[edzer]

library(sf)
o1 = st_read("doc_txt.osm", "points")
o2 = st_read("doc_txt.osm", "lines")
o3 = st_read("doc_txt.osm", "multilinestrings")
o4 = st_read("doc_txt.osm", "multipolygons")
o5 = st_read("doc_txt.osm", "other_relations")

plot(o1)
plot(o2, add = TRUE, col = 'red')
plot(o3, add = TRUE, col = 'blue')
plot(o4, add = TRUE, col = 'green')
plot(o5, add = TRUE, col = 'orange')

(the last two are empty)