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