[R-sig-Geo] Height transformation (ellipsoidal to EGM2008)

[Roger Bivand]

On Tue, 26 Jul 2022, inti luna wrote:

> Hi all,
> I have a table with data points with CRS EPSG 4326 and ellipsoidal height
> (m) and would like to transform it to another CRS. Target horizontal is
> EPSG: 25831 and target height is based on EGM2008 (EPSG 3855).
>
> I have carried the horizontal transformation with the sf package but I
> don´t know how to perform the height transformation. Any guidance is
> appreciated. I have been checking
> https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fkeen-swartz-3146c4.netlify.app%2Fsf.html%23st_transform-sf_project&data=05%7C01%7CRoger.Bivand%40nhh.no%7Cccb34858b29945e7edb508da6f71f01d%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C637944834193544665%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C2000%7C%7C%7C&sdata=tSwZgNFNWpnn5%2BeFWpyUUDSxws7QGYbSi%2Bn3G7GywkI%3D&reserved=0 but
> still have not found information about height transformation

No, this is much more complex. Please read

https://proj.org/apps/cs2cs.html

to begin with (not self-explanatory but accurate). R does not provide PROJ 
apps, but running on my Linux machine (PROJ 9.0.1, but 7.2.1 is fine):

$ export PROJ_NETWORK=ON
$ echo 2.32152788 41.75758799 1235.577 | cs2cs EPSG:4326 
EPSG:25831+EPSG:3855
5182922.95	329657.82 1264.66
$ export PROJ_NETWORK=OFF
$ echo 2.32152788 41.75758799 1235.577 | cs2cs EPSG:4326 
EPSG:25831+EPSG:3855
5182922.95	329657.82 1235.58

When the CDN network is ON, us_nga_egm08_25.tif is available for grid 
transformation in 3D. When off, it is not:

$ projinfo -o PROJ -s EPSG:4326 -t EPSG:25831+EPSG:3855

gives a lot of output  describing all of the transformation operations 
that might be possible if known grids are available.

Starting R with the CDN ON:

$ PROJ_NETWORK=ON R
...
library(sf)
o <- sf_proj_pipelines("EPSG:4326", "EPSG:25831+EPSG:3855")
o

o[1,] is the best than can be instantiated with:

+proj=pipeline
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +inv +proj=vgridshift +grids=us_nga_egm08_25.tif +multiplier=1
  +step +proj=utm +zone=31 +ellps=GRS80

so now making a vertical shift using a grid from the CDN.

Creating your data.points (thanks for a good reprex!) but naming 
"Ellipsoidal_height" as the third coordinate:

data.points <- st_as_sf(data1, coords = c("Longitude", "Latitude",
   "Ellipsoidal_height"))

I get:

data.points2<-st_transform(data.points, "EPSG:25831")
data.points2$geometry
First 5 geometries:
POINT Z (443597.6 4623084 1235.577)

which is wierd, but using the pipeline found above:

data.points2a<-st_transform(data.points, "EPSG:25831",
  pipeline=o$definition[1])
data.points2a$geometry

giving

POINT Z (5182923 329657.8 1264.658)

as in the PROJ app.

One point is that you need either to install the grid locally, or enable 
the CDN *before* you start R and load sf.

Another is that sf_proj_pipelines() uses PROJ directly, so can identify 
combined pipelines, but st_transform() uses PROJ via GDAL, and doesn't 
handle this correctlly, so one must pass the pipeline. The output CRS will 
also be wrong, it should be:

COMPOUNDCRS["ETRS89 / UTM zone 31N + EGM2008 height", ...

not

PROJCRS["ETRS89 / UTM zone 31N", ... .

Finally, your sf object needed the third coordinate defined, otherwise it 
had no idea that it wasn't 2D.

Consider creating an sf issue, and if possible contact others needing 
vertical transformation and consider contributing a blog or vignette. Both 
sf and terra need user input to document important topics like this.

Best wishes,

Roger


>
> Repex
>
> library(sf)
>
> #create sample dataset
>
> data1<-"Name,Longitude,Latitude,Ellipsoidal_height
> gcp1,2.32152788,41.75758799,1235.577
> gcp2,2.32129905,41.75757489,1234.38
> gcp3,2.3212163,41.75760007,1233.49
> gcp4,2.32165884,41.75765011,1235.785
> gcp5,2.32163868,41.75764588,1235.741
> gcp6,2.32167394,41.75793979,1236.763
> gcp7,2.32118573,41.75769993,1232.364
> gcp8,2.32115307,41.75770638,1231.897
> gcp9,2.32097163,41.75779261,1228.454"
>
> data1<-read.table(text = data1, header = T,sep = ",")
> head(data1)
> str(data1)
>
> sf::sf_extSoftVersion()
>
> crs_epsg25831 <- st_crs(25831)
> crs_epsg25831
>
> data.points <- st_as_sf(data1, coords = c("Longitude", "Latitude"))
> st_crs(data.points) <- 4326
> data.points
>
> #transform from EPSG 4326 to EPSG 25831(ETRS89/UTM 31N)
> data.points2<-st_transform(data.points,crs_epsg25831)
>
> st_crs(data.points2)
> data.points2
> data.points2$geometry
> write.csv(data.points2, "sample.points_epsg25831.csv")
>
> # height transformation ?
>
>
> Machine and software details:
>
> OS: Windows
> R: version 4.2.1 (2022-06-23 ucrt) -- "Funny-Looking Kid"
> Rstudio: RStudio 2022.07.1+554 "Spotted Wakerobin" Release
> (7872775ebddc40635780ca1ed238934c3345c5de, 2022-07-22) for Windows
> Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like
> Gecko) QtWebEngine/5.12.8 Chrome/69.0.3497.128 Safari/537.36
>
> GEOS: 3.9.1
> GDAL 3.4.3
> PROJ 7.2.1
> sf_use_s2: TRUE
>
>
>
> Kind regards,
> Inti Luna
>
>

-- 
Roger Bivand
Emeritus Professor
Department of Economics, Norwegian School of Economics,
Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway.
e-mail: Roger.Bivand using nhh.no
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en