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

Fri Jul 29 16:56:47 CEST 2022

On Thu, 28 Jul 2022, inti luna wrote:

> 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.

Good, thanks.

It appears that terra does not accept POINT Z geometries:

> library(sf)
Linking to GEOS 3.11.0, GDAL 3.5.1, PROJ 9.0.1; sf_use_s2() is TRUE
> data.points <- st_as_sf(data1, coords = c("Longitude", "Latitude",
+    "Ellipsoidal_height"))
> dp_sv <- vect(data.points)
> geom(dp_sv, wkt=TRUE)
[1] "POINT (2.321528 41.757588)" "POINT (2.321299 41.757575)"
[3] "POINT (2.321216 41.7576)"   "POINT (2.321659 41.75765)"
[5] "POINT (2.321639 41.757646)" "POINT (2.321674 41.75794)"
[7] "POINT (2.321186 41.7577)"   "POINT (2.321153 41.757706)"
[9] "POINT (2.320972 41.757793)"
> st_geometry(data.points)
Geometry set for 9 features
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 2.320972 ymin: 41.75757 xmax: 2.321674 ymax: 41.75794
CRS:           NA
First 5 geometries:
POINT Z (2.321528 41.75759 1235.577)
POINT Z (2.321299 41.75757 1234.38)
POINT Z (2.321216 41.7576 1233.49)
POINT Z (2.321659 41.75765 1235.785)
POINT Z (2.321639 41.75765 1235.741)

so exploring vertical grids is not productive for now. This is worth 
addressing, to provide someting to check sf with.

Roger

>
> 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
>
>
>
>

-- 
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