search

eblondel / ows4R

R Interface for OGC Web-Services (OWS)
https://eblondel.github.io/ows4R/
Other
36 stars 8 forks source link

`OWSBoundingBox` - `getBBOX` returning incorrect crs #81

Closed annakrystalli closed 1 year ago

annakrystalli commented 2 years ago

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 2 years ago

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 2 years ago

Great, I'm back on the project tomorrow so will take a look and let you know!

annakrystalli commented 2 years ago

Just checked and all is working as expected!

annakrystalli commented 2 years ago

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 ``` r 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 1 year ago

@annakrystalli i'm afraid the info is given in two distinct SRS and not EPSG:3857 for both:

image

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

In R, if you type summary$getDescription()$boundedBy$attrs$srsName you will not get EPSG:3857