`OWSBoundingBox` - `getBBOX` returning incorrect crs
annakrystalli opened this issue · 5 comments
annakrystalli commented
Hi @eblondel ,
I noticed that the crs information appended to bbox objects returned by getBBOX are incorrect. The crs seems to default to EPSG:4326 so when the coverage crs deviates from that it returns the correct coordinates but attaches the wrong crs to the crs bbox attribute. In the example below, when using getBBOX on the BoundingBox returned by getBoundingBox(), the coordinates are different but the crs attached is the same as when using getWGS84BoundingBox() (which correctly returns EPSG:4326). The correct crs should be EPSG:3857, as indicated when interrogating the BoundingBox attrs themselves.
library(ows4R)
#> Loading required package: geometa
#> Loading ISO 19139 XML schemas...
#> Loading ISO 19115 codelists...
#> Loading IANA mime types...
#> No encoding supplied: defaulting to UTF-8.
#> Loading required package: keyring
seabed <- WCSClient$new(url = "https://ows.emodnet-seabedhabitats.eu/geoserver/emodnet_open_maplibrary/wcs", serviceVersion = "2.0.1", logger = "DEBUG")
#> [ows4R][INFO] OWSGetCapabilities - Fetching https://ows.emodnet-seabedhabitats.eu/geoserver/emodnet_open_maplibrary/wcs?service=WCS&version=2.0.1&request=GetCapabilities
cov <- seabed$getCapabilities()$findCoverageSummaryById("emodnet_open_maplibrary__mediseh_maerl", exact = TRUE)
bbox <- cov$getBoundingBox()$BoundingBox$getBBOX()
bbox
#> xmin ymin xmax ymax
#> -729957.1 3431706.8 4782842.9 5991306.8
attr(bbox, "crs")
#> Coordinate Reference System:
#> User input: EPSG:4326
#> wkt:
#> GEOGCRS["WGS 84",
#> ENSEMBLE["World Geodetic System 1984 ensemble",
#> MEMBER["World Geodetic System 1984 (Transit)"],
#> MEMBER["World Geodetic System 1984 (G730)"],
#> MEMBER["World Geodetic System 1984 (G873)"],
#> MEMBER["World Geodetic System 1984 (G1150)"],
#> MEMBER["World Geodetic System 1984 (G1674)"],
#> MEMBER["World Geodetic System 1984 (G1762)"],
#> ELLIPSOID["WGS 84",6378137,298.257223563,
#> LENGTHUNIT["metre",1]],
#> ENSEMBLEACCURACY[2.0]],
#> PRIMEM["Greenwich",0,
#> ANGLEUNIT["degree",0.0174532925199433]],
#> CS[ellipsoidal,2],
#> AXIS["geodetic latitude (Lat)",north,
#> ORDER[1],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> AXIS["geodetic longitude (Lon)",east,
#> ORDER[2],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> USAGE[
#> SCOPE["Horizontal component of 3D system."],
#> AREA["World."],
#> BBOX[-90,-180,90,180]],
#> ID["EPSG",4326]]
cov$getBoundingBox()$BoundingBox$attrs
#> $crs
#> [1] "http://www.opengis.net/def/crs/EPSG/0/EPSG:3857"
bboxWGS84 <- cov$getWGS84BoundingBox()$WGS84BoundingBox$getBBOX()
bboxWGS84
#> xmin ymin xmax ymax
#> -6.557316 29.439519 42.965009 47.300773
attr(bboxWGS84, "crs")
#> Coordinate Reference System:
#> User input: EPSG:4326
#> wkt:
#> GEOGCRS["WGS 84",
#> ENSEMBLE["World Geodetic System 1984 ensemble",
#> MEMBER["World Geodetic System 1984 (Transit)"],
#> MEMBER["World Geodetic System 1984 (G730)"],
#> MEMBER["World Geodetic System 1984 (G873)"],
#> MEMBER["World Geodetic System 1984 (G1150)"],
#> MEMBER["World Geodetic System 1984 (G1674)"],
#> MEMBER["World Geodetic System 1984 (G1762)"],
#> ELLIPSOID["WGS 84",6378137,298.257223563,
#> LENGTHUNIT["metre",1]],
#> ENSEMBLEACCURACY[2.0]],
#> PRIMEM["Greenwich",0,
#> ANGLEUNIT["degree",0.0174532925199433]],
#> CS[ellipsoidal,2],
#> AXIS["geodetic latitude (Lat)",north,
#> ORDER[1],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> AXIS["geodetic longitude (Lon)",east,
#> ORDER[2],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> USAGE[
#> SCOPE["Horizontal component of 3D system."],
#> AREA["World."],
#> BBOX[-90,-180,90,180]],
#> ID["EPSG",4326]]Created on 2022-08-09 by the reprex package (v2.0.1)
eblondel commented
Thanks @annakrystalli for the comprehensive report. I think that's something for which i forgot to add a 'TODO', and EPSG:4326 was still hardcoded in the getBBOX() of OWSBoundingBox (common to all OGC WxS capabilities). I've pushed some code for you to test. let me know
annakrystalli commented
Great, I'm back on the project tomorrow so will take a look and let you know!
annakrystalli commented
Just checked and all is working as expected!
annakrystalli commented
I've done some more work and noticed that there might now be a problem with the coordinates returned now through getBBOX(). They now do not match the boundedBy limits even though the crs seems correct:
library(ows4R)
#> Loading required package: geometa
#> Loading ISO 19139 XML schemas...
#> Loading ISO 19115 codelists...
#> Loading IANA mime types...
#> No encoding supplied: defaulting to UTF-8.
#> Loading required package: keyring
wcs <- ows4R::WCSClient$new(url = "https://ows.emodnet-humanactivities.eu/wcs",
serviceVersion = "2.0.1", logger = "DEBUG")
#> [ows4R][INFO] OWSGetCapabilities - Fetching https://ows.emodnet-humanactivities.eu/wcs?service=WCS&version=2.0.1&request=GetCapabilities
summary <- wcs$getCapabilities()$findCoverageSummaryById("emodnet__2017_01_st_00")
bbox <- summary$BoundingBox$BoundingBox$getBBOX()
bbox
#> xmin ymin xmax ymax
#> -9763441 1684719 10868459 22822360
attr(bbox, "crs")
#> Coordinate Reference System:
#> User input: EPSG:3857
#> wkt:
#> PROJCRS["WGS 84 / Pseudo-Mercator",
#> BASEGEOGCRS["WGS 84",
#> ENSEMBLE["World Geodetic System 1984 ensemble",
#> MEMBER["World Geodetic System 1984 (Transit)"],
#> MEMBER["World Geodetic System 1984 (G730)"],
#> MEMBER["World Geodetic System 1984 (G873)"],
#> MEMBER["World Geodetic System 1984 (G1150)"],
#> MEMBER["World Geodetic System 1984 (G1674)"],
#> MEMBER["World Geodetic System 1984 (G1762)"],
#> ELLIPSOID["WGS 84",6378137,298.257223563,
#> LENGTHUNIT["metre",1]],
#> ENSEMBLEACCURACY[2.0]],
#> PRIMEM["Greenwich",0,
#> ANGLEUNIT["degree",0.0174532925199433]],
#> ID["EPSG",4326]],
#> CONVERSION["Popular Visualisation Pseudo-Mercator",
#> METHOD["Popular Visualisation Pseudo Mercator",
#> ID["EPSG",1024]],
#> PARAMETER["Latitude of natural origin",0,
#> ANGLEUNIT["degree",0.0174532925199433],
#> ID["EPSG",8801]],
#> PARAMETER["Longitude of natural origin",0,
#> ANGLEUNIT["degree",0.0174532925199433],
#> ID["EPSG",8802]],
#> PARAMETER["False easting",0,
#> LENGTHUNIT["metre",1],
#> ID["EPSG",8806]],
#> PARAMETER["False northing",0,
#> LENGTHUNIT["metre",1],
#> ID["EPSG",8807]]],
#> CS[Cartesian,2],
#> AXIS["easting (X)",east,
#> ORDER[1],
#> LENGTHUNIT["metre",1]],
#> AXIS["northing (Y)",north,
#> ORDER[2],
#> LENGTHUNIT["metre",1]],
#> USAGE[
#> SCOPE["Web mapping and visualisation."],
#> AREA["World between 85.06°S and 85.06°N."],
#> BBOX[-85.06,-180,85.06,180]],
#> ID["EPSG",3857]]
summary$getDescription()$boundedBy
#> [ows4R][INFO] WCSDescribeCoverage - Fetching https://ows.emodnet-humanactivities.eu/wcs?service=WCS&version=2.0.1&coverageId=emodnet__2017_01_st_00&request=DescribeCoverage
#> Loading required package: sf
#> Linking to GEOS 3.9.1, GDAL 3.4.0, PROJ 8.1.1; sf_use_s2() is TRUE
#> <GMLEnvelope>
#> ....|-- lowerCorner: -622999.999999999 604000.000000006
#> ....|-- upperCorner: 6885000 7035000.00000001
summary$getBoundingBox()$BoundingBox$attrs$crs
#> [1] "http://www.opengis.net/def/crs/EPSG/0/EPSG:3857"Created on 2022-08-17 by the reprex package (v2.0.1)
Session info
sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#> setting value
#> version R version 4.1.3 (2022-03-10)
#> os macOS Big Sur/Monterey 10.16
#> system x86_64, darwin17.0
#> ui X11
#> language (EN)
#> collate en_GB.UTF-8
#> ctype en_GB.UTF-8
#> tz Europe/Athens
#> date 2022-08-17
#> pandoc 2.18 @ /Applications/RStudio.app/Contents/MacOS/quarto/bin/tools/ (via rmarkdown)
#>
#> ─ Packages ───────────────────────────────────────────────────────────────────
#> package * version date (UTC) lib source
#> askpass 1.1 2019-01-13 [1] CRAN (R 4.1.0)
#> assertthat 0.2.1 2019-03-21 [1] CRAN (R 4.1.0)
#> class 7.3-20 2022-01-16 [1] CRAN (R 4.1.3)
#> classInt 0.4-7 2022-06-10 [1] CRAN (R 4.1.2)
#> cli 3.3.0 2022-04-25 [1] CRAN (R 4.1.2)
#> codetools 0.2-18 2020-11-04 [1] CRAN (R 4.1.3)
#> curl 4.3.2 2021-06-23 [1] CRAN (R 4.1.0)
#> DBI 1.1.3 2022-06-18 [1] CRAN (R 4.1.2)
#> digest 0.6.29 2021-12-01 [1] CRAN (R 4.1.0)
#> dplyr 1.0.9 2022-04-28 [1] CRAN (R 4.1.2)
#> e1071 1.7-11 2022-06-07 [1] CRAN (R 4.1.2)
#> evaluate 0.15 2022-02-18 [1] CRAN (R 4.1.2)
#> fansi 1.0.3 2022-03-24 [1] CRAN (R 4.1.2)
#> fastmap 1.1.0 2021-01-25 [1] CRAN (R 4.1.0)
#> fs 1.5.2 2021-12-08 [1] CRAN (R 4.1.0)
#> generics 0.1.2 2022-01-31 [1] CRAN (R 4.1.2)
#> geometa * 0.6-6 2022-01-26 [1] CRAN (R 4.1.2)
#> glue 1.6.2 2022-02-24 [1] CRAN (R 4.1.2)
#> highr 0.9 2021-04-16 [1] CRAN (R 4.1.0)
#> htmltools 0.5.3 2022-07-18 [1] CRAN (R 4.1.2)
#> httr 1.4.3 2022-05-04 [1] CRAN (R 4.1.2)
#> jsonlite 1.8.0 2022-02-22 [1] CRAN (R 4.1.2)
#> KernSmooth 2.23-20 2021-05-03 [1] CRAN (R 4.1.3)
#> keyring * 1.3.0 2021-11-29 [1] CRAN (R 4.1.0)
#> knitr 1.39 2022-04-26 [1] CRAN (R 4.1.2)
#> lifecycle 1.0.1 2021-09-24 [1] CRAN (R 4.1.0)
#> magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.1.2)
#> openssl 2.0.2 2022-05-24 [1] CRAN (R 4.1.2)
#> ows4R * 0.3 2022-08-17 [1] Github (eblondel/ows4R@1cd9dbc)
#> pillar 1.8.0 2022-07-18 [1] CRAN (R 4.1.2)
#> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.1.0)
#> proxy 0.4-27 2022-06-09 [1] CRAN (R 4.1.2)
#> purrr 0.3.4 2020-04-17 [1] CRAN (R 4.1.0)
#> R.cache 0.16.0 2022-07-21 [1] CRAN (R 4.1.2)
#> R.methodsS3 1.8.2 2022-06-13 [1] CRAN (R 4.1.2)
#> R.oo 1.25.0 2022-06-12 [1] CRAN (R 4.1.2)
#> R.utils 2.12.0 2022-06-28 [1] CRAN (R 4.1.2)
#> R6 2.5.1 2021-08-19 [1] CRAN (R 4.1.0)
#> Rcpp 1.0.9 2022-07-08 [1] CRAN (R 4.1.2)
#> reprex 2.0.1 2021-08-05 [1] CRAN (R 4.1.0)
#> rlang 1.0.4 2022-07-12 [1] CRAN (R 4.1.2)
#> rmarkdown 2.14 2022-04-25 [1] CRAN (R 4.1.2)
#> rstudioapi 0.13 2020-11-12 [1] CRAN (R 4.1.0)
#> sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.1.0)
#> sf * 1.0-7 2022-03-07 [1] CRAN (R 4.1.2)
#> stringi 1.7.6 2021-11-29 [1] CRAN (R 4.1.0)
#> stringr 1.4.0 2019-02-10 [1] CRAN (R 4.1.2)
#> styler 1.7.0 2022-03-13 [1] CRAN (R 4.1.2)
#> terra 1.5-21 2022-02-17 [1] CRAN (R 4.1.2)
#> tibble 3.1.8 2022-07-22 [1] CRAN (R 4.1.2)
#> tidyselect 1.1.2 2022-02-21 [1] CRAN (R 4.1.2)
#> units 0.8-0 2022-02-05 [1] CRAN (R 4.1.2)
#> utf8 1.2.2 2021-07-24 [1] CRAN (R 4.1.0)
#> vctrs 0.4.1 2022-04-13 [1] CRAN (R 4.1.2)
#> withr 2.5.0 2022-03-03 [1] CRAN (R 4.1.2)
#> xfun 0.31 2022-05-10 [1] CRAN (R 4.1.2)
#> XML 3.99-0.10 2022-06-09 [1] CRAN (R 4.1.2)
#> yaml 2.3.5 2022-02-21 [1] CRAN (R 4.1.2)
#>
#> [1] /Library/Frameworks/R.framework/Versions/4.1/Resources/library
#>
#> ──────────────────────────────────────────────────────────────────────────────Here's also the header of the downloaded coverage which agrees with the boundaries returned from boundedBy:
class : SpatRaster
dimensions : 6431, 7508, 1 (nrow, ncol, nlyr)
resolution : 1000, 1000 (x, y)
extent : -623000, 6885000, 604000, 7035000 (xmin, xmax, ymin, ymax)
coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035)
source : emodnet__2017_01_st_00.tif
name : emodnet__2017_01_st_00
eblondel commented
@annakrystalli i'm afraid the info is given in two distinct SRS and not EPSG:3857 for both:
- the bounding box retrived through
summary$BoundingBox$BoundingBox$getBBOX()is retrieved through the WCS Capabilities document, and indeed srs is EPSG:3857
Still if you do summary$getBoundingBox()$BoundingBox$attrs$crs you get the BoundingBox SRS associated to the coverage summary as defined in the capabilities: EPSG:3857
- the envelope retrieved through
summary$getDescription()$boundedByis got from the WCS coverage description at https://ows.emodnet-humanactivities.eu/wcs?service=WCS&version=2.0.1&coverageId=emodnet__2017_01_st_00&request=DescribeCoverage If you look at the envelope, we have a different srsName:
In R, if you type summary$getDescription()$boundedBy$attrs$srsName you will not get EPSG:3857