Comment calculer l'aire d'un polygone sur la surface de la terre à l'aide de python?

Le titre, fondamentalement, tout est dit. J'ai besoin de calculer la superficie à l'intérieur d'un polygone sur la surface de la Terre à l'aide de Python. Calcul de la surface délimitée par l'arbitraire d'un polygone sur la surface de la Terre dit quelque chose sur elle, mais reste vague sur les détails techniques:

Si vous voulez le faire avec un plus
"SIG" la saveur, vous devez alors sélectionner
une unité de mesure de votre région et de
trouver la bonne projection
préserve de la zone (mais pas tous). Puisque vous
parlons de calcul d'un
arbitraire polygone, je voudrais utiliser
quelque chose comme un Lambert Azimuthal
L'égalité de la Zone de projection. Définir le
origine/centre de la projection de l'être
le centre de votre polygone, projet
le polygone de nouvelles coordonnées
système, puis de calculer la surface à l'aide de
standard plane techniques.

Alors, comment dois-je le faire en Python?

  • Si vous choisissez une zone de préservation de la projection, les côtés du polygone de ne plus être tout droit
InformationsquelleAutor andreas-h | 2011-01-13
  1. 34

    Disons que vous avez une représentation de l'état du Colorado en format GeoJSON

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

    Toutes les coordonnées de longitude et de latitude. Vous pouvez utiliser pyproj de projet les coordonnées et Galbées pour trouver la surface de tout projet de polygone:

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

    C'est une égalité de la zone de projection centrée sur et bracketing de la zone d'intérêt. Maintenant, faire de nouvelles projetées représentation GeoJSON, se transformer en un Galbe objet géométrique, et de prendre la région:

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

    C'est une très bonne approximation de la zone étudiée. Pour des fonctions plus complexes, vous aurez besoin de l'échantillon sur les bords, entre les sommets, pour obtenir des valeurs précises. Toutes les mises en garde ci-dessus sur les délais, etc, s'appliquent. Si vous êtes seulement intéressés par la région, vous pouvez traduire votre fonction à l'écart de la ligne de changement de date avant la projection.

    • à proprement parler, qui GeoJSON devrait avoir un cinquième point de fermeture, [-102.05, 41.0]. La spec est un peu vague sur ces choses, parfois.
    • quelle est la part de la zone de résultat? est-il en mètres carrés ou en pouces carrés?
    • Dans le cas où quelqu'un d'autre se pose la question, c'est en mètres carrés. Vous pouvez google région du colorado et obtenir 104,185 milles carrés. De conversion que pour mètres carrés donne à peu près le même résultat.
    • Cela fonctionne tout aussi bien pour les polygones qui pourrait ne pas être simple recherche quadrilatères?
    • Donc, en gros, je peux utiliser toute l'égalité de la zone de projection pour convertir les coordonnées et l'utilisation de Green théorème pour calculer l'aire, ai-je le droit? Si oui, alors on peut alternativement utiliser basemap de faire de la projection (par exemple 'aea', 'ace'). En termes de Vert du théorème, un one-liner area=np.abs(0.5*np.sum(y[:-1]*np.diff(x) - x[:-1]*np.diff(y))) vous libère de la nécessité de la shapely module
    InformationsquelleAutor sgillies
  2. 24

    La façon la plus simple de le faire (à mon avis), c'est projeter les choses en (très simple) de l'égalité de la zone de projection et d'utiliser l'un de l'habitude plane techniques pour le calcul de la superficie.

    Tout d'abord, je vais supposer que d'une terre sphérique est assez proche de vos besoins, si vous vous posez cette question. Si non, alors vous avez besoin pour reprojeter vos données à l'aide d'un ellipsoïde approprié, dans ce cas, vous allez avoir à utiliser une projection réelle de la bibliothèque (tout ce qui est utilise proj4 derrière les coulisses, ces jours-ci) tels que les liaisons python pour GDAL/OGR ou (ce qui est beaucoup plus convivial) pyproj.

    Toutefois, si vous êtes d'accord avec une forme sphérique de la terre, qu'il est assez simple de le faire sans aucun des bibliothèques spécialisées.

    La plus simple de l'égalité de la zone de projection pour le calcul est un projection sinusoïdale. Fondamentalement, il suffit de multiplier la latitude de la longueur d'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

    Bon... Maintenant tout ce que nous avons à faire est de calculer l'aire de l'arbitraire d'un polygone dans un avion.

    Il y a un certain nombre de 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

    J'espère que ça va vous diriger dans la bonne direction, de toute façon...

    • Notez que les lignes sera déformé par une projection sinusoïdale. Il pourrait être sage pour interpoler le long de grands cercles entre les points, de sorte que les polygones ne sont pas déformés.
    • Bien sûr! J'essayais juste de montrer un simple rapprochement. Il n'est pas conçue pour être tout à fait exacte. La distorsion des côtés du polygone ne importe si il est en train de calculer un polygone avec très grands côtés et quelques verticies, cependant.
    • Je vous remercie beaucoup pour que la reprojection de la fonction, vous m'avez sauvé après une journée de perdue! Note pour les autres: après la projection, vous pouvez également utiliser Python shapely package pour calculer la zone avec shapely.geometry.Polygon(list(zip(x, y))).area pour obtenir la zone de la sq. mètres.
    InformationsquelleAutor Joe Kington
  3. 6

    Un peu tard peut-être, mais ici, c'est une autre méthode, à l'aide de Girard du théorème. Il indique que l'aire d'un polygone de grands cercles est de R**2 fois la somme des angles entre les polygones moins (N-2)*pi où N est le nombre de coins.

    Je pensais que ce serait la peine de poster, car il ne repose pas sur toutes les autres bibliothèques de numpy, et c'est une tout autre méthode que les autres. Bien sûr, cela ne fonctionne que sur une sphère, donc il y aura une certaine imprécision lors de l'application sur la Terre.

    Tout d'abord, je définir une fonction pour calculer l'angle d'azimut du point 1 le long d'un grand cercle 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 l'utiliser pour trouver les angles, et puis la région (Dans la suite, lons-le-saunier et le lats doit être spécifié, 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 le Colorado coordonnées données dans une autre réponse, et avec le rayon de la Terre 6371 km, j'obtiens que la zone est 268930758560.74808

    • J'ai essayé votre méthode et a obtenu une valeur négative. En outre, son abs (139013699.103) est assez différente de la valeur ( 809339.212 ), j'ai obtenu de l'aide d'une "grille" de la méthode, qui est de prendre toutes les cellules de la grille à l'intérieur du polygone à partir d'un bâtiment rectangulaire de Mercator de la grille, et en additionnant les cellules de grille domaines. La surface de chaque cellule est le produit de il zonale et méridienne intervalles (converti à partir lat,lon km).
    • Pourriez-vous clarifier votre commentaire? Où sont les numéros que vous référer à venir? J'ai essayé le code que j'ai posté à nouveau, et ça donne le résultat déclaré...
    • Je l'ai appliquée sur le contour trouvé à l'aide de matplotlib est contour fonction, en fait, j'ai quelques dizaines d'entre eux et tout le monde a signalé une zone négative à l'aide de votre méthode.
    • Ok. Je ne comprends toujours pas où les chiffres que vous mentionnez viennent. Je serais heureux de corriger toute erreur si vous pouvez être plus précis. Pouvez-vous poster un exemple de lons-le-saunier/lats qui donnent un négatif?
    • J'ai exporté un contour coordonnées ici pastebin.com/2fMkFgvu vous pourriez peut-être avoir un essai. Une chose que j'ai remarqué que si je sous-échantillon de 135 points de contour, en prenant toutes les 5 points (x=x[::5]; y=y[::5]), j'ai pu obtenir une valeur positive qui est proche de la "vraie" valeur.
    • Oui, je reçois zone négative trop avec votre exemple. Je suis assez sûr que la raison est que les erreurs d'arrondi lors du calcul des angles, depuis les points sont très proches. C'est pourquoi cela fonctionne si vous sous-échantillon. Il n'est donc pas une méthode robuste à haute précision, ce qui n'est pas très surprenant. Merci pour elle, mais

    InformationsquelleAutor sulkeh
  4. 4

    Parce que la terre est une surface fermée fermée polygone dessiné sur sa surface crée DEUX polygonale domaines. Vous devez également définir qui l'on est à l'intérieur et à l'extérieur!

    La plupart du temps les gens vont être aux prises avec des petits polygones, et donc il est "évident", mais une fois que vous avez des choses à la taille des océans et des continents, vous feriez mieux de vous assurer d'obtenir ce dans le bon sens.

    Aussi, n'oubliez pas que les lignes peuvent aller de l' (-179,0) à (+179,0) de deux façons différentes. L'un est beaucoup plus longue que l'autre. Encore une fois, surtout que vous allez faire l'hypothèse que c'est une ligne qui va de (-179,0) à (-180,0) qui est (+180,0) et ensuite (+179,0), mais un jour... il ne sera pas.

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

    • Très vrai! Je viens de donner un simple exemple de la plus basique des cas dans ma réponse... Ce n'est pas même à distance d'essayer de traiter avec passage au méridien d'origine (en 0-360) ou le antimeridian (de -180 à 180).
    • Ouais. La vie était tellement plus facile quand la terre était plate.
    InformationsquelleAutor Spacedman
  5. 3

    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.

    • J'aime que cette méthode me permet d'utiliser la latitude et la longitude des points directement, et les comptes de la courbure de la terre sans avoir besoin de projeter les points sur une grille Cartésienne. J'ai comparé les résultats de cette méthode et le standard Vert, le théorème de la zones des différents états sur Google et trouvé que cette bibliothèque a donné les résultats les plus précis. Merci!
    InformationsquelleAutor Ikar Pohorský
  6. 2

    Ici est 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 bien. NOTEZ que j'ai ajouté le 5ème point, de sorte que le chemin d'accès est 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.

    • Pour plus de précisions, quelle est la source pour que l'équation utilisée pour calculer l'aire?
    • C'est Vert du théorème.
    InformationsquelleAutor Jason