Message-ID: <a24b11a9-d531-88cf-7234-2b0d18855c7@reclus.nhh.no>
Date: 2021-03-26T14:47:01Z
From: Roger Bivand
Subject: Not perfect match between the bounding box for a EPSG when projected (st_bbox from sf)
In-Reply-To: <CABkCotSWqOf_uT31DmoVyKixbg0qFL6M+O7KYo-A+mWaMbnOvQ@mail.gmail.com>
On Thu, 25 Mar 2021, Manuel Sp?nola wrote:
> Dear list members,
>
> I am working with the epsg 5367.
>
> The information for this epsg according to this website (
> https://epsg.io/5367) is:
>
> WGS84 bounds:
> -90.45 2.15
> -81.43 11.77
>
> Projected bounds:
> -218963.01 239235.14
> 780173.56 1302793.77
>
> When I create a box with st_bbox from the package sf and transform the
> object I don't get the same bounding box that I was expecting. Is it a
> fair expectation?
No. Firstly, epsg.io is not an EPSG service. Secondly, two different
definitions are offered. I think that this works (sf 0.9-8):
pts1 <- st_sfc(list(st_point(c(-90.45, 2.15)), st_point(c(-81.43,
11.77))), crs="OGC:CRS84")
st_transform(pts1, "EPSG:5367")
# POINT (-218962 239233.5)
# POINT (780173.7 1302793)
o <- sf_proj_pipelines("OGC:CRS84", "EPSG:5367")
o$definition
st_transform(pts1, "EPSG:5367", pipeline = o$definition[1])
# POINT (-218962 239233.5)
# POINT (780173.7 1302793)
# Costa Rica - onshore and offshore, accuracy 1.0 m, code 8969 [grid]
st_transform(pts1, "EPSG:5367", pipeline = o$definition[2])
# POINT (-218963 239235.1)
# POINT (780173.6 1302794)
# Costa Rica - onshore and offshore, accuracy 1.0 m, code 5376 (default)
# [grid]
So, as in https://github.com/r-spatial/sf/issues/1634, treating epsg.io as
authoritative (especially when there are different transformation
pipelines in play) may be unhelpful.sf/PROJ chooses a different pipeline
as default than epsg.io.
This roundtrips:
st_transform(st_transform(pts1, "EPSG:5367"), "OGC:CRS84")
# POINT (-90.45 2.15)
# POINT (-81.43 11.77)
# Costa Rica - onshore and offshore, accuracy 1.0 m, code 8969 [grid]
but for epsg.io's default, one needs to intervene manually to choose
pipelines both ways that are not sf's default choices.
Of course, here I am not using the bounding box polygon, because its
projected version is not a rectangle, so its bounding box expands the
range, especially further from the equator.
Hope this clarifies,
Roger
>
> b <- st_bbox(c(xmin = -90.45, xmax = -81.43, ymax = 11.77, ymin = 2.15),
> crs = st_crs(4326))
>
> b <- st_as_sfc(b)
>
> b_5367 <- st_transform(b, crs = 5367)
>
> b_5367
>
> Geometry set for 1 feature
> Geometry type: POLYGON
> Dimension: XY
> Bounding box: xmin: -218962 ymin: 237951.9 xmax: 785959.7 ymax: 1309622
> Projected CRS: CR05 / CRTM05POLYGON ((-218962 239233.5, 785959.7 237951.9, ...
>
>
>
>
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway.
e-mail: Roger.Bivand at nhh.no
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en