Skip to content

Instantly share code, notes, and snippets.

@gka
Created October 6, 2011 14:03
Show Gist options
  • Save gka/1267470 to your computer and use it in GitHub Desktop.
Save gka/1267470 to your computer and use it in GitHub Desktop.
Some GIS utils for Python, like area computation for geo-polygons
"""
implementation taken from
http://forum.worldwindcentral.com/showthread.php?t=20724
algorithm expects a list of [lng,lat] pairs in degrees
"""
def haversine(x):
from math import cos
return ( 1.0 - cos(x) ) / 2.0
def area(pts, earthrad=6371):
from math import radians as rad, sin, cos, asin, sqrt, pi, tan, atan
pihalf = pi * .5
n = len(pts)
sum = 0
for j in range(0,n):
k = (j+1)%n
if j == 0:
lam1 = rad(pts[j][0])
beta1 = rad(pts[j][1])
cosB1 = cos(beta1)
else:
lam1 = lam2
beta1 = beta2
cosB1 = cosB2
lam2 = rad(pts[k][0])
beta2 = rad(pts[k][1])
cosB2 = cos(beta2)
if lam1 != lam2:
hav = haversine( beta2 - beta1 ) + cosB1 * cosB2 * haversine(lam2 - lam1)
a = 2 * asin(sqrt(hav))
b = pihalf - beta2
c = pihalf - beta1
s = 0.5 * (a+b+c)
t = tan(s*0.5) * tan((s-a)*0.5) * tan((s-b)*0.5) * tan((s-c)*0.5)
excess = abs(4*atan(sqrt(abs(t))))
if lam2 < lam1:
excess = -excess
sum += excess
return abs(sum)*earthrad*earthrad
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment