Esta es la solución más completa que encomtrarás dar, pero mírala pausadamente y analiza si se puede adaptar a tu trabajo.
Solución:
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. Puede 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 áreas iguales centrada y entre corchetes el área de interés. Ahora haga una nueva representación GeoJSON proyectada, conviértala en un objeto geométrico Shapely 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. Se aplican todas las advertencias anteriores sobre fechas, etc. Si solo está interesado en el área, puede traducir su característica fuera de la línea de fecha antes de proyectar.
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. De lo contrario, debe reproyectar sus datos usando un elipsoide apropiado, en cuyo caso querrá usar una biblioteca de proyección real (todo 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á de acuerdo con una tierra esférica, es bastante simple hacerlo sin bibliotecas especializadas.
La proyección de áreas iguales más simple de 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
Con suerte, eso le indicará la dirección correcta, de todos modos …
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
… devuelve el área en metros cuadrados.
Calificaciones y reseñas
Recuerda algo, que tienes la opción de glosar si diste con el hallazgo.