¿Cómo calcula ArcGIS la distancia entre dos puntos con una proyección no equidistante?

10

Esta es una pregunta complementaria a la anterior. ¿Puede sugerir algunos textos introductorios bien escritos sobre proyecciones de sistemas de coordenadas?


Supongamos que estoy trabajando con la proyección del mapa CH1903, que por lo que sé es conforme, pero no equidistante. El significado, los ángulos (forma) se han conservado, pero no las áreas, distancias o escalas. (Al menos estos no se han conservado exactamente ). Hasta ahora tan bueno.

Me pregunto qué tipo de cálculo realiza ArcGIS cuando ahora quiero calcular la distancia entre dos puntos. En ArcObjects, podría usar la IProximityOperatorinterfaz de la siguiente manera:

IPoint a = ...,
       b = ...;

double distance = ((IProximityOperator)a).ReturnDistance(b);

Pregunta: Cuando estoy trabajando con un sistema de referencia que no preserva con precisión las distancias, ¿qué haría ArcGIS cuando le pregunto la distancia entre dos puntos (como se muestra arriba)?

  • ¿Simplemente hace algunas matemáticas de Pitágoras (a 2 + b 2 = c 2 ) para obtener la distancia, lo que significa que la distancia devuelta solo será tan precisa como lo permita la proyección?

  • ¿O hará algo más complicado, como alguna forma de reproyección, para obtener una distancia más precisa?

( La misma pregunta, pero en general: una vez que se han proyectado esas geometrías, ¿realiza ArcGIS todos los cálculos simplemente en el espacio euclidiano, o la proyección del mapa utilizado todavía influye en los cálculos de distancias, ángulos, áreas, etc.?)

stakx
fuente
2
Cree una nueva pregunta en lugar de modificar su original. De lo contrario, subvierte todos los mecanismos en este sitio: ¿qué significan las calificaciones cuando hay dos o más preguntas en juego en un hilo? ¿Qué significaría marcar una respuesta como correcta? Etc.
whuber
1
@whuber: Aunque todo lo que se ha escrito en este hilo todavía está en el tema WRT, la pregunta original se hizo, estoy de acuerdo en que ahora hay realmente dos preguntas aquí. Es demasiado tarde para cambiar eso ahora, pero tendremos en cuenta sus consejos para la próxima vez.
stakx

Respuestas:

10

Si desea un método estable para calcular distancias geodésicas, le recomiendo el contenedor de Richie Carmichael para el motor de proyección de ESRI .

Actualización: Acabo de probar el código de Richie con ArcGIS 10.0 en Vista64 y obtengo una excepción después de llamar LoadLibrary. Lo investigaré más tarde.

Por ahora, sin embargo, aquí hay un código en respuesta a las preguntas en los comentarios de otra respuesta.

El código compara IProximityOperator para puntos con y sin referencias espaciales. Luego muestra cómo usar una proyección equidistante azimutal (con el primer punto como punto de tangencia) para encontrar la distancia del gran círculo.

private void Test()
{
    IPoint p1 = new PointClass();
    p1.PutCoords(-98.0, 28.0);

    IPoint p2 = new PointClass();
    p2.PutCoords(-78.0, 28.0);

    Debug.Print("Euclidian Distance {0}", EuclidianDistance(p1, p2));
    Debug.Print("Distance with no spatialref {0}", GetDistance(p1, p2));

    ISpatialReferenceFactory srf = new SpatialReferenceEnvironmentClass();
    IGeographicCoordinateSystem gcs =
    srf.CreateGeographicCoordinateSystem((int)esriSRGeoCSType.esriSRGeoCS_WGS1984);

    p1.SpatialReference = gcs;
    p2.SpatialReference = gcs;

    Debug.Print("Distance with spatialref {0}", GetDistance(p1, p2));
    Debug.Print("Great Circle Distance {0}", GreatCircleDist(p1, p2));

}
private double GetDistance(IPoint p1, IPoint p2)
{
    return ((IProximityOperator)p1).ReturnDistance(p2);
}

private double EuclidianDistance(IPoint p1, IPoint p2)
{
    return Math.Sqrt(Math.Pow((p2.X - p1.X),2.0) + Math.Pow((p2.Y - p1.Y), 2.0));
}

private double GreatCircleDist(IPoint p1, IPoint p2)
{
    ISpatialReferenceFactory srf = new SpatialReferenceEnvironmentClass();
    IProjectedCoordinateSystem pcs =
    srf.CreateProjectedCoordinateSystem((int)esriSRProjCSType.esriSRProjCS_WGS1984N_PoleAziEqui);
    pcs.set_CentralMeridian(true, p1.X);
    ((IProjectedCoordinateSystem2)pcs).LatitudeOfOrigin = p1.Y;
    p1.SpatialReference = pcs.GeographicCoordinateSystem;
    p1.Project(pcs);
    p2.SpatialReference = pcs.GeographicCoordinateSystem;
    p2.Project(pcs);
    return EuclidianDistance(p1, p2);
}

Aquí está la salida:

Euclidian Distance 20
Distance with no spatialref 20
Distance with spatialref 20
Great Circle Distance 1965015.61318737

Creo que sería interesante probar esto contra el motor de proyección dll (pe.dll). Publicaré resultados si alguna vez hago que el código de Richie funcione.

