查找 R 中的投影与 ArcMap 中的投影的差异
Finding differences in projections in R and projections in ArcMap
我当前的问题是创建一个 R 程序或函数，它创建一个从纬度和经度点到 NAD-1983 2011 StatePlane Pennsylvania South FIPS - 3702 数据的数据帧投影。目前，我遇到了障碍。这是我一直在使用的输入数据：
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
我一直在尝试的当前函数使用 rgdal 和 sf 函数。我正在尝试执行以下操作：
- 在lat/long
中定义data_x数据投影
- 将数据投影到 NAD-1983 2011 StatePlane Pennsylvania South FIPS - 3702 数据
StatePlane 信息如下：
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
这是我想出的代码
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)
}
这是我发现的一些输出结果的示例。
来自 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
来自 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
虽然我没有遇到任何破译错误，但这些结果是距离不一致的因素。有没有我遗漏的包裹？我是否缺少命令或步骤？
AM 在 运行 处理你的代码时遇到问题，因为如果我传递一个数据框，它会出错
Error in UseMethod("st_crs<-") :
no applicable method for 'st_crs<-' applied to an object of class "data.frame"
>
我想您一定是在以某种方式将其读入空间数据框中。此外，您的目标投影字符串似乎已损坏。我将向您展示我将如何做到这一点：
首先，使用这样的文本文件中的数据（long-lat）：
-75.35843 40.12232
-75.36189 40.12347
-75.36404 40.12456
-75.37228 40.1287
我读了它：
> 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
然后使用sf包函数将其变成EPSG:4326坐标系的点——常用GPS坐标：
> 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)
现在我可以转换到任何其他坐标系。您的目标似乎是 EPSG:6564 (WKID: 6564 Authority: EPSG) (http://epsg.io/6564) 尽管您还提到了 3702，这似乎与怀俄明州有关...)。
> 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)
>
并且这些点与您为 R 代码或 ARCGis 提供的点不匹配，假设您提供的 lat-longs 应该对应于您提供的转换后的 X-Y 坐标。
我非常有信心地说我上面给出的数字是源中 EPSG:4326 lat-long 坐标的 EPSG:6564 坐标。
我还设法修复了损坏的投影字符串：
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"
产生与 EPSG:6564 相同的结果。
无法 运行 您的代码并重现您的结果，我不确定问题出在哪里，但我的代码就是我这样做的方式，我非常有信心它是正确的。
我当前的问题是创建一个 R 程序或函数，它创建一个从纬度和经度点到 NAD-1983 2011 StatePlane Pennsylvania South FIPS - 3702 数据的数据帧投影。目前，我遇到了障碍。这是我一直在使用的输入数据：
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
我一直在尝试的当前函数使用 rgdal 和 sf 函数。我正在尝试执行以下操作：
- 在lat/long 中定义data_x数据投影
- 将数据投影到 NAD-1983 2011 StatePlane Pennsylvania South FIPS - 3702 数据
StatePlane 信息如下：
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
这是我想出的代码
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)
}
这是我发现的一些输出结果的示例。
来自 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
来自 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
虽然我没有遇到任何破译错误，但这些结果是距离不一致的因素。有没有我遗漏的包裹？我是否缺少命令或步骤？
AM 在 运行 处理你的代码时遇到问题，因为如果我传递一个数据框，它会出错
Error in UseMethod("st_crs<-") :
no applicable method for 'st_crs<-' applied to an object of class "data.frame"
>
我想您一定是在以某种方式将其读入空间数据框中。此外，您的目标投影字符串似乎已损坏。我将向您展示我将如何做到这一点：
首先，使用这样的文本文件中的数据（long-lat）：
-75.35843 40.12232
-75.36189 40.12347
-75.36404 40.12456
-75.37228 40.1287
我读了它：
> 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
然后使用sf包函数将其变成EPSG:4326坐标系的点——常用GPS坐标：
> 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)
现在我可以转换到任何其他坐标系。您的目标似乎是 EPSG:6564 (WKID: 6564 Authority: EPSG) (http://epsg.io/6564) 尽管您还提到了 3702，这似乎与怀俄明州有关...)。
> 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)
>
并且这些点与您为 R 代码或 ARCGis 提供的点不匹配，假设您提供的 lat-longs 应该对应于您提供的转换后的 X-Y 坐标。
我非常有信心地说我上面给出的数字是源中 EPSG:4326 lat-long 坐标的 EPSG:6564 坐标。
我还设法修复了损坏的投影字符串：
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"
产生与 EPSG:6564 相同的结果。
无法 运行 您的代码并重现您的结果，我不确定问题出在哪里，但我的代码就是我这样做的方式，我非常有信心它是正确的。