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


#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"))

#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

ESRI Project tool


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")
| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.