st_read() errors with GML StringList field

[jlehtoma]

I'm trying to read GML data using st_read(), which is causing an error.

library(sf)
library(rgdal)

# Define a query URL to Finnish Meteorological Institute WFS. The response
# contains data from 22 weather stations in Southern Finland within a single day.
# Data is not tidy. The URL below uses a personal API-key I'll delete later on.
> wfs_url <- 'http://data.fmi.fi/fmi-apikey/0253711c-056f-439b-af9d-a5a3a2920faf/wfs?request=getFeature&storedquery_id=fmi::observations::weather::daily::timevaluepair&starttime=2014-01-01&endtime=2014-01-01&bbox=23.5,59.7,26.5,60.7&parameters=tday'

# Download as a local file
> fmi_gml_file <- "fmi_test.gml"
> download.file(wfs_url, fmi_gml_file)

Reading fails with st_read():

> sf_dat <- sf::read_sf(fmi_gml_file)
Error in CPL_read_ogr(dsn, layer, as.character(options), quiet, type,  : 
  Unrecognized field type

Reading with rgdal::readOGR() works:

> ogr_dat <- rgdal::readOGR(fmi_gml_file)
OGR data source with driver: GML 
Source: "fmi_test.gml", layer: "PointTimeSeriesObservation"
with 22 features
It has 10 fields, of which 1 list fields

The 1 list field could be the culprit:

> rgdal::ogrInfo(fmi_gml_file)
Source: "fmi_test.gml", layer: "PointTimeSeriesObservation"
Driver: GML; number of rows: 22 
Feature type: wkbPoint with 2 dimensions
Extent: (23.57309 59.82076) - (26.44457 60.64668)
CRS: +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs  
Number of fields: 10 
Number of list fields: 1 
                                                      name type length   typeName maxListCount
1                                                   gml_id    4      0     String            0
2                                            beginPosition    4     20     String            0
3                                              endPosition    4     20     String            0
4                                             timePosition    4     20     String            0
5                                                    value    4     31     String            0
6                                               identifier    0      0    Integer            0
7                                                     name    5      0 StringList            3
8                                                   region    4     11     String            0
9                                                     time    4     20     String            0
10 result|MeasurementTimeseries|point|MeasurementTVP|value    2      0       Real            0

Field name is of type stringList, which have been supported since OGR 1.8.0. Perhaps there's some low-level, GML-driver specific missing in sf which is causing st_read() to trip over?

[edzer]

The field is of type OFTStringList; currently, sf does not support any of the OFT*List types.

[jlehtoma]

Understood, will revert back to rgdal in that case.

[edzer]

I can't find good pointers to how to go through field lists in the OGR API, and find the relevant code in rgdal quite obscure. What does it mean that the name field (of type StringList) above has length zero? Are these things variable-length (width), and supposed to map into R list-columns? @rsbivand

[jlehtoma]

Can't comment on the zero-length, but FWIW, rgdal splits the StringList field into separate columns. Continuing from the code examples above:

# Write GML using rgdal
ogr_gml_file <- "ogr_out.gml"
rgdal::writeOGR(ogr_dat, ogr_gml_file, "layer1", driver = "GML")
rgdal::ogrInfo(ogr_gml_file)
Source: "ogr_out.gml", layer: "layer1"
Driver: GML; number of rows: 22 
Feature type: wkbPoint with 2 dimensions
Extent: (23.57309 59.82076) - (26.44457 60.64668)
Number of fields: 13 
                                                      name type length typeName
1                                                      fid    4      0   String
2                                                   gml_id    4      0   String
3                                            beginPosition    4      0   String
4                                              endPosition    4      0   String
5                                             timePosition    4      0   String
6                                                    value    4      0   String
7                                               identifier    0     16  Integer
8                                                    name1    4      0   String
9                                                    name2    4      0   String
10                                                   name3    4      0   String
11                                                  region    4      0   String
12                                                    time    4      0   String
13 result.MeasurementTimeseries.point.MeasurementTVP.value    2      0     Real

For example, the original StringList value (3:Raasepori Jussarö,-16000041,2757) for name is split into 3 cols: Raasepori Jussarö (name1), -16000041 (name2) and 2757 (name3).

I haven't really thought this through, but I think I'd prefer this behavior over a list-column.

[edzer]

This reads OFT*List types into a list-column, you can get these into regular fields e.g. by

x = st_read("fmi_test.gml")
cbind(x, do.call(rbind, x$name))

and you may want to rename the newly added columns.

I can follow your preference for having these in columns by default (rgdal behavior), but I don't see the point of these types when you know ahead the length and type of the list elements - why not make them single-field columns in the first place, then? When they have unequal length, having them distributed over columns with new names would seem terribly inconvenient. Also, checking all lengths before reading might be very expensive.

[jlehtoma]

Excellent, this works well for me. Thanks a lot for a very fast fix!

I can follow your preference for having these in columns by default (rgdal behavior), but I don't see the point of these types when you know ahead the length and type of the list elements - why not make them single-field columns in the first place, then? When they have unequal length, having them distributed over columns with new names would seem terribly inconvenient. Also, checking all lengths before reading might be very expensive.

I have no idea on what grounds the FMI has decided the data types for this particular data or more generally, how common the type is. To me, single-field columns would make much more sense too. Since splitting the field afterwards is fairly trivial, I would certainly prioritize performance here and leave the implementation as it is currently until more compelling reasons come along.