Confusion on how to use spatial ECQL predicates about emodnetwfs

[maelle]

Closing this as reading the issue thread, it seems to me that #132 was enough. @salvafern please re-open if not.

[annakrystalli]

Test cql filtering on physics which is not using geoserver

[annakrystalli]

So the hunch was correct!!

library(magrittr)

wfs_gs <- EMODnetWFS::emodnet_init_wfs_client(service = "geology_seabed_substrate_maps")
#> Loading ISO 19139 XML schemas...
#> Loading ISO 19115 codelists...
#> Loading IANA mime types...
#> No encoding supplied: defaulting to UTF-8.
#> ✓ WFS client created succesfully
#> ℹ Service: 'https://drive.emodnet-geology.eu/geoserver/gtk/wfs'
#> ℹ Version: '2.0.0'

ref_sf <- EMODnetWFS::emodnet_get_layers(wfs = wfs_gs, layers = "seabed_substrate_1m",
                       cql_filter ="id=10848", reduce_layers = TRUE) 


# create reference map
ref_map <- ref_sf %>% 
    sf::st_cast(to = "GEOMETRYCOLLECTION") %>% # mapview doesn't seem to work with 
    # EWKT so MULTIURFACE needs casting to a geometry collection
mapview::mapview(zcol = "folk_7cl_txt", burst = TRUE)
#> Warning: multiple methods tables found for 'direction'
#> Warning: multiple methods tables found for 'gridDistance'
#> Warning: multiple methods tables found for 'crop'
#> Warning: multiple methods tables found for 'extend'

Transform to GEOMETRYCOLLECTION.

ref_geom_colllection <- ref_sf %>% 
    sf::st_cast(to = "GEOMETRYCOLLECTION") %>% 
  sf::st_geometry() %>% sf::st_as_text()

Use WKT to construct DWITHIN query -

geom_colllection_res <- EMODnetWFS::emodnet_get_layers(wfs = wfs_gs, layers = "seabed_substrate_1m",
                       cql_filter = paste0("DWITHIN(geom,", ref_geom_colllection, ", 0.5, kilometers)"),
                       reduce_layers = TRUE ) 


#  returns very strange results!
ref_map + geom_colllection_res %>% 
    sf::st_cast(to = "GEOMETRYCOLLECTION") %>%
mapview::mapview(zcol = "folk_7cl_txt", burst = TRUE)

ref_geom_colllection
#> [1] "GEOMETRYCOLLECTION (POLYGON ((3184146 2544100, 3184141 2544099, 3184108 2544110, 3184055 2544148, 3184014 2544192, 3183974 2544250, 3183940 2544317, 3183914 2544387, 3183899 2544453, 3183890 2544510, 3183864 2544808, 3183730 2546928, 3183967 2546943, 3184058 2545528, 3184057 2545520, 3184054 2545456, 3184060 2545400, 3184057 2545336, 3184054 2545271, 3184058 2545205, 3184058 2545118, 3184090 2545026, 3184103 2544826, 3184102 2544806, 3184107 2544751, 3184098 2544557, 3184103 2544360, 3184103 2544317, 3184139 2544137, 3184146 2544100)))"

Take the first point of geom_colllection_res. Keep coordinates as returned by original request.

wkt <- wellknown::geojson2wkt(list(Point = c(3184146, 2544100)))

# SAME RESULT
EMODnetWFS::emodnet_get_layers(wfs = wfs_gs, layers = "seabed_substrate_1m",
                       cql_filter = paste0("DWITHIN(geom,", wkt, ", 0.5, kilometers)"),
                       reduce_layers = TRUE )  %>% 
    sf::st_cast(to = "GEOMETRYCOLLECTION") %>%
mapview::mapview(zcol = "folk_7cl_txt", burst = TRUE)

swap the coordinates

wkt <- wellknown::geojson2wkt(list(Point = c(2544100, 3184146)))

# WORKS CORRECTLY!! Although the distance now feels off
EMODnetWFS::emodnet_get_layers(wfs = wfs_gs, layers = "seabed_substrate_1m",
                       cql_filter = paste0("DWITHIN(geom,", wkt, ", 0.5, kilometers)"),
                       reduce_layers = TRUE )  %>% 
    sf::st_cast(to = "GEOMETRYCOLLECTION") %>%
