查找 R 中的投影和 ArcMap 中的投影的差异答案

【问题标题】：Finding differences in projections in R and projections in ArcMap查找 R 中的投影和 ArcMap 中的投影的差异
【发布时间】：2019-05-29 14:55:36
【问题描述】：

我当前的问题是创建一个 R 程序或函数，它创建从纬度和经度点到 NAD-1983 2011 StatePlane Pennsylvania South FIPS - 3702 数据的数据帧投影。目前，我遇到了一个障碍。这是我一直在使用的输入数据：

纬度/经度：

-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 函数。我正在尝试执行以下操作：

以纬度/经度定义 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

虽然我没有遇到任何代码破解错误，但这些结果是距离不一致的因素。有我丢失的包裹吗？我是否缺少命令或步骤？

【问题讨论】：

标签： r projection sf map-projections arcmap

【解决方案1】：

我在运行您的代码时遇到问题，因为如果我传递一个数据框，它会出错

Error in UseMethod("st_crs<-") : 
  no applicable method for 'st_crs<-' applied to an object of class "data.frame"
>

我猜你一定是以某种方式将它读入空间数据框。此外，您的目标投影字符串似乎已损坏。我会告诉你我会怎么做：

首先，将数据放在这样的文本文件中（长纬度）：

-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 提供的点不匹配，假设您提供的经纬度应该与您提供的转换后的 X-Y 坐标相对应。

我非常有信心地说我上面给出的数字是源中 EPSG:4326 经纬度坐标的 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 相同的结果。

无法运行您的代码并重现您的结果，我不确定问题出在哪里，但我的代码就是我这样做的方式，我非常相信它是正确的。

【讨论】：