st_read(..., drivers = "GMLAS") fails to read CityGML data

[dimfalk]

I was messing around with CityGML LoD1 data lately and was looking for a way to deal with this data using R. A blog post suggests using GDAL's GMLAS driver to read this data. Since GMLAS does not seem to be included in the GDAL distribution being shipped with sf on Windows - hope I got the dependencies right, verbally - I switched to Ubuntu 22.04 to check this out.

(By the way, any chance GMLAS will be available for sf on Windows sometime?)

Making use of GDAL directly seems to work since the file structure is being recognized:

ogrinfo --version
> GDAL 3.8.4, released 2024/02/08

# per default, GML is being used and detects only one layer
ogrinfo LoD1_32_280_5657_1_NW.gml 
> INFO: Open of `LoD1_32_280_5657_1_NW.gml'
>       using driver `GML' successful.
> Metadata:
>   NAME=LoD1_32_280_5657_1_NW
> 1: Building (3D PolyhedralSurface)

# enforcing GMLAS gets you 150 layers
ogrinfo GMLAS:LoD1_32_280_5657_1_NW.gml 
> INFO: Open of `GMLAS:LoD1_32_280_5657_1_NW.gml'
>      using driver `GMLAS' successful.
> 1: buildingpart (PolyhedralSurface, Multi Curve)
> 2: buildingpart_metadataproperty (None)
> 3: buildingpart_name (None)
> 4: building (PolyhedralSurface, Multi Curve)
> 5: building_metadataproperty (None)
> [...]

This part can be reproduced executing the commands from {sf}:

library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.0; sf_use_s2() is TRUE

# get some data
f <- "LoD1_32_280_5657_1_NW.gml"
utils::download.file(paste0("https://www.opengeodata.nrw.de/produkte/geobasis/3dg/lod1_gml/lod1_gml/", f),
                     f)

sf::gdal_utils(util = "ogrinfo", source = f)
#> INFO: Open of `LoD1_32_280_5657_1_NW.gml'
#>       using driver `GML' successful.
#> Metadata:
#>   NAME=LoD1_32_280_5657_1_NW
#> 1: Building (3D PolyhedralSurface)

sf::gdal_utils(util = "ogrinfo", source = paste0("GMLAS:", f))
#> INFO: Open of `GMLAS:LoD1_32_280_5657_1_NW.gml'
#>      using driver `GMLAS' successful.
#> 1: buildingpart (PolyhedralSurface, Multi Curve)
#> 2: buildingpart_metadataproperty (None)
#> 3: buildingpart_name (None)
#> 4: building (PolyhedralSurface, Multi Curve)
#> 5: building_metadataproperty (None)

Describe the bug
Trying to read this file using st_read(), the GMLAS driver does not seem to be supported.

To Reproduce

# GML is used per default, like above; seems to work, but only recognizes one layer
sf::st_read(f)
#> Reading layer `Building' from data source `/root/LoD1_32_280_5657_1_NW.gml' using driver `GML'
#> Simple feature collection with 81 features and 22 fields
#> Geometry type: POLYHEDRALSURFACE
#> Dimension:     XYZ
#> Bounding box:  xmin: 280851.8 ymin: 5657120 xmax: 281004.9 ymax: 5657728
#> z_range:       zmin: 35.74 zmax: 47.632
#> Projected CRS: ETRS89 / UTM zone 32N

# enforcing GMLAS, reading fails
sf::st_read(f, drivers = "GMLAS")
#> options:        GMLAS 
#> Error: Cannot open "/root/LoD1_32_280_5657_1_NW.gml"; The source could be corrupt or not supported. See `st_drivers()` for a list of supported formats.

# the driver should be supported, as far as I understand
"GMLAS" %in% sf::st_drivers()$name
#> TRUE

R version 4.1.2 (2021-11-01) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Ubuntu 22.04.4 LTS

[edzer]

@rouault looking at the problem above, I get a segfault in GetFeatureCount with this minimal C program:

#include <stdio.h>
#include <ogrsf_frmts.h>

int main() {
	OGRRegisterAll();
	auto poDS = (GDALDataset *) GDALOpenEx("GMLAS:LoD1_32_280_5657_1_NW.gml", GDAL_OF_VECTOR | GDAL_OF_READONLY,
        NULL, NULL, NULL );
	if (poDS == NULL) {
		printf("cannot open file\n");
		return 0;
	}
	auto poLayer =   poDS->GetLayerByName("address1");
	if (poLayer == NULL) {
		printf("cannot open layer\n");
		return 0;
	}
	size_t n = poLayer->GetFeatureCount();
	printf("size: %ld\n", n);
	return 0;
}

The file is found here. Is this a logic problem on my side, or a problem on GDAL's side?

[rouault]

or a problem on GDAL's side?

Yes, a bug in the OGR GMLAS driver. Fix in OSGeo/gdal#9590.
A workaround is to call poLayer->GetLayerDefn() before GetFeatureCount()

[edzer]

Great, thanks!

[dimfalk]

Thank you very much!

@edzer Allow me this concluding beginner question since I'm not aware of the general workflow here:

When can we roughly expect the release of a new GDAL version that includes this bug fix (just out of curiousity)? And is it generally expected that the GMLAS driver will sometime be included in GDAL, which is used by sf on Windows?

[rouault]

@dimfalk you're lucky. This was the last fix to be in the queue for the imminent GDAL 3.8.5 release

[dimfalk]

@rouault Hooray, lucky me! Thanks!

[edzer]

@dimfalk I implemented the workaround hinted to by @rouault in 43aec36 so things should now work when you install sf directly from github (or r-universe).

[dimfalk]

@edzer Yeah, I noticed that - will give it a try, thanks again! I was just curious how fast dependencies are updated in general, just to get an idea about lifecycle management. Not crucial regarding this issue. But what about GMLAS for GDAL on Windows? 😏

[dimfalk]

Just tried to repeat the steps described above using sf 1.0.17, but still failing, unfortunately. Am I missing something?

devtools::install_github("r-spatial/sf")

packageVersion("sf")
#> [1] '1.0.17'

library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE

sf::st_read(f, drivers = "GMLAS")
#> options:        GMLAS 
#> Error: Cannot open "/root/LoD1_32_280_5657_1_NW.gml"; The source could be corrupt or not supported. See `st_drivers()` for a list of supported formats.

[edzer]

Looks like you have to prepend GMLAS: to the file name.

[dimfalk]

Oh, I assumed this is done underneath making use of drivers = "GMLAS". Thanks!

[StefanCiesilski]

@dimfalk Thanks for testing with CityGML LoD1. I saw the same blog post but didn't get any further. Have you already tried it with CityGML LoD2 data, too?

[dimfalk]

@StefanCiesilski: Not yet, sorry!

[StefanCiesilski]

@dimfalk OK!