math - pick - Cálculo del área encerrada por un polígono arbitrario en la superficie de la Tierra

shape length arcgis

Question

Digamos que tengo un conjunto arbitrario de pares de latitud y longitud que representan puntos en alguna curva cerrada simple. En el espacio cartesiano, podría calcular fácilmente el área delimitada por dicha curva usando el Teorema de Green. ¿Cuál es el enfoque análogo para calcular el área en la superficie de una esfera? Creo que lo que busco es (incluso una aproximación) del algoritmo detrás de la función areaint de Matlab .

Answer 1

No sé nada sobre la función de Matlab, pero aquí vamos. Considere dividir su polígono esférico en triángulos esféricos, por ejemplo dibujando diagonales desde un vértice. El área de superficie de un triángulo esférico está dada por

R^2 * ( A + B + C - /pi)

donde R es el radio de la esfera, y A , B y C son los ángulos interiores del triángulo (en radianes). La cantidad entre paréntesis se conoce como el "exceso esférico".

Su polígono de n lados se dividirá en n-2 triángulos. Sumando todos los triángulos, extrayendo el factor común de R^2 y juntando todos los /pi , el área de su polígono es

R^2 * ( S - (n-2)/pi )

donde S es la suma de ángulo de su polígono. La cantidad entre paréntesis es nuevamente el exceso esférico del polígono.

[edit] Esto es verdadero independientemente de si el polígono es convexo o no. Todo lo que importa es que se puede diseccionar en triángulos.

Puedes determinar los ángulos a partir de un poco de matemática vectorial. Supongamos que tiene tres vértices A , B , C y está interesado en el ángulo en B Por lo tanto, debemos encontrar dos vectores tangentes (sus magnitudes son irrelevantes) a la esfera desde el punto B largo de los segmentos del gran círculo (los bordes del polígono). Vamos a resolverlo para BA . El gran círculo se encuentra en el plano definido por OA y OB , donde O es el centro de la esfera, por lo que debe ser perpendicular al vector normal OA x OB . También debería ser perpendicular al OB ya que es tangente allí. Tal vector es, por lo tanto, dado por OB x (OA x OB) . Puede usar la regla de la mano derecha para verificar que esté en la dirección adecuada. Tenga en cuenta también que esto se simplifica a OA * (OB.OB) - OB * (OB.OA) = OA * |OB| - OB * (OB.OA) OA * (OB.OB) - OB * (OB.OA) = OA * |OB| - OB * (OB.OA) .

A continuación, puede utilizar el producto de buen viejo punto para buscar el ángulo entre los lados: BA''.BC'' = |BA''|*|BC''|*cos(B) , donde BA'' y BC'' son los vectores tangentes de B a los lados de A y C

[editado para que quede claro que estos son vectores tangentes, no literales entre los puntos]

Answer 2

Hay varias formas de hacer esto.

1) Integrar las contribuciones de las bandas latitudinales. Aquí el área de cada tira será (Rcos (A) (B1-B0)) (RdA), donde A es la latitud, B1 y B0 son las longitudes inicial y final, y todos los ángulos están en radianes.

2) Rompa la superficie en triángulos esféricos , y calcule el área usando el Teorema de Girard , y sume estos.

3) Como sugiere James Schek, en el trabajo GIS utilizan un área que preserva la proyección en un espacio plano y calculan el área allí.

De la descripción de sus datos, en sonidos como el primer método podría ser el más fácil. (Por supuesto, puede haber otros métodos más fáciles que no conozco).

Hay un documento sobre esto , si tiene acceso, pero me topé con él buscando "bandas longitudinales", así que supongo que están usando el método 1 anterior.

Editar - comparando estos dos métodos:

En la primera inspección, puede parecer que el enfoque del triángulo esférico es más fácil, pero, en general, este no es el caso. El problema es que uno no solo necesita dividir la región en triángulos, sino en triángulos esféricos , es decir, triángulos cuyos lados son grandes arcos de círculo. Por ejemplo, los límites latitudinales no califican , por lo que estos límites deben dividirse en aristas que se aproximen mejor a los arcos del gran círculo. Y esto se vuelve más difícil de hacer para bordes arbitrarios donde los círculos grandes requieren combinaciones específicas de ángulos esféricos. Considere, por ejemplo, cómo se rompería una banda intermedia alrededor de una esfera, por ejemplo, toda el área entre lat 0 y 45deg en triángulos esféricos.

