WFS query returns unhandy MULTISURFACE geometries

[fdetsch]

I was wondering why some WFS queries in {ows4R} return MULTISURFACE geometries. As documented numerously (e.g. r-spatial/sf#748, How to simplify a MULTISURFACE [...] geometry? - Stack Overflow), these can be tedious to handle in R as even very basic operations like

• sf::st_coordinates(),
• as(..., "Spatial"),
• terra::vect() and
• plot()

are not supported. For demonstration purposes, take the following example:

cli = ows4R::WFSClient$new(
  "https://geodatendienste-landesplanung-hessen.de/geoserver/suedhessen/wfs"
  , serviceVersion = "2.0.0"
)

features = cli$getFeatures(
  "suedhessen:gruenzug"
)

sf::st_geometry_type(features[1L, ])
# [1] MULTISURFACE

plot(features[1L, ])

Error in xy.coords(x, y, xlabel, ylabel, log) :
  'x' is a list, but does not have components 'x' and 'y'

For reference, here is the corresponding {owslib} workflow in Python that produces regular MultiPolygon features:

from geopandas import read_file
from owslib.wfs import WebFeatureService

cli = WebFeatureService(
    'https://geodatendienste-landesplanung-hessen.de/geoserver/suedhessen/wfs',
    version='2.0.0'
)

response = cli.getfeature('suedhessen:gruenzug')
result = read_file(response)

result.iloc[[0]].geom_type
# 0    MultiPolygon
# dtype: object

result.iloc[[0]].plot()

(Note: my intention is not to pit Python and R against each other, and {owslib} just serves for demonstration purposes)

[eblondel]

Answer is simple: WFS 2.0 delivers multisurface and not multipolygons, so that is default behaviour that sticks with the OGC standard. If the python library delivers multiPolygons this is a deviation from the standard they might have developed an additional logic to revert to multipolygons in the response. Le jeu. 19 déc. 2024, 11:24, Florian Detsch ***@***.***> a écrit :

…

I was wondering why some WFS queries in {ows4R} return MULTISURFACE geometries. As documented numerously (e.g. r-spatial/sf#748 <r-spatial/sf#748>, How to simplify a MULTISURFACE [...] geometry? - Stack Overflow <https://stackoverflow.com/questions/68810009/how-to-simplify-a-multisurface-curvepolygon-or-compoundcurve-geometry>), these can be tedious to handle in R as even very basic operations like - sf::st_coordinates(), - as(..., "Spatial"), - terra::vect() and - plot() are not supported. For demonstration purposes, take the following example: cli = ows4R::WFSClient$new( "https://geodatendienste-landesplanung-hessen.de/geoserver/suedhessen/wfs" , serviceVersion = "2.0.0" ) features = cli$getFeatures( "suedhessen:gruenzug" ) sf::st_geometry_type(features[1L, ])# [1] MULTISURFACE plot(features[1L, ]) Error in xy.coords(x, y, xlabel, ylabel, log) : 'x' is a list, but does not have components 'x' and 'y' For reference, here is the corresponding {owslib} workflow in Python that produces regular MultiPolygon features: from geopandas import read_filefrom owslib.wfs import WebFeatureService cli = WebFeatureService( 'https://geodatendienste-landesplanung-hessen.de/geoserver/suedhessen/wfs', version='2.0.0' ) response = cli.getfeature('suedhessen:gruenzug')result = read_file(response) result.iloc[[0]].geom_type# 0 MultiPolygon# dtype: object result.iloc[[0]].plot() image.png (view on web) <https://github.com/user-attachments/assets/c9c64641-309c-4900-8604-42c882991668> (Note: my intention is not to pit Python and R against each other, and {owslib} just serves for demonstration purposes) — Reply to this email directly, view it on GitHub <#132>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAKDK3AII25EISCQOYIJYTL2GKNFPAVCNFSM6AAAAABT4U5OJWVHI2DSMVQWIX3LMV43ERDJONRXK43TNFXW4OZXG4YTGOJRGI> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[fdetsch]

Is this specific to the WFS from above? For example, the code in #131 produces MULTIPOLYGON geometries in R although the service version is also 2.0.

[eblondel]

It is part of the WFS specification, version 2.0, that returns GML 3 as default format. For complex geometries, you will get geometries in this format. For the 'multipolygons', that will be MultiSurface. In the example you mentioned you do not have MULTIPOLYGON but POLYGON which is a simple geometry.

[fdetsch]

Forgot to mention, you need to set maxfeatures = 5000L in order for the #131 code to return MULTIPOLYGON geometries:

features = cli$getFeatures(
  "adv:AX_Bodenschaetzung"
  , maxfeatures = 5000L
)

unique(
  sf::st_geometry_type(features)
)
# [1] MULTIPOLYGON

According to your logic, shouldn't this be MULTISURFACE?

[eblondel]

To my logic, no, to the WFS 2.0 standard. I don't know what is the WFS implementation behind this server, it may be that set it more permissive and returning GML2 instead of GML3 by default.
I looked at the server, and it seems it accepts other output formats, and this is supported by ows4R . You can try application/json it will work, and you will not depend on GML3 /MULTISURFACE format.

Just add the outputFormat as argument to your getFeatures call.

features = cli$getFeatures(
  "suedhessen:gruenzug",
 outputFormat = "application/json"
)

[fdetsch]

Marvellous, thanks! Start tag expected, '<' not found message aside that appears when running your code, this produces exactly what I was looking for.

[eblondel]

Yes, I get the same message. It's not blocking but I will look at it.

[eblondel]

FYI, this small issue has been fixed with #135