Features returned from WFS are unexpectedly of geometry type CURVEPOLYGON about ows4r HOT 6 CLOSED

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

from ows4r.

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

from ows4r.

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

from ows4r.

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.

from ows4r.

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

from ows4r.

@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

from ows4r.