Al final, si uno debe hacer esto correctamente con errores similares para cada método, el método 2 dará menos triángulos, pero serán más difíciles de determinar. El método 1 da más tiras, pero son triviales para determinar. Por lo tanto, sugiero el método 1 como el mejor enfoque.

Answer 3

Usted menciona "geografía" en una de sus etiquetas, por lo que solo puedo suponer que está detrás del área de un polígono en la superficie de un geoide. Normalmente, esto se hace usando un sistema de coordenadas proyectadas en lugar de un sistema de coordenadas geográficas (es decir, lon / lat). Si tuviera que hacerlo en lon / lat, entonces asumiría que la unidad de medida devuelta sería el porcentaje de la superficie de la esfera.

Si desea hacer esto con un sabor más "GIS", entonces necesita seleccionar una unidad de medida para su área y encontrar una proyección apropiada que preserve el área (no todas). Ya que está hablando de calcular un polígono arbitrario, usaría algo así como una proyección de Área Igualitaria Azimutal de Lambert . Establezca el origen / centro de la proyección para que sea el centro de su polígono, proyecte el polígono al nuevo sistema de coordenadas y luego calcule el área utilizando técnicas planas estándar.

Si necesita hacer muchos polígonos en un área geográfica, es probable que haya otras proyecciones que funcionen (o que estén lo suficientemente cerca). UTM, por ejemplo, es una excelente aproximación si todos sus polígonos están agrupados alrededor de un solo meridiano.

No estoy seguro de si algo de esto tiene algo que ver con el funcionamiento de la función areaint de Matlab.

Answer 4

Reescribí la función "areaint" de MATLAB en java, que tiene exactamente el mismo resultado. "areaint" calcula la "superficie por unidad", por lo que multipliqué la respuesta por el área de la superficie de la Tierra (5.10072e14 metros cuadrados).

private double area(ArrayList<Double> lats,ArrayList<Double> lons)
{       
    double sum=0;
    double prevcolat=0;
    double prevaz=0;
    double colat0=0;
    double az0=0;
    for (int i=0;i<lats.size();i++)
    {
        double colat=2*Math.atan2(Math.sqrt(Math.pow(Math.sin(lats.get(i)*Math.PI/180/2), 2)+ Math.cos(lats.get(i)*Math.PI/180)*Math.pow(Math.sin(lons.get(i)*Math.PI/180/2), 2)),Math.sqrt(1-  Math.pow(Math.sin(lats.get(i)*Math.PI/180/2), 2)- Math.cos(lats.get(i)*Math.PI/180)*Math.pow(Math.sin(lons.get(i)*Math.PI/180/2), 2)));
        double az=0;
        if (lats.get(i)>=90)
        {
            az=0;
        }
        else if (lats.get(i)<=-90)
        {
            az=Math.PI;
        }
        else
        {
            az=Math.atan2(Math.cos(lats.get(i)*Math.PI/180) * Math.sin(lons.get(i)*Math.PI/180),Math.sin(lats.get(i)*Math.PI/180))% (2*Math.PI);
        }
        if(i==0)
        {
             colat0=colat;
             az0=az;
        }           
        if(i>0 && i<lats.size())
        {
            sum=sum+(1-Math.cos(prevcolat  + (colat-prevcolat)/2))*Math.PI*((Math.abs(az-prevaz)/Math.PI)-2*Math.ceil(((Math.abs(az-prevaz)/Math.PI)-1)/2))* Math.signum(az-prevaz);
        }
        prevcolat=colat;
        prevaz=az;
    }
    sum=sum+(1-Math.cos(prevcolat  + (colat0-prevcolat)/2))*(az0-prevaz);
    return 5.10072E14* Math.min(Math.abs(sum)/4/Math.PI,1-Math.abs(sum)/4/Math.PI);
}