Comment calculer la superficie d'un polygone sur la surface de la terre en utilisant python?

la façon la plus facile de faire cela (à mon avis), est de projeter les choses dans (une projection de zone) égale (très simple) et d'utiliser l'une des techniques planaires habituelles pour calculer la zone.

tout d'abord, je vais supposer qu'une Terre sphérique est assez proche pour vos fins, Si vous posez cette question. Dans le cas contraire, vous devez alors reprogrammer vos données en utilisant un ellipsoïde approprié, auquel cas vous allez vouloir utiliser une bibliothèque de projection réelle (tout utilise proj4 dans les coulisses, ces jours-ci) tels que les fixations de python à GDAL/OGR ou (le beaucoup plus amical) pyproj .

cependant, si vous êtes d'accord avec une Terre sphérique, il est assez simple de le faire sans bibliothèques spécialisées.

la projection de surface égale la plus simple à calculer est une projection sinusoïdale . En gros, on multiplie la latitude par la longueur de un degré de latitude et la longitude de la longueur d'un degré de latitude et le cosinus de la latitude.

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... Maintenant, tout ce que nous avons à faire est de calculer la surface d'un polygone arbitraire dans un plan.

il y a plusieurs façons de le faire. Je vais utiliser ce qui est probablement le plus commun ici.

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

espérons que cela vous mènera dans la bonne direction, de toute façon...

un peu tard peut-être, Mais voici une méthode différente, en utilisant le théorème de Girard. Il indique que la surface d'un polygone de grands cercles est R**2 fois la somme des angles entre les polygones moins (n-2)*pi où N est le nombre de coins.

j'ai pensé que cela vaudrait la peine de poster, car il ne s'appuie sur aucune autre bibliothèque que numpy, et c'est une méthode très différente que les autres. Bien sûr, cela ne fonctionne que sur une sphère, donc il y aura une certaine inexactitude quand en l'appliquant à la Terre.

tout d'abord, je définit une fonction pour calculer l'angle de roulement à partir du point 1 le long d'un grand cercle jusqu'au point 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)

maintenant je peux utiliser ceci pour trouver les angles, et puis la zone (dans ce qui suit, lons et lats devraient bien sûr être spécifiés, et ils devraient être dans le bon ordre. Aussi, le rayon de la sphère doit être spécifié.)

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

avec les coordonnées du Colorado indiquées dans une autre réponse, et avec le rayon de la Terre 6371 km, je comprends que la zone est 268930758560.74808

parce que la Terre est une surface fermée un polygone fermé dessiné sur sa surface crée deux zones polygonales. Vous devez également définir qui l'on est à l'intérieur et à l'extérieur!

la plupart du temps, les gens ont affaire à de petits polygones, et donc c'est "évident", mais une fois que vous avez des choses de la taille d'un océan ou d'un continent, vous feriez mieux de vous assurer que vous obtenez cela de la bonne façon.

aussi, rappelez-vous que les lignes peuvent aller de (-179,0) à (+179,0) de deux façons différentes. L'un est beaucoup plus longue que l'autre. Encore une fois, la plupart du temps vous supposerez qu'il s'agit d'une ligne qui va de (-179,0) à (-180,0) qui est (+180,0) et puis à (+179,0), mais un jour... ça n'arrivera pas.

traiter lat-long comme un simple (x,y) système de coordonnées, ou même négliger le fait que toute projection de coordonnées va avoir des distorsions et des ruptures, peut vous faire échouer à grand-temps sur les sphères.

Voici une solution qui utilise basemap , au lieu de pyproj et shapely , pour la conversion de coordonnées. L'idée est la même que celle suggérée par @sgillies. Notez que j'ai ajouté le 5ème point pour que le chemin soit une boucle fermée.

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

le résultat est 268993.609651 en km^2.

ou tout simplement utiliser une bibliothèque: 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

...renvoie la surface en mètres carrés.