mapview::mapview(zcol = "folk_7cl_txt", burst = TRUE)

^{Created on 2022-03-01 by the reprex package (v2.0.1)}

[annakrystalli]

The result of BEYOND filter

[annakrystalli]

Question in email

We’ve spend the last couple of days working on the EMODnetWFS package, myself in particular on figuring out why the spatial cql filtering wasn't correctly. There’s complete details in this issue: #30 but the conclusion was that at least in the first service I tested, using a WKT geom (extracted from an sf object created from a response form the server itself) in spatial filters, the coordinates seem to be interpreted as y x in the request and flipping them returns the correct result.

I don’t know if this is common practice, whether it is a setting of the databases or geoserver and how widespread it is across servers. I know you mentioned it’s a trial and error path in your previous response and I’ve written a script to try every service and figure out whether x,y or y,x coordinates return the correct results but I’m having to handle a lot of exceptions and this seems like an unnecessarily laborious way to figure this out.

Do you have any advice on how I might be able to get access to this information in a more authoritative manner? Also, if we're having such issues, it is likely that other developers trying to build on the web services might have similar problems. I was therefore wondering wether a more formal logging of the tech stack and configuration for each would make the services much easier for developers to build on therefore increasing their impact?

Any advice would be greatly appreciated.

[annakrystalli]

@LennertSchepers 's

If it can be of any help: depending on the WFS version, the axis ordering can be different. So WFS 1.0.0 vs WFS 1.1.0. I think this explains it quite well: https://docs.geoserver.org/latest/en/user/services/wfs/axis_order.html

