Features returned from WFS are unexpectedly of geometry type CURVEPOLYGON

[odgersn]

Hi,

I am trying to retrieve polygons from a WFS server using the following:

# Connect to WFS
wfs <- WFSClient$new("https://lris.scinfo.org.nz/services;key=[personal token]/wfs",
                     serviceVersion="2.0.0",
                     logger = "INFO")

# Query for polygons within a bounding box
polys <- wfs$getFeatures(typeName = "layer-48076",
                         cql_filter = "bbox(GEOMETRY,5138900,1568632,5172362,1610673)")

This works, but the geometry type of the returned features in my polys object is CURVEPOLYGON, which I did not expect. This is inconvenient because curvepolygons are not fully supported by sf (therefore I can't cast the features to a supported type).

When I test the same query in a web browser, the geometry type of the features in the result is gml:Polygon, so I assumed the geometry type of my polys object would be POLYGON or MULTIPOLYGON.

Are curvepolygons really the correct geometry type in this case, or did I make a mistake somewhere? Am I just misunderstanding something?

Thanks,

Nathan

[eblondel]

Hi @odgersn do you have the possibility to send me a test token by email so I could test your WFS and troubleshoot this? Thanks

[eblondel]

Hi @odgersn thanks for this, i had a look. Unfortunately, the management as CURVEPOLYGON comes from upstream from sf package / gdal driver for GML. The WFS 2.0.0 handles by default GML 3.2 as format, hence this geometry type, which I think is the correct behavior from the standard WFS/GML viewpoint, but then indeed we can't handle it properly with sf.

I tried to use a lower version of WFS with your service (1.0.0 or 1.1.1) but with no effect (1.0.0 seems unsupported, and 1.1.1 seems to return a 2.0.0 response).

However I found a simple workaround to make sure you can retrieve polygons, and this consists in asking the WFS to return result in GML v2, by adding outputFormat=GML2 to the getFeatures request:

polys <- wfs$getFeatures(
    typeName = "layer-48076",
    cql_filter = "bbox(GEOMETRY,5138900,1568632,5172362,1610673)", 
    outputFormat="GML2"
)

Let me know if it works for you

[eblondel]

@odgersn did you have a chance to test the above solution?

[odgersn]

Hi Emmanuel, I'm sorry for not getting back to you sooner. Yes your query works, but not quite as I expected it to work. I ran your query with the following amendment in order to set the crs correctly:

polys <- wfs$getFeatures(
  typeName = "layer-48076",
  cql_filter = "bbox(GEOMETRY,5138900,1568632,5172362,1610673)", 
  outputFormat="GML2"
) %>% 
  st_set_crs(2193)

It returns an sf data frame of POLYGON features. Excellent. But when I plot the geometry, the results are odd. This is what I see when I run plot(st_geometry(polys)):

But this is what I should see:

I took a closer look at head(polys) and noticed that the x's and y's in the sf bbox are around the wrong way. For example xmin should be 1552588 and ymin should be 5141942:

Simple feature collection with 6 features and 17 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 5141942 ymin: 1552588 xmax: 5160510 ymax: 1608106
projected CRS:  NZGD2000 / New Zealand Transverse Mercator 2000
             gml_id OBJECTID TYPE LEGEND LUC1C LUC1S LUC1UU LUCPL LUC2C LUC2S LUC2UU LCORRC LCORRS LCORRUU LCORRQ  LUC LCORR
1  layer-48076.8681     8681    n     00     6     e     10  <NA>  <NA>  <NA>   <NA>      6      e      10   <NA> 6e10  6e10
2  layer-48076.8824     8824    n     00     6     e     11  <NA>  <NA>  <NA>   <NA>      6      e      11   <NA> 6e11  6e11
3 layer-48076.12306    12306    n     00     6     s     10  <NA>  <NA>  <NA>   <NA>      6      s      10   <NA> 6s10  6s10
4  layer-48076.6668     6668    n     00     4     s      8  <NA>  <NA>  <NA>   <NA>      4      s       8   <NA> 4s 8  4s 8
5   layer-48076.843      843    l     00     l     a     ke  <NA>  <NA>  <NA>   <NA>      l      a      ke   <NA> lake  lake
6  layer-48076.1592     1592    n     00     2     w      1  <NA>  <NA>  <NA>   <NA>      2      w       1   <NA> 2w 1  2w 1
                        GEOMETRY
1 POLYGON ((5146474 1589339, ...
2 POLYGON ((5150556 1606905, ...
3 POLYGON ((5150857 1573000, ...
4 POLYGON ((5151114 1574463, ...
5 POLYGON ((5152045 1582189, ...
6 POLYGON ((5159722 1571210, ...

I took another look at the BBOX in the WFS query and I noticed that the coordinates were in the form ymin,xmin,ymax,xmax. I thought I must have got the x's and y's around the wrong way. So I tried again with the following query, where I rearranged the order of the BBOX parameters to be xmin,ymin,xmax,ymax:

polys <- wfs$getFeatures(
  typeName = "layer-48076",
  cql_filter = "bbox(GEOMETRY,1568632,5138900,1610673,5172362)", 
  outputFormat="GML2"
) %>% 
  st_set_crs(2193)

But this just returns an empty data frame, which I would expect to see if my original WFS BBOX had the coordinates in the correct order after all (because that would mean that the coordinates of the new BBOX were out of order, and don't intersect with the target area):

Simple feature collection with 0 features and 17 fields
bbox:           xmin: NA ymin: NA xmax: NA ymax: NA
projected CRS:  NZGD2000 / New Zealand Transverse Mercator 2000
 [1] gml_id   OBJECTID TYPE     LEGEND   LUC1C    LUC1S    LUC1UU   LUCPL    LUC2C    LUC2S    LUC2UU   LCORRC   LCORRS  
[14] LCORRUU  LCORRQ   LUC      LCORR    GEOMETRY

I don't know how else to diagnose what's going on. I've since deleted the test token; I'd have to reissue you another one if you want to test for yourself.

[eblondel]

@odgersn back from my holidays. yes please send me a test token by email I will test it

[eblondel]

@odgersn in #42 I've added the possibility to manage other WFS output formats, such as json or csv. This feature can improve a lot data retrieval as the default format (GML) is more verbose and time consuming to handle.

I've tried to use another format in your case, and it seems to fix the issue. To add an outputFormat, just add it to your getFeatures request, like this:

polys <- wfs$getFeatures(
    typeName = "layer-48076",
    cql_filter = "bbox(GEOMETRY,5138900,1568632,5172362,1610673)", outputFormat = "json"
)

BTW, we still have an issue with CRS which is not detected in this request. I will investigate this.

Let me know