¿Cómo se calcula el área de un polígono en la superficie de la tierra con python?

Question

El título básicamente lo dice todo. Necesito calcular el área dentro de un polígono en la superficie de la Tierra usando Python. Calculating area enclosed by arbitrary polygon on Earth's surface dice algo al respecto, pero sigue siendo vago en los detalles técnicos:

Si desea hacer esto con un sabor más "SIG", entonces usted necesita para seleccionar una unidad de medida para su área y encuentran una proyección adecuada que conserva el área (no todas lo hacen). Dado que están hablando de calcular un polígono arbitrario , usaría algo así como una proyección de área igual de Lambert Azimuthal . 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 , luego calcule el área usando las técnicas planas estándar .

Entonces, ¿cómo hago esto en Python?

Answer 1

Digamos que usted tiene una representación del estado de Colorado en formato GeoJSON

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

Todas las coordenadas son la longitud, latitud. Puede utilizar pyproj para proyectar las coordenadas y Shapely para encontrar el área de cualquier polígono proyectado:

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") 

Esa es un área igual de proyección centrada en y entre paréntesis el área de interés. Ahora hacer una nueva representación GeoJSON proyectada, se convierten en un objeto geométrico bien proporcionado, y tomar la zona:

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

Es una aproximación muy cerca de la zona estudiada. Para características más complejas, tendrá que muestrear a lo largo de los bordes, entre los vértices, para obtener valores precisos. Se aplican todas las advertencias anteriores sobre datelines, etc. Si solo está interesado en el área, puede traducir su función fuera de la fecha antes de proyectar.

Answer 2

La manera más fácil de hacer esto (en mi opinión), es proyectar cosas en (una muy simple) proyección de área igual y usar una de las técnicas planas usuales para calcular el área.

En primer lugar, voy a suponer que una tierra esférica es lo suficientemente cerca para sus propósitos, si está haciendo esta pregunta. De lo contrario, deberá reproyectar sus datos utilizando un elipsoide apropiado, en cuyo caso querrá usar una biblioteca de proyección real (todo lo que usa proj4 detrás de escena, actualmente) como las vinculaciones de Python al GDAL/OGR o (el mucho más amigable) pyproj.

Sin embargo, si está de acuerdo con una tierra esférica, es bastante simple hacerlo sin bibliotecas especializadas.

La proyección de área igual más simple para calcular es sinusoidal projection. Básicamente, simplemente multiplique la latitud por la longitud de un grado de latitud, y la longitud por la longitud de un grado de latitud y el coseno de la latitud.

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 

Bien ... Ahora todo lo que tenemos que hacer es calcular el área de un polígono arbitrario en un plano.

Hay varias formas de hacerlo. Voy a usar lo que es probably the most common one aquí.

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 

suerte que le apuntan en la dirección correcta, de todos modos ...

Answer 3

Debido a que la tierra es una superficie cerrada un polígono cerrado dibujado en su superficie crea DOS áreas poligonales. ¡También necesita definir cuál está adentro y cuál afuera!

La mayoría de las veces las personas lidiarán con polígonos pequeños, por lo que es 'obvio', pero una vez que tenga cosas del tamaño de océanos o continentes, es mejor que se asegure de hacerlo bien.

Además, recuerde que las líneas pueden ir desde (-179,0) hasta (+179,0) de dos maneras diferentes. Uno es mucho más largo que el otro. Una vez más, en su mayoría supondrá que se trata de una línea que va de (-179,0) a (-180,0) que es (+180,0) y luego a (+179,0), pero uno día ... no lo hará.

Tratar lat-long como un sistema de coordenadas simple (x, y), o incluso descuidar el hecho de que cualquier proyección de coordenadas va a tener distorsiones y roturas, puede hacer que falles a gran velocidad en las esferas.

Answer 4

Un poco tarde quizás, pero aquí hay un método diferente, usando el teorema de Girard. Establece que el área de un polígono de círculos grandes es R ** 2 veces la suma de los ángulos entre los polígonos menos (N-2) * pi, donde N es el número de esquinas.

Pensé que esto valdría la pena publicarlo, ya que no depende de ninguna otra biblioteca que numpy, y es un método bastante diferente al resto. Por supuesto, esto solo funciona en una esfera, por lo que habrá cierta inexactitud al aplicarlo a la Tierra.

Primero, definir una función para calcular el ángulo de marcación desde el punto 1 a lo largo de un gran círculo 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) 

Ahora puedo usar esto para encontrar los ángulos, y luego el área (en el siguiente, Lons y dorsales deben por supuesto ser especificados, y que debe estar en el orden correcto. Además, el radio de la esfera debe ser especificado.)

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 el Colorado coordenadas dadas en otra respuesta, y con Radio de la Tierra 6371 km, entiendo que el área es 268930758560.74808

Answer 5

o simplemente utilizar una biblioteca: 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 

... devuelve el área en metros cuadrados.

Answer 6

Aquí hay una solución que usa basemap, en lugar de pyproj y shapely, para la conversión de coordenadas. La idea es la misma sugerida por @sgillies. TEN EN CUENTA que agregué el quinto punto para que la ruta sea un ciclo cerrado.

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 

El resultado es 268993.609651 en km^2.