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:
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?
AM having trouble running your code because if I pass a data frame it errors with
Error in UseMethod("st_crs<-") :
no applicable method for 'st_crs<-' applied to an object of class "data.frame"
>
I guess you must be reading it into a spatial data frame somehow. Also your target projection string seems broken. I'll show you how I'd do this:
First, with the data in a text file like this (which is long-lat):
-75.35843 40.12232
-75.36189 40.12347
-75.36404 40.12456
-75.37228 40.1287
I'd read it in with:
> pts = read.table("latlong.txt",head=FALSE)
> head(pts)
V1 V2
1 -75.35843 40.12232
2 -75.36189 40.12347
3 -75.36404 40.12456
4 -75.37228 40.12870
5 -75.37835 40.13243
6 -75.38479 40.13835
Then use the sf package functions to make it into points with the EPSG:4326 coordinate system - the common GPS coordinates:
> pts2 = st_as_sf(pts, coords=c(1,2), crs=4326)
> pts2
Simple feature collection with 14 features and 0 fields
geometry type: POINT
dimension: XY
bbox: xmin: -75.43579 ymin: 40.12232 xmax: -75.35843 ymax: 40.17576
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
First 10 features:
geometry
1 POINT (-75.35843 40.12232)
2 POINT (-75.36189 40.12347)
3 POINT (-75.36404 40.12456)
4 POINT (-75.37228 40.1287)
5 POINT (-75.37835 40.13243)
6 POINT (-75.38479 40.13835)
7 POINT (-75.38961 40.14198)
8 POINT (-75.39536 40.14724)
9 POINT (-75.40018 40.15457)
10 POINT (-75.40614 40.15849)
Now I can transform to any other coordinate system. Your target appears to be EPSG:6564 (WKID: 6564 Authority: EPSG) (http://epsg.io/6564) although you also mention 3702 which is seems to be related to Wyoming...).
> st_transform(pts2, 6564)
Simple feature collection with 14 features and 0 fields
geometry type: POINT
dimension: XY
bbox: xmin: 797083.4 ymin: 90364.93 xmax: 803830.7 ymax: 96120.91
epsg (SRID): 6564
proj4string: +proj=lcc +lat_1=40.96666666666667 +lat_2=39.93333333333333 +lat_0=39.33333333333334 +lon_0=-77.75 +x_0=600000 +y_0=0 +ellps=GRS80 +units=m +no_defs
First 10 features:
geometry
1 POINT (803830.7 90364.93)
2 POINT (803532.4 90484.59)
3 POINT (803345.9 90600.62)
4 POINT (802631.5 91041.2)
5 POINT (802103.2 91441.3)
6 POINT (801536.9 92083.67)
7 POINT (801115.5 92475.59)
8 POINT (800610.2 93046.34)
9 POINT (800177.9 93849)
10 POINT (799658.8 94270.61)
>
and those points don't match the ones you gave for either your R code or ARCGis, assuming the lat-longs you gave are supposed to be corresponding to the transformed X-Y coords you gave.
I'd be very confident in saying that the numbers I've given above are the EPSG:6564 coordinates of the EPSG:4326 lat-long coordinates in the source.
Also I managed to fix your broken projection string:
proj <- "+lat_1=39.93333333333333 +lat_2=40.96666666666667 +lat_0=39.33333333333334 +lon_0=-77.75 +x_0=600000.0 +y_0=0 +proj=lcc +ellps=GRS80 +datum=NAD83 +units=m"
and that produces the same results as EPSG:6564.
Not being able to run your code and reproduce your results, I'm not sure where the problem lies, but my code is how I'd do this and I'm pretty confident its right.