ArcGIS Runtime:计算不同空间参考中相同地图点之间的角度时得到不同的结果

ArcGIS Runtime: Getting different results when calculating angles between same map points in different spatial references

我正在使用 ArcGIS Runtime SDK for dotNet 10.2.6 进行开发。

我正在尝试计算地图上两点之间的方位角(与北方的度数)。

地图的空间参考为 3857 (WGS_1984_Web_Mercator_Auxiliary_Sphere),我使用空间参考 4326 (GCS_WGS_1984) 供用户显示。

为了使用户 POV 向北,我让用户单击他在地图上的位置以及他当前正在查看的地图上的另一个点。然后我使用简单的正切来计算两个给定点之间的方位角。

当我使用地图空间参考 (3857) 时,我得到了合乎逻辑的结果。 当我将这两个点投影到用户空间参考 (4326) 并进行相同的计算时,当方位角约为 45 度时(所有象限都相同,即 135、225 和 315),我得到约 5 度的误差。

示例:

Azimuth 0:
[spatial reference Wkid=3857] dAzimuth=00.066 
x1=3910430.46448709 y1=3847292.56056299 
x2=3910431.15764925 y2=3847888.00648926
[spatial reference Wkid=4326] dAzimuth=00.079 
x1=35.1279945373535 y1=32.6375465393066 
x2=35.1280007641351 y2=32.6420507989248

Azimuth 45:
[spatial reference Wkid=3857] dAzimuth=45.002 
x1=3910430.46448709 y1=3847292.56056299 
x2=3910850.25858266 y2=3847712.32276465
[spatial reference Wkid=4326] dAzimuth=49.901 (*** ~5 degrees error ***)
x1=35.1279945373535 y1=32.6375465393066 
x2=35.1317656118759 y2=32.6407218603483

示例中各点之间的距离约为 500 米。 角度计算为:

dAzimuth = Math.Atan(dDeltaX / dDeltaY) * RADIAN_TO_DEGREE_FACTOR;

投影是用几何引擎完成的:

projectedPoint = GeometryEngine.Project(inputPoint, targetSpatialReference) as MapPoint;

什么可能导致此错误?

"error"可能是由于预测的差异。地理坐标系 (4326) 和投影坐标系 (3857) 将给出不同的结果。

另请参阅:

EPSG:4326不是真正的投影,如果你把它看成投影,它肯定不是共形的。

所以你不应该依赖那个 "projection" 的角度,但你必须依赖地球 bearing formulas 来获得精确的结果。


所以,地球上两点之间的方位(纬度,经度)(R代码):

R> vectorBearing <- function(a_lat, a_lon, b_lat, b_lon) {
    a_phi    <- a_lat * pi / 180
    a_lambda <- a_lon * pi / 180
    b_phi    <- b_lat * pi / 180
    b_lambda <- b_lon * pi / 180
    dL <- b_lambda - a_lambda
    cosBLat <- cos(b_phi)
    Y <- cosBLat*sin(dL)
    X <- (cos(a_phi)*sin(b_phi))-(sin(a_phi)*cosBLat*cos(dL))
    atan2(Y, X) * 180 / pi                                     
}

在几何上,平面上 3 点 ABC 之间的角度(在 A 上):

R> angleABC <- function(x1, y1, x2, y2, x3, y3) {
    P12 <- sqrt((x1-x2)^2 + (y1-y2)^2)
    P13 <- sqrt((x1-x3)^2 + (y1-y3)^2)
    P23 <- sqrt((x2-x3)^2 + (y2-y3)^2)

    if (P12==0 || P23==0 || P13 ==0) return(0)
    acos((P12**2 + P13**2 - P23**2) / (2*P12*P13) ) / pi * 180
}

你的具体例子,在 latitude/longitude 点上使用轴承公式给出了正确的结果:

R> vectorBearing(32.6375465393066, 35.1279945373535, 32.6407218603483, 35.1317656118759)
[1] 45.00116

同样在投影点上使用几何公式:

R> x1<-3910430.46448709
R> y1<-3847292.56056299
R> x2<-3910430.46448709
R> y2<-4000000.00000000  # any
R> x3<-3910850.25858266
R> y3<-3847712.32276465
R> angleABC(x1,y1,x2,y2,x3,y3)
[1] 45.00218

(结果不完全相同,原因 1)可能是一些有限算术错误 2)EPSG:3857(Web Mercator)是 slightly non-conformal)。

然而,直接在 lat/lon 值上使用平面函数:

R> x1<-35.1279945373535
R> y1<-32.6375465393066
R> x2<-35.1317656118759
R> y2<-40.0000000000000  # any
R> x2<-35.1279945373535
R> x3<-35.1317656118759
R> y3<-32.6407218603483
R> angleABC(x1,y1,x2,y2,x3,y3)
[1] 49.90194

...给出您所看到的不精确的结果。


另请参阅:https://gis.stackexchange.com/a/21363/6998