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.biv...@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.biv...@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]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
