Obliczanie powierzchni ogrodzonej dowolnym wielokątem na powierzchni Ziemi

Przepisałem ponownie funkcję "areaint" Matlaba w Javie, która ma dokładnie taki sam wynik. "areaint" oblicza "suface na jednostkę", więc pomnożyłem odpowiedź przez powierzchnię Ziemi (5. 10072e14 m2).

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);
}

W jednym ze swoich tagów wspominasz o "geografii", więc mogę tylko założyć, że chodzi Ci o obszar wielokąta na powierzchni geoidy. Zwykle odbywa się to za pomocą rzutowanego układu współrzędnych, a nie geograficznego układu współrzędnych (tj. lon / lat). Jeśli miałbyś to zrobić w lon / lat, to założyłbym, że zwrócona jednostka miary byłaby procentami powierzchni kuli.

Jeśli chcesz to zrobić z bardziej" GIS " smak, następnie trzeba wybrać jednostkę miary dla swojego obszaru i znaleźć odpowiedni rzut, który zachowuje obszar(nie wszystkie). Skoro mówisz o obliczaniu dowolnego wielokąta, użyłbym czegoś w rodzaju Azymutalnej równej powierzchni projekcji Lamberta. Ustaw początek / środek projekcji jako środek wielokąta, przekaż wielokąt do nowego układu współrzędnych, a następnie Oblicz obszar za pomocą standardowych technik planarnych.

Jeśli trzeba zrobić wiele wielokątów w obszarze geograficznym, są prawdopodobnie inne projekcje, które będą pracy (lub będzie wystarczająco blisko). UTM, na przykład, jest doskonałym przybliżeniem, jeśli wszystkie Twoje wielokąty są skupione wokół pojedynczego południka.

Nie jestem pewien, czy cokolwiek z tego ma coś wspólnego z tym, jak działa funkcja areaint Matlaba.

Nie wiem nic o funkcji Matlaba, ale zaczynamy. Rozważ podzielenie wielokąta sferycznego na Trójkąty sferyczne, np. rysując przekątne z wierzchołka. Pole powierzchni trójkąta sferycznego jest podane przez

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

Gdzie R jest promieniem kuli, oraz A, B, i C są wewnętrznymi kątami trójkąta (w radianach). Ilość w nawiasach jest znana jako "nadmiar kulisty".

Twój n-wielokąt jednostronny zostanie podzielony w n-2 Trójkąty. Sumując wszystkie trójkąty, wydobywając wspólny czynnik R^2 i łącząc wszystkie \pi razem, obszar twojego wielokąta wynosi

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

Gdzie S jest sumą kątów twojego wielokąta. Ilość w nawiasach jest ponownie sferycznym nadmiarem wielokąta.

[edytuj] to prawda, czy wielokąt jest wypukły. Liczy się tylko to, że można ją podzielić na Trójkąty.

Możesz określić kąty z bitu matematyki wektorowej. Załóżmy, że masz trzy wierzchołki A,B,C i są zainteresowani kątem w B. Musimy zatem znaleźć dwa wektory styczne (ich wielkość nie ma znaczenia) do sfery z punktu B wzdłuż segmentów Wielkiego okręgu (krawędzi wielokąta). Rozwiążmy to dla BA. Wielki okrąg leży w płaszczyźnie określonej przez OA i OB, gdzie O jest środkiem kuli, więc powinien być prostopadły do wektora normalnego OA x OB. Powinna być również prostopadła do OB, ponieważ jest tam styczna. Taki wektor jest zatem dany przez OB x (OA x OB). Możesz użyć reguły prawej ręki, aby sprawdzić, czy jest to w odpowiednim kierunku. Zauważ również, że upraszcza to OA * (OB.OB) - OB * (OB.OA) = OA * |OB| - OB * (OB.OA).

Możesz następnie użyć starego dobrego produktu do znalezienia kąta między bokami: BA'.BC' = |BA'|*|BC'|*cos(B), gdzie BA' i BC' są wektorami stycznymi od B wzdłuż boków do A i C.

[edytowane, aby było jasne, że są to wektory styczne, a nie literalne między punktami]

Oto implementacja Pythona 3, luźno inspirowana powyższymi odpowiedziami:

def polygon_area(lats, lons, algorithm = 0, radius = 6378137):
    """
    Computes area of spherical polygon, assuming spherical Earth. 
    Returns result in ratio of the sphere's area if the radius is specified.
    Otherwise, in the units of provided radius.
    lats and lons are in degrees.
    """
    from numpy import arctan2, cos, sin, sqrt, pi, power, append, diff, deg2rad
    lats = np.deg2rad(lats)
    lons = np.deg2rad(lons)

    # Line integral based on Green's Theorem, assumes spherical Earth

    #close polygon
    if lats[0]!=lats[-1]:
        lats = append(lats, lats[0])
        lons = append(lons, lons[0])

    #colatitudes relative to (0,0)
    a = sin(lats/2)**2 + cos(lats)* sin(lons/2)**2
    colat = 2*arctan2( sqrt(a), sqrt(1-a) )

    #azimuths relative to (0,0)
    az = arctan2(cos(lats) * sin(lons), sin(lats)) % (2*pi)

    # Calculate diffs
    # daz = diff(az) % (2*pi)
    daz = diff(az)
    daz = (daz + pi) % (2 * pi) - pi

    deltas=diff(colat)/2
    colat=colat[0:-1]+deltas

    # Perform integral
    integrands = (1-cos(colat)) * daz

    # Integrate 
    area = abs(sum(integrands))/(4*pi)

    area = min(area,1-area)
    if radius is not None: #return in units of radius
        return area * 4*pi*radius**2
    else: #return in ratio of sphere total area
        return area

Proszę znaleźć nieco bardziej wyraźną wersję (i wiele więcej odniesień i TODOs...) tutaj .

Możesz również spojrzeć na ten kod spherical_geometry pakietu: Tutaj i tutaj . Podaje dwie różne metody obliczania powierzchni wielokąta sferycznego.