My current issue is creating an R program or function which creates a dataframe projection from latitude and longitude points to NAD-1983 2011 StatePlane Pennsylvania South FIPS - 3702 data. Currently, I've hit a snag. Here is the input data I've been using:
Latitude/Longitude:
-75.35843 40.12232
-75.36189 40.12347
-75.36404 40.12456
-75.37228 40.1287
-75.37835 40.13243
-75.38479 40.13835
-75.38961 40.14198
-75.39536 40.14724
-75.40018 40.15457
-75.40614 40.15849
-75.40939 40.16014
-75.41906 40.16572
-75.43034 40.17133
-75.43579 40.17576
The current function I've been trying uses rgdal and sf functions. I'm trying to do the following:
- Define the data_x data projection in lat/long
- Project the data to NAD-1983 2011 StatePlane Pennsylvania South FIPS - 3702 data
StatePlane info below:
NAD_1983_2011_StatePlane_Pennsylvania_South_FIPS_3702
WKID: 6564 Authority: EPSG
Projection: Lambert_Conformal_Conic
False_Easting: 600000.0
False_Northing: 0.0
Central_Meridian: -77.75
Standard_Parallel_1: 39.93333333333333
Standard_Parallel_2: 40.96666666666667
Latitude_Of_Origin: 39.33333333333334
Linear Unit: Meter (1.0)
Geographic Coordinate System: GCS_NAD_1983_2011
Angular Unit: Degree (0.0174532925199433)
Prime Meridian: Greenwich (0.0)
Datum: D_NAD_1983_2011
Spheroid: GRS_1980
Semimajor Axis: 6378137.0
Semiminor Axis: 6356752.314140356
Inverse Flattening: 298.257222101
And here's the code I came up with
pckgs <- c("rgdal", "sf")
install.packages(pckgs)
library(rgdal)
library(sf)
project_coordinates <- function(df="data_x", x="Long", y="Lat"){
# returns a dataframe with two new columns containing projected coordinates to NAD83
## df(str): dataframe to use
## x (str): name of column containing Longitude in degrees
## y (str): name of column containing Latitude in degrees
df <- st_set_crs(df, "+proj=longlat +ellps=GRS80 +datum=NAD83")
proj <- "+lat_1=39.93333333333333 +lat_2=40.96666666666667 +lat_0=39.33333333333334 +lon_0=-77.75 +x_0= 600000.0 +y_0=0 +k_0=lcc +ellps=GRS80 +datum=NAD83 +units=m no_defs"
lat_long_proj <- project(as.matrix(cbind(df[x],df[y])),proj)
df["Long_proj"] <-lat_long_proj[,1]
df["Lat_proj"] <- lat_long_proj[,2]
return(df)
}
Here is an example of some of the results from what output I've found.
From R:
2637234.459 296472.2704
2636255.877 296864.8632
2635644.121 297245.5394
2633300.071 298690.992
2631566.848 300003.65
2629709.03 302111.1583
2628326.55 303396.9907
2626668.485 305269.5429
2625250.476 307902.9271
2623547.213 309286.1584
2622623.208 309862.9293
2619867.744 311823.4605
2616662.666 313783.4043
2615097.884 315356.6863
And from ArcGIS:
103492.4449 778962.9451
103492.4652 778963.7957
103492.4652 778963.7957
103728.4336 778891.792
103785.0208 778469.1788
103675.0824 778236.0495
104262.4265 777852.7858
104271.2263 777849.1735
104276.7765 777849.042
104276.8168 777850.743
104277.9268 777850.7168
104279.057 777851.541
104279.057 777851.541
104280.1671 777851.5147
While I haven't encountered any code-breaking errors, these results are an inconsistent factor of distance apart. Is there a package I'm missing? Am I missing a command or a step?