Actualización: una vez que cambié el código de Richies para compilar para x86, lo ejecuté. Interesante ... la gran distancia del círculo que me da es 1960273.80162999, una diferencia significativa con respecto al método equidistante azimutal anterior.

Kirk Kuykendall
fuente
La razón de la discrepancia es probable porque el segmento de línea (en el PCS) que conecta los puntos no es una geodésica proyectada, que sería curvilínea cuando se proyecta. En consecuencia, obtienes un valor menor del que deberías. Una prueba de esta teoría es simple de hacer: tome una geodésica simple (como el ecuador) y compare dos cálculos de la distancia entre dos puntos ampliamente separados en la geodésica. Uno es el cálculo directo, como en su código; el otro divide la geodésica en segmentos, calcula directamente las longitudes de los segmentos y los suma. Este último debería ser más preciso.
whuber
9

En ArcGIS 10, vea IGeometryServer2 que ahora tiene GetDistanceGeodesic (distancia geodésica entre dos geometrías), GetLengthsGeodesic (devuelve las longitudes geodésicas de cada polilínea) y DensifyGeodesic (densifica una polilínea trazando puntos a lo largo de las líneas geodésicas que conectan los vértices, usa IPolycurve4: : GeodesicDensify) métodos.

Como se mencionó en las otras respuestas, ArcGIS todavía usa principalmente cálculos planos.

Melita Kennedy


Algunos comentarios sobre las otras respuestas (¡no hay suficientes representantes para comentar directamente!).

La proyección equidistante azimutal de Esri admite elipsoides. El código GreatCircleDist está creando un PCS que utiliza un GCS basado en elipsoide / esferoide, por lo tanto, las distancias desde el centro / punto de origen serán distancias geodésicas, no grandes distancias circulares. También podría simplificarse. Conocemos las coordenadas proyectadas del primer punto porque es el centro de la proyección: 0,0. Por lo tanto, solo es necesario proyectar el segundo punto. Entonces se podría usar una función EuclidianDistance simplificada.

Verifiqué los resultados con las funciones geodésicas de pe.dll y coincidió. Parece que la aplicación de Richie está usando una esfera, por lo que está devolviendo distancias / coordenadas de gran círculo en su aplicación de prueba. Es por eso que los resultados no coinciden. No reconocí los valores del radio; ¡Creo que necesito hablar con él al respecto!

mkennedy
fuente
2
Melita - ¡Qué bueno verte aquí!
Kirk Kuykendall el
1
Estoy de acuerdo, bienvenido a bordo!
Matt Wilkie
8

La precisión de cualquier respuesta sobre ArcGIS está sujeta a cambios en cualquier momento: por lo que sabemos, se introducirán nuevos procedimientos en el próximo paquete de servicios sin advertencia ni documentación. Dicho esto, el software ESRI ha utilizado durante mucho tiempo cálculos euclidianos ( por ejemplo , la fórmula de Pitágoras para distancias) siempre que se utilizan coordenadas proyectadas. A menudo, en cálculos como el que ilustra, el software ni siquiera tiene acceso a la información de proyección, entonces, ¿qué más se puede hacer?

Su pregunta en sí parece sugerir que los cálculos de distancia euclidiana para una proyección equidistante son correctos. Nada mas lejos de la verdad. Para una proyección equidistante de un punto, se garantiza que la distancia euclidiana al punto base sea ​​igual a la distancia geodésica; Para una proyección equidistante de dos puntos, se garantiza que la distancia euclidiana a cualquier punto base sea igual a la distancia geodésica. A cambio de esas garantías, la distorsión métrica entre todos los otros pares de puntos generalmente aumenta considerablemente en comparación con otras proyecciones que uno podría elegir.

whuber
fuente
@whuber: Gracias por responder. Con respecto al primer párrafo: pensé que ArcGIS podría ver que se usa la proyección del mapa CH1903 (que usa el elipsoide Bessel 1841), y luego proyectar los puntos nuevamente en ese elipsoide a través del dato, y luego hacer cálculos de distancia en el elipsoide. De su respuesta, supongo que ArcGIS no hará todo eso y permanecerá en el espacio Euclidiano XY para hacer cálculos. (¿Qué tal otro software SIG?) - 2º párrafo: Tienes razón, por supuesto, gracias por aclarar este punto.
stakx
Un mecanismo de reproyección oculto es posible solo si los objetos puntuales mantienen referencias a una proyección. No creo que lo hagan.
whuber
@whuber: ¿Sería suficiente (para cálculos más precisos) conocer el elipsoide utilizado para la proyección? AFAIK, ArcGIS almacena una referencia a la proyección utilizada con cada clase de entidad (capa de datos).
stakx
1
En realidad, IPoint, que se deriva de IGeometry, tiene SpatialReference como una propiedad. help.arcgis.com/en/sdk/10.0/arcobjects_net/componenthelp/… Sin embargo, no creo que ReturnDistance lo use. Sin embargo, podría valer la pena probar para ver si eso ha cambiado.
Kirk Kuykendall
1
@stakx He actualizado mi respuesta para incluir código que muestre que la configuración de espacial no tiene ningún impacto en ReturnDistance.
Kirk Kuykendall el