programa - ¿Cómo calcular el área de un polígono en la superficie de la tierra usando python?

formula para calcular areas irregulares

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. El área de cálculo delimitada por un polígono arbitrario en la superficie de la Tierra dice algo al respecto, pero sigue siendo vaga en los detalles técnicos:

Si desea hacer esto con un sabor más "GIS", entonces debe seleccionar una unidad de medida para su área y encontrar una proyección adecuada que conserve el área (no todas lo hacen). Ya que está hablando de calcular un polígono arbitrario, usaría algo como una proyección de Lambert Azimuthal Equal Area. 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.

Entonces, ¿cómo hago esto en Python?

Answer 1

Aquí hay una solución que usa el basemap , en lugar de pyproj y shapely , para la conversión de coordenadas. Sin embargo, la idea es la misma que sugiere @sgillies. TENGA EN CUENTA que he agregado 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.

Answer 2

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

La mayoría de las veces la gente tratará con polígonos pequeños, por lo que es ''obvio'', pero una vez que tenga cosas del tamaño de océanos o continentes, será mejor que se asegure de hacerlo de la manera correcta.

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

Tratar el lat-long como un sistema de coordenadas simple (x, y), o incluso ignorar el hecho de que cualquier proyección de coordenadas tendrá distorsiones y roturas, puede hacer que falle en grandes esferas.

Answer 3

Digamos que 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 longitud, latitud. Puedes usar 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 una proyección de área igual centrada en el área de interés. Ahora haga una nueva representación GeoJSON proyectada, conviértala en un objeto geométrico bien formado y tome el área:

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 cercana al área encuestada. Para características más complejas, necesitará muestrear a lo largo de los bordes, entre los vértices, para obtener valores precisos. Todas las advertencias anteriores sobre líneas de datos, etc., se aplican. Si solo le interesa el área, puede traducir su función fuera de la línea de datos antes de proyectar.

Answer 4

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

En primer lugar, voy a suponer que una tierra esférica está lo suficientemente cerca para sus propósitos, si está haciendo esta pregunta. Si no, entonces necesita reproyectar sus datos usando un elipsoide apropiado, en cuyo caso usted va a querer usar una biblioteca de proyección real (todo lo que se usa proj4 detrás de escena, en estos días) como los enlaces de Python a GDAL/OGR o (el mucho más amigable) pyproj .

Sin embargo, si estás de acuerdo con una tierra esférica, es bastante simple hacerlo sin ninguna biblioteca especializada.

La proyección más simple de igual área para calcular es una proyección sinusoidal . Básicamente, simplemente multiplica 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 muchas maneras de hacer esto. Voy a usar el que probablemente sea el más común 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

Esperemos que eso te dirija en la dirección correcta, de todos modos ...

Answer 5

O simplemente use 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

... vuelve el área en metros cuadrados.

Answer 6

Quizás un poco tarde, pero aquí hay un método diferente, usando el teorema de Girard. Indica 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 valdría la pena publicar esto, ya que no se basa en ninguna otra biblioteca que no sea numpy, y es un método bastante diferente a los demás. Por supuesto, esto solo funciona en una esfera, por lo que habrá algo de inexactitud al aplicarlo a la Tierra.

Primero, defino una función para calcular el ángulo del rodamiento desde el punto 1 a lo largo de un gran círculo hasta el 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 lo que sigue, por supuesto, se deben especificar lones y tiempos, y deben estar en el orden correcto. Además, se debe especificar el radio de la esfera).

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