Inconsistencies in EPSG of returned layers

[tomjwebb]

Hi - I've noticed that the EPSG of layers is not always 4326 by default. For instance, running the example in the documentation here https://emodnet.github.io/emodnet.wfs/articles/emodnet.wfs.html#:~:text=chr%3E%2C%20layer_namespace%20%3Cchr%3E-,We,-are%20now%20ready

seabed_wfs_client <-
  emodnet_init_wfs_client(service = "seabed_habitats_general_datasets_and_products")

habitats_directive_layer_names <- c("art17_hab_1110", "art17_hab_1140", "art17_hab_1150")

habitats_directive_layers <- emodnet_get_layers(
  wfs = seabed_wfs_client,
  layers = habitats_directive_layer_names,
  reduce_layers = TRUE
)

Gives me the following:

> st_crs(habitats_directive_layers, parameters = TRUE)$srid
[1] "EPSG:3035"

Other seabed habitat layers are giving me a different projection, e.g.

eu_seamap23_finemud <- emodnet_get_layers(
  wfs = seabed_wfs_client,
  layers = "eusm2023_subs_full",
  cql_filter = "substrate = 'Fine mud'",
  reduce_layers = TRUE
)

Results in

> st_crs(eu_seamap23_finemud, parameters = TRUE)$srid
[1] "EPSG:3857"

I don't think this is a bug per se because as you say in the documentation, you can set the CRS when you load the layer - so this works:

eu_seamap23_finemud <- emodnet_get_layers(
  wfs = seabed_wfs_client,
  layers = "eusm2023_subs_full",
  cql_filter = "substrate = 'Fine mud'",
  crs = 4326,
  reduce_layers = TRUE
)
st_crs(eu_seamap23_finemud, parameters = TRUE)$srid
[1] "EPSG:4326"

So it may just be a case of recommending users explicitly set the CRS?

I'm using emodnet.wfs_2.0.2.9000 in R version 4.4.1 (2024-06-14), Platform: x86_64-apple-darwin20, Running under: macOS Sonoma 14.6.1.

[maelle]

Thank you @tomjwebb!

@salvafern I first thought this was a documentation issue but now I wonder:

• why do the layers from the same service have different CRS?
• should the default of the crs argument not be NULL. For instance it could be getOption("emodnet.wfs.crs", "EPSG:4326"). I'm unsure whether using an option is ok but at the same time it'd be a pain for users to input it for every call... unless we show an example with purrr::partial() in the docs.

Also, now that the package is under review, we should refrain from changes, so this could be done in a branch. Maybe the reviewers will have feedback on this aspect.

[salvafern]

This inconsistency is because the crs is defined for each layer, not for the whole server. In each layer you have a native crs (if available in the stored data) but also a declared crs, that we configure manually. See https://docs.geoserver.org/stable/en/user/data/webadmin/layers.html#coordinate-reference-systems

In practice, the native and declared crs are the same in most cases. I assume you get the declared crs by default.

I think this behaviour is the most correct, as you are retrieving the data as the providers have configure in geoserver. getOption("emodnet.wfs.crs", "EPSG:4326") sounds good but as something the users have to activate themselves. Also I would include an info message saying that the layer is being transformed. Let me know what you think!

Beware this applies only to geoserver - not sure how it works in OGC webservices made via python. But I believe the behaviour is similar.