I’ve done some detailed google searches in the past so it’s really not so easy. I remember that years ago the community was looking for a successor of WFS, the ‘OGC API – Features’, to have a more standardised version compared to WFS. I think specifications are out (https://ogcapi.ogc.org/features/) but I haven’t seen it in practice - but note that I haven’t opened a GIS system (other than R :p) in 2 years so my knowledge might be a bit outdated.

[annakrystalli]

What is the best way to find out which versions each server supports?

/wfs?request=getcapabilities

[LennertSchepers]

What is the best way to find out which versions each server supports?

/wfs?request=getcapabilities

yes, correct

[annakrystalli]

Thanks @LennertSchepers ! Sorry if I caused confusion, I was just pasting responses from an email thread for our records here.

@maelle I was thinking that our users don't really need to know about different WFS versions. What we could do as the simplest solution is fix the version the package uses to communicate with the servers (I believe the current default version is available on all servers), remove the version argument and build our backend around a predictable response.

@LennertSchepers one question I had is whether service administrators are able to override the standard version behaviour and whether these configurations are accessible through the service?

[LennertSchepers]

@annakrystalli could you give an example of what you mean?

[annakrystalli]

Hey @LennertSchepers

So, what we want to do is for an R user to be able to supply an sf object as the geom to any spatial filter. Currently, at least with WFS version 2.0.1, the coordinates in the geom are flipped compared to the server, so any spatial filters which use a WKT representation of an R sf object return incorrect results of what the user expects (hence this issue). However according to your note in the email thread, this behaviour could differ if different WFS versions are used.

My initial idea was to build a small package that could translate an sf object to a correct WKT representation for any given version and be able to correctly plug it in to a spatial cql_filter. On the other hand, a simpler solution described in my comment above, would be to fix the version the package uses to communicate with the server, making the expectation of the server more predictable and simple to deal with within the package.

Does this make sense?

[annakrystalli]

@LennertSchepers so my question above was whether server administrators could override the default axis ordering of a version through configs and if so, whether those are available as metadata through request to the server.

[LennertSchepers]

I'm not sure about this - but I would assume not, as this seems to be specified in the WFS standard. From this page it seems that the axis ordering can only be doubtful for WFS 1.0.0.

Just wondering about it, do you know if the problem is situated at

• the axis ordering of the layer hosted on the WFS, or
• the axis ordering of the geometry that you provide to filter on (the provided wkt string)? It's mentioned here that (E)CQL provides no mechanism for specifying the projection of a geometry and always assumes that it is in the same projection as the data being queried. The default projection is mentioned in

[maelle]

@LennertSchepers is your comment incomplete? The last sentence has no end. 😸 (thanks for all the information in any case!)

[maelle]

I wonder how to get the supported versions for the marine litter service. The others all seem to support the same versions?

get_version <- function(service_url) {
    if (service_url == "https://www.ifremer.fr/services/wfs/emodnet_chemistry2") {
        return(NA_character_)
    }
    httr2::request(service_url) |>
        httr2::req_url_query(request = "getCapabilities") |>
        httr2::req_perform() |>
        httr2::resp_body_xml() |>
        xml2::xml_find_first("//ows:OperationsMetadata") |>
        xml2::xml_find_first("//ows:Operation[@name='GetCapabilities']") |>
        xml2::xml_find_first("//ows:Parameter[@name='AcceptVersions']") |>
        xml2::xml_find_first("ows:AllowedValues") |>
        xml2::xml_find_all("ows:Value") |>
        purrr::map_chr(xml2::xml_text)
}

services <- EMODnetWFS::emodnet_wfs()

tibble::tibble(
    service = services$service_name,
    versions = purrr::map(services$service_url, get_version)
) |>
    knitr::kable()

service	versions
bathymetry	1.0.0, 1.1.0, 2.0.0
biology	1.0.0, 1.1.0, 2.0.0
biology_occurrence_data	1.0.0, 1.1.0, 2.0.0
chemistry_cdi_data_discovery_and_access_service	1.0.0, 1.1.0, 2.0.0
chemistry_cdi_distribution_observations_per_category_and_region	1.0.0, 1.1.0, 2.0.0
chemistry_contaminants	1.0.0, 1.1.0, 2.0.0
chemistry_marine_litter	NA
geology_coastal_behavior	1.0.0, 1.1.0, 2.0.0
geology_events_and_probabilities	1.0.0, 1.1.0, 2.0.0
geology_marine_minerals	1.0.0, 1.1.0, 2.0.0
geology_sea_floor_bedrock	1.0.0, 1.1.0, 2.0.0
geology_seabed_substrate_maps	1.0.0, 1.1.0, 2.0.0
geology_submerged_landscapes	1.0.0, 1.1.0, 2.0.0
human_activities	1.0.0, 1.1.0, 2.0.0
physics	1.0.0, 1.1.0, 2.0.0
seabed_habitats_general_datasets_and_products	1.0.0, 1.1.0, 2.0.0
seabed_habitats_individual_habitat_map_and_model_datasets	1.0.0, 1.1.0, 2.0.0

^{Created on 2022-10-14 with reprex v2.0.2}

[maelle]

Related #126

[maelle]

I can't find my notes but I remember @annakrystalli and I discussed it'd be easier to call all services with one single version as it'd mean the coordinates would be in a known order. 🤔

If that is correct then

• how to get the missing versions for a service #30 (comment)
• what version to choose

[annakrystalli]

I think the best option would be to fix the latest version 2.0.0 (which is also currently the default) which I believe all servers can use.

[maelle]

So instead of surfacing service_version we'd set it to 2.0.0? That's what you mean by fix, making constant, as opposed to repairing, correct? Thank you!

[salvafern]

I'm a bit hesitant but if it makes the development easier go ahead :) we can always reconsider later upon request.

I guess the only issue could be if someone wants to integrate EMODnetWFS with some older system that use a WFS version lower than 2.0.0.

But this type of user will probably go directly to ows4R rather than using the EMODnetWFS client.

[maelle]

However fixing the version is probably only part of the fix? I'm still confused (:sweat_smile:) as to whether we need to flip coordinates somewhere for something.

[maelle]

i.e. what Anna meant by "build our backend around a predictable response.".

[salvafern]

In principle if we have been always using 2.0.0 and the axis order defined by the EPGS definition it shouldn't raise much issues? e.g see the order of coordinates for each version: https://docs.geoserver.org/2.22.x/en/user/services/wfs/axis_order.html

WFS 1.0.0: Provides geographic coordinates in east/north and may not be trusted to respect the EPSG definition axis order.
WFS 1.1.0: Respects the axis order defined by the EPSG definition.
WFS 2.0.0: Respects the axis order defined by the EPSG definition.

[maelle]

so #132 would be enough?