Come calcolare l'area di un poligono sulla superficie terrestre usando python?

Question

Il titolo dice praticamente tutto. Ho bisogno di calcolare l'area all'interno di un poligono sulla superficie terrestre usando Python. Calculating area enclosed by arbitrary polygon on Earth's surface dice qualcosa al riguardo, ma rimane vago sui dettagli tecnici:

Se si vuole fare questo con un più "GIS" sapore, quindi è necessario selezionare un'unità-di-misura per la vostra zona e trovare una proiezione appropriata che conserva l'area (non tutti lo fanno). Dal momento che l'utente sta parlando di calcolare un poligono arbitrario , utilizzerei lo qualcosa di simile a una proiezione di area uguale di Lambert Azimuthal . Impostare l'origine /centro della proiezione di essere il centro del poligono, progetto il poligono per il nuovo sistema di coordinate , quindi calcolare l'area utilizzando tecniche planari standard.

Quindi, come faccio a farlo in Python?

Answer 1

Diciamo che si dispone di una rappresentazione dello stato del Colorado in formato GeoJSON

{"type": "Polygon", 
"coordinates": [[ 
    [-102.05, 41.0], 
    [-102.05, 37.0], 
    [-109.05, 37.0], 
    [-109.05, 41.0] 
]]} 

Tutte le coordinate sono la longitudine, latitudine. È possibile utilizzare pyproj per proiettare le coordinate e Shapely per trovare l'area di un poligono proiettata:

co = {"type": "Polygon", "coordinates": [ 
    [(-102.05, 41.0), 
    (-102.05, 37.0), 
    (-109.05, 37.0), 
    (-109.05, 41.0)]]} 
lon, lat = zip(*co['coordinates'][0]) 
from pyproj import Proj 
pa = Proj("+proj=aea +lat_1=37.0 +lat_2=41.0 +lat_0=39.0 +lon_0=-106.55") 

Questa è una proiezione superficie pari centrata sulla e bracketing l'area di interesse. Ora fate nuova rappresentazione GeoJSON proiettata, si trasformano in un oggetto geometrico ben fatto, e prendere l'area:

x, y = pa(lon, lat) 
cop = {"type": "Polygon", "coordinates": [zip(x, y)]} 
from shapely.geometry import shape 
shape(cop).area # 268952044107.43506 

E 'una buona approssimazione della zona esaminata. Per caratteristiche più complesse, è necessario campionare lungo i bordi, tra i vertici, per ottenere valori accurati. Si applicano tutte le avvertenze sopra relative alle linee dati, ecc. Se sei interessato solo all'area, puoi tradurre la tua funzione dalla linea di demarcazione prima di proiettarla.

Answer 2

Il modo più semplice per farlo (secondo me) è di proiettare le cose in una (molto semplice) proiezione di area uguale e utilizzare una delle consuete tecniche planari per calcolare l'area.

Prima di tutto, assumerò che una terra sferica è abbastanza vicina ai tuoi scopi, se stai facendo questa domanda. In caso contrario, è necessario riproiettare i dati utilizzando un ellissoide appropriato, nel qual caso si vorrà utilizzare una vera e propria libreria di proiezione (tutto utilizza proj4 dietro le quinte, in questi giorni) come i collegamenti Python a GDAL/OGR o (il molto più amichevole) pyproj.

Tuttavia, se stai bene con una terra sferica, è abbastanza semplice farlo senza librerie specializzate.

La proiezione di area uguale più semplice da calcolare è una sinusoidal projection. Fondamentalmente, si moltiplica la latitudine per la lunghezza di un grado di latitudine e la longitudine per la lunghezza di un grado di latitudine e il coseno della latitudine.

def reproject(latitude, longitude): 
    """Returns the x & y coordinates in meters using a sinusoidal projection""" 
    from math import pi, cos, radians 
    earth_radius = 6371009 # in meters 
    lat_dist = pi * earth_radius/180.0 

    y = [lat * lat_dist for lat in latitude] 
    x = [long * lat_dist * cos(radians(lat)) 
       for lat, long in zip(latitude, longitude)] 
    return x, y 

Ok ... Ora tutto quello che dobbiamo fare è quello di calcolare l'area di un poligono arbitrario in un piano.

Ci sono diversi modi per farlo. Userò ciò che è probably the most common one qui.

def area_of_polygon(x, y): 
    """Calculates the area of an arbitrary polygon given its verticies""" 
    area = 0.0 
    for i in range(-1, len(x)-1): 
     area += x[i] * (y[i+1] - y[i-1]) 
    return abs(area)/2.0 

Speriamo che si punta nella direzione giusta, in ogni caso ...

Answer 3

