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

[inti luna]

Dear Professor Roger,
Thank you so much for the clarification. I will try to contact others
interested in vertical transformations too and make a post.

Kind regards,
Inti Luna

El jue, 28 jul 2022 a las 13:38, Roger Bivand (<Roger.Bivand using nhh.no>)
escribió:

> 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



-- 
+ 34 6042 19917
Skype: inti_luna
Zoom *ID* : 901-077-6684
*https://evolors.com <https://evolors.com/>*

	[[alternative HTML version deleted]]