I have some spatial data recorded by a GPS/GNSS that includes Lat/Lon and Height Above Ellipsoid (HAE, meters) altitude values, all are referenced to NAD83(2011), EPSG:6318. I'm trying to convert the horz. to State Plane NY East (EPSG: 6537, which is successful) but how do I convert the altitude/height data to NAVD88, US-ft (EPSG: 6360)? I can apply the horizontal using sf::st_transform(6537) but I don't know how to apply the vertical. I tried just stringing another sf::st_transform(6360) after the first transform but that didn't work. Below is a reproducible example and for reference, the output from NOAA VDatum online that shows the correct value (transformed elevation should be ~5.899 ft)
library(sf)
library(tidyverse)
pt <- data.frame(Lat = 41.1578110483, Lon = -73.8716163883, Alt = -29.3619984455) %>%
st_as_sf(coords = c("Lon", "Lat", "Alt"), crs = 6318, agr = "constant", remove = FALSE) %>%
st_transform(6537) %>%
mutate(x = st_coordinates(.)[,1],
y = st_coordinates(.)[,2],
z = st_coordinates(.)[,3])
pt
Simple feature collection with 1 feature and 6 fields
Attribute-geometry relationship: 3 constant, 0 aggregate, 0 identity, 3 NA's
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: 665148.732436 ymin: 847314.319632 xmax: 665148.732436 ymax: 847314.319632
z_range: zmin: -29.3619984455 zmax: -29.3619984455
Projected CRS: NAD83(2011) / New York East (ftUS)
Lat Lon Alt geometry x y z
1 41.1578110483 -73.8716163883 -29.3619984455 POINT Z (665148.732436 8473... 665148.732436 847314.319632 -29.3619984455
Update:
sf_extSoftVersion()
GEOS GDAL proj.4 GDAL_with_GEOS USE_PROJ_H PROJ
"3.9.0" "3.2.1" "7.2.1" "true" "true" "7.2.1"
