Using specific transformation in R sf::st_transform() to reproduce ESRI’s Project tool

Question

I am reading in a CSV with latitude/longitude points (WGS 84, EPSG 4326), transforming into the NAD_1983_HARN_StatePlane_Maryland_FIPS_1900_Feet projection (WKID 2893) and saving as a feature class. I wrote an R script to automate work that was being done in ESRI ArcMap.
When performing this operation in ArcMap using the Project tool, the default Geographic Transformation is NAD_1983_HARN_To_WGS_1984_2 (WKID 1900, accuracy 0.10 m).
When I use sf::st_transform(), the results do not match.  R appears to be using the transformation NAD_1983_HARN_To_WGS_1984 (WKID 1580, accuracy 1.0 m). It results in my points being off by ~2.88 feet.
Is there a way to force st_transform() to use the more accurate transformation?
I can't find a way to reference a specific transformation with this package. I even tried building a custom pipeline to pass to st_transform(), but am having trouble understanding the syntax and order of operations.
Here are the State Plane coordinates from my two example data points that show the discrepancy:
R Data:
POINT (1454828 413134.1)
POINT (1454429 485972.8)

ESRI Data:
POINT (1454829 413131.3)
POINT (1454430 485970)

And example code showing how I'm doing the transformation in R
library(sf)
library(arcgisbinding)
arc.check_product()

#create example data and save to CSV
df = data.frame(X_COORD=c(-76.5,-76.5),Y_COORD=c(38.8,39),DESC=c("Pt_A","Pt_B"))
write.csv(df,"df.csv")

#Read CSV as sf object and define the CRS
R_WGS84 <- st_read("df.csv",options=c("X_POSSIBLE_NAMES=X_COORD","Y_POSSIBLE_NAMES=Y_COORD"))
st_crs(R_WGS84) <- 4326

#Transform the data into State Plane projection. This is where things deviate from ArcMap.
R_StatePlane <- st_transform(R_WGS84,2893)

#Write feature class to file GDB
arc.write("PathToExampleGDB.gdb/R_StatePlane"), R_StatePlane,shape_info=list(type='Point',WKID=2893),overwrite = TRUE)

And an example of the "Project" Tool in ArcMap

Brennan · Accepted Answer

Use rgdal::list_coordOps() to get more information about what transformations are available.
> list_coordOps("EPSG:4326", "EPSG:2893")
Candidate coordinate operations found:  2 
Strict containment:  FALSE 
Visualization order:  TRUE 
Source: EPSG:4326 
Target: EPSG:2893 
Best instantiable operation has accuracy: 1 m
Description: axis order change (2D) + Inverse of NAD83(HARN) to WGS 84 (1) + SPCS83 Maryland zone (US Survey feet) 
Definition:  +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=lcc +lat_0=37.6666666666667 +lon_0=-77 +lat_1=39.45 +lat_2=38.3
             +x_0=399999.9998984 +y_0=0 +ellps=GRS80 +step +proj=unitconvert +xy_in=m +xy_out=us-ft

What is shown above corresponds to the less accurate NAD_1983_HARN_To_WGS_1984 method. The other candidate operation corresponds to NAD_1983_HARN_To_WGS_1984_3.
> list_coordOps("EPSG:4326", "EPSG:2893")$definition[2]
[1] "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 +step +proj=cart +ellps=WGS84 +step +inv +proj=helmert +x=-0.991 +y=1.9072 +z=0.5129 +rx=-0.0257899075194932 +ry=-0.0096500989602704 +rz=-0.0116599432323421 +s=0 +convention=coordinate_frame +step +inv +proj=cart +ellps=GRS80 +step +proj=pop +v_3 +step +proj=lcc +lat_0=37.6666666666667 +lon_0=-77 +lat_1=39.45 +lat_2=38.3 +x_0=399999.9998984 +y_0=0 +ellps=GRS80 +step +proj=unitconvert +xy_in=m +xy_out=us-ft"

So take that +proj string, add in the axisswaps that were missing for some reason, and update the helmert parameters from NAD_1983_HARN_To_WGS_1984_2.  The results match the ESRI output within 0.0003 feet.
R_StatePlane <-st_transform(R_WGS84, pipeline ="+proj=pipeline 
                      +step +proj=axisswap +order=2,1 
                      +step +proj=unitconvert +xy_in=deg +xy_out=rad 
                      +step +proj=push +v_3 
                      +step +proj=cart +ellps=WGS84 
                      +step +inv +proj=helmert +x=-0.9738 +y=1.9453 +z=0.5486 +rx=-0.02755079017042466 +ry=-0.01004922136035853 +rz=-0.0113590028800276 +s=0 +convention=coordinate_frame 
                      +step +inv +proj=cart +ellps=GRS80 
                      +step +proj=pop +v_3 
                      +step +proj=lcc +lat_0=37.6666666666667 +lon_0=-77 +lat_1=39.45 +lat_2=38.3 +x_0=399999.9998984 +y_0=0 +ellps=GRS80 
                      +step +proj=unitconvert +xy_in=m +xy_out=us-ft
                      +step +proj=axisswap +order=2,1")

Using specific transformation in R sf::st_transform() to reproduce ESRI’s Project tool

One Answer

Add your own answers!

Ask a Question