Some GIS utils for Python, like area computation for geo-polygons

gisutils.py

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