Geographic Information Algorithms

gis.py

# -*- coding: utf-8 -*-
import math
import re
from decimal import Decimal

# SAD-69 World Datum
A = 6378160
E2 = 0.00669454185


def deg2rad(deg, m=0.0, s=0.0):
    dd = Decimal(deg) + Decimal(m / 60.0) + Decimal(s / 3600.0)
    return Decimal(math.radians(dd))


def humanize(fn):
    def radfmt(rad):
        "Format radians to degrees, minutes and seconds"
        dd = Decimal(math.degrees(rad))
        d = int(dd)
        m = int((dd - d) * 60)
        s = (dd - Decimal(d) - Decimal(m) / Decimal(60)) * 3600
        return dd
        return u"%d:%d:%.4f" % (d, abs(m), abs(s))

    def process(*args, **kwargs):
        r = list(fn(*args, **kwargs))
        r = map(radfmt, r)
        return tuple(r)

    return process


@humanize
def cart2geo(x, y, z, precision=5):
    P = math.sqrt(x ** 2 + y ** 2)
    N = 1
    h = 0

    lng = math.atan2(y, x)
    print "lng = %f" % lng

    for i in range(precision):
        lat = math.atan(z / (P * (1 - E2 * N / (N + h))))
        print "lat = %f, N = %f, h = %f" % (lat, N, h)

        N = A / math.sqrt(1 - E2 * math.sin(lat) ** 2)
        h = P / math.cos(lat) - N

    return lat, lng


def geo2cart(lat, lng, h):
    lat = math.radians(lat)
    lng = math.radians(lng)

    N = A / math.sqrt(1 - E2 * math.sin(lat) ** 2)
    X = (N + h) * math.cos(lat) * math.cos(lng)
    Y = (N + h) * math.cos(lat) * math.sin(lng)
    Z = ((1 - E2) * N + h) * math.sin(lat)

    return X, Y, Z


def k_distortion(lat, lng, lng0):
    k0 = Decimal(0.9996)
    k = k0 / Decimal(math.sqrt(1 - (math.cos(lat) * math.sin(lng - lng0)) ** 2))

    return k


def distance(p1, p2):
    dx = Decimal(p1[0] - p2[0])
    dy = Decimal(p1[1] - p2[1])
    return Decimal(math.sqrt(dx ** 2 + dy ** 2))


def gps_area(pontos):
    """Returns the area given a list os GPS points (E, N).
            Area = 1/2 * sum(n^i * e^i+1) - sum(e^i * n^i+1)
    """
    soma1, soma2 = Decimal(0.0), Decimal(0.0)

    for i in range(len(pontos)):
        j = i + 1
        if j > len(pontos) - 1:
            j = 0

        soma1 += Decimal(pontos[i][1]) * Decimal(pontos[j][0])
        soma2 += Decimal(pontos[i][0]) * Decimal(pontos[j][1])

    return Decimal(0.5) * (soma1 - soma2)


print "\nConvert cartesian coordinates to geodesic:"
print cart2geo(4010210.546, -4260166.288, -2533008.133)

print "\nConvert geodesic coordinates to cartesian:"
print geo2cart(-20.58143311, -47.38050533, 1013.8147)

print "\nCalculate the UTM distortion (k) between 2 points:"
print "lat:", math.degrees(deg2rad(25,45,56))
print "lon:", math.degrees(deg2rad(46,41,49))
print k_distortion(deg2rad(25,45,56), deg2rad(46,41,49), deg2rad(45))

print "\nCalculate the area inside the given points:"
print gps_area([
    (323321.0,7393941.0),
    (323377.0,7393912.0),
    (323345.0,7393863.0),
    (323295.0,7393888.0),
])