Perché la terra è una superficie chiusa un poligono chiuso disegnato sulla sua superficie crea DUE aree poligonali. È inoltre necessario definire quale è all'interno e quale all'esterno!

La maggior parte delle volte le persone avranno a che fare con piccoli poligoni, quindi è "ovvio", ma una volta che si hanno le dimensioni di oceani o continenti, è meglio assicurarsi di ottenere questo nel modo giusto.

Inoltre, ricorda che le linee possono passare da (-179,0) a (+179,0) in due modi diversi. Uno è molto più lungo dell'altro. Ancora una volta, per lo più si suppone che questa sia una linea che va da (-179,0) a (-180,0) che è (+180,0) e quindi a (+179,0), ma uno giorno ... non lo farà.

Trattare lat-long come un semplice sistema di coordinate (x, y), o anche trascurare il fatto che qualsiasi proiezione di coordinate avrà distorsioni e interruzioni, può farti fallire a grandi passi nelle sfere.

Answer 4

Un po 'tardi, forse, ma qui c'è un metodo diverso, usando il teorema di Girard. Afferma che l'area di un poligono di grandi cerchi è R ** 2 volte la somma degli angoli tra i poligoni meno (N-2) * pi dove N è il numero di angoli.

Ho pensato che sarebbe valso la pena di postare, dal momento che non si basa su nessun'altra libreria che su numpy, ed è un metodo piuttosto diverso rispetto agli altri. Naturalmente, questo funziona solo su una sfera, quindi ci sarà qualche imprecisione quando si applica alla Terra.

Innanzitutto, definire una funzione per calcolare l'angolo di rilevamento dal punto 1 lungo un cerchio al punto 2:

import numpy as np 
from numpy import cos, sin, arctan2 

d2r = np.pi/180 

def greatCircleBearing(lon1, lat1, lon2, lat2): 
    dLong = lon1 - lon2 

    s = cos(d2r*lat2)*sin(d2r*dLong) 
    c = cos(d2r*lat1)*sin(d2r*lat2) - sin(lat1*d2r)*cos(d2r*lat2)*cos(d2r*dLong) 

    return np.arctan2(s, c) 

Ora posso usare questo per trovare gli angoli, e quindi l'area (Nella seguente, lons ei dorsali dovrebbero naturalmente essere specificati, e dovrebbero essere nel giusto ordine. Inoltre, il raggio della sfera deve essere specificata.)

N = len(lons) 

angles = np.empty(N) 
for i in range(N): 

    phiB1, phiA, phiB2 = np.roll(lats, i)[:3] 
    LB1, LA, LB2 = np.roll(lons, i)[:3] 

    # calculate angle with north (eastward) 
    beta1 = greatCircleBearing(LA, phiA, LB1, phiB1) 
    beta2 = greatCircleBearing(LA, phiA, LB2, phiB2) 

    # calculate angle between the polygons and add to angle array 
    angles[i] = np.arccos(cos(-beta1)*cos(-beta2) + sin(-beta1)*sin(-beta2)) 

area = (sum(angles) - (N-2)*np.pi)*R**2 

con il Colorado coordinate fornite in un'altra risposta, e con Raggio di terra 6371 km, ho capito che l'area è 268930758560.74808

Answer 5

o semplicemente utilizzare una libreria: https://github.com/scisco/area

from area import area 
>>> obj = {'type':'Polygon','coordinates':[[[-180,-90],[-180,90],[180,90],[180,-90],[-180,-90]]]} 
>>> area(obj) 
511207893395811.06 

... restituisce l'area in metri quadrati.

Answer 6

Ecco una soluzione che utilizza basemap, anziché pyproj e shapely, per la conversione delle coordinate. L'idea è la stessa suggerita da @sgillies però. NOTA che ho aggiunto il 5 ° punto in modo che il percorso sia un ciclo chiuso.

import numpy 
from mpl_toolkits.basemap import Basemap 

coordinates=numpy.array([ 
[-102.05, 41.0], 
[-102.05, 37.0], 
[-109.05, 37.0], 
[-109.05, 41.0], 
[-102.05, 41.0]]) 

lats=coordinates[:,1] 
lons=coordinates[:,0] 

lat1=numpy.min(lats) 
lat2=numpy.max(lats) 
lon1=numpy.min(lons) 
lon2=numpy.max(lons) 

bmap=Basemap(projection='cea',llcrnrlat=lat1,llcrnrlon=lon1,urcrnrlat=lat2,urcrnrlon=lon2) 
xs,ys=bmap(lons,lats) 

area=numpy.abs(0.5*numpy.sum(ys[:-1]*numpy.diff(xs)-xs[:-1]*numpy.diff(ys))) 
area=area/1e6 

print area 

Il risultato è 268993.609651 in km^2.