mercoledì 23 febbraio 2011

La mia prima volta... con python

Ho iniziato a studiare il Pyhton; mi sembra un passaggio obbligato per un utente GIS...: No more excuses to not use Python (James Fee)!

Così ho iniziato con il tutorial e con Dive Into Python.

Dopo qualche esercizio ho avuto bisogno di qualcosa di più concreto da fare: quale occasione migliore per ripescare nella memoria qualche confuso ricordo di cartografia?

Quel che segue è un brutto codice, sicuramente ridondante... benvenuti i suggerimenti e le correzioni!

Quasi tutte le fonti sono citate nei commenti.

#! /usr/bin/env python
# -*- coding: utf-8 -*-
from __future__ import  division
import math

class ellipsoid():
    """
        b semi-minor axis (polar)
        a semi-major axis (equatorial)
        f flattering
        rf reverse (inverse) flattering
        """
    def __init__(self, a,**kargs):    
        self.a=a
        self.b=None
        ## a sloppy and shallow parsing... 
        if len(kargs)<1:
            exit
        elif 'b' in kargs.keys():
            self.b = kargs['b']
        elif 'f' in kargs.keys():
            self.b =self.a*(1-kargs['f'])
        elif 'rf' in kargs.keys():
            self.b=self.a*(1-1/kargs['rf'])
            
        if 'name' in kargs.keys():
            self.name=kargs['name']
        else:
            self.name=None
            
    def getFirstEccentricity(self):
        """
        first eccentricity = sqrt(1-(b2/a2)**2)
        """
        return math.sqrt(1 - (self.b / self.a) ** 2)

    def getFirstEccentricitySquared(self):
        """
        first eccentricity squared = (a**2-b**2)/a**2
        """
        return (self.a ** 2 - self.b ** 2) / self.a ** 2

    def getSecondEccentricity(self):
        """
        second eccentricity = sqrt((a**2-b**2)/b**2)
        """
        return math.sqrt((self.a ** 2 - self.b ** 2) / self.b ** 2)

    def getSecondEccentricitySquared(self):
        """
        second eccentricity squared (a**2-b**2)/b***2
        """
        return (self.a ** 2 -self.b ** 2) / self.b ** 2

    def getFlattering(self):
        """
        flattering  =(a-b)/a
        """
        return (self.a - self.b) / self.a

    def getInverseFlattering(self):
        """
        inverse flattering =1/f = a/(a-b)
        """
        return 1 / self.getFlattering()

    def getNormal(self, lat):
        """
        prime vertical, normal at latitude lat
        """
        #        return self.a**2/math.sqrt( \    
        #                (self.a*math.cos(lat))**2 + \
        #                (self.b*math.sin(lat))**2)
        return self.a / math.sqrt(\
                1 - self.getFirstEccentricitySquared()*(math.sin(lat)) ** 2)

    def getMeridional(self, lat):
        """
        radius of curvature in the (north-south) meridian at latitude lat 
        """
        #-------------------------------- return (self.a*self.b)**2/math.pow(( \
        #--------------------------------- (self.a*math.cos(lat))**2 + \
        #--------------------------------- (self.b*math.sin(lat))**2), \
        #-------------------------------------------------------- (3/2))
        return (self.a * (1 - self.getFirstEccentricitySquared())) / (\
            math.pow(1 - self.getFirstEccentricitySquared()*(math.sin(lat)) ** 2, 3 / 2))

    def getPolarRadius(self):
        """
        Distance from the Earth's center to the North (or South) Pole
        """
        return self.b

    def getEquatorialRadius(self):
        """
        Distance from the Earth's center to the Equator
        """
        return self.a

    def getEllipsoidArea(self):
        """
        Surface area S of an oblate ellipsoid 
        (generated by an ellipse rotating around its minor axis)
        """
        #        formula from:
        #        http://www.numericana.com/answer/geometry.htm#ovalarea
        area=2 * math.pi * self.a ** 2 * (\
            1 + (self.a / self.b) ** 2 * \
            math.atan(self.getFirstEccentricity()) / \
            (self.getFirstEccentricity()))
        return area

    def getEllipsoidVolume(self):
        """
        Volume v of an oblate ellipsoid 
        (generated by an ellipse rotating around its minor axis)
        """
        vol=(4 / 3) * math.pi * self.a ** 2 * self.b
        return vol

    def getAuthalicRadius(self):
        """
        radius of a hypothetical perfect sphere which has 
        the same ellipsoid surface
        """
        #        formula from:
        #        http://en.wikipedia.org/wiki/Earth_radius#Notable_radii
        return math.sqrt(self.getEllipsoidArea() / (math.pi * 4))
    
    def getVolumetricRadius(self):
        """
        radius of a hypothetical perfect sphere which has 
        the same ellipsoid volume
        """
        #        formula from:
        #        http://en.wikipedia.org/wiki/Earth_radius#Notable_radii
        return math.pow(self.a**2*self.b, 1/3)

    def getGeodeticRadius(self, lat):
        """
        Distance from the Earth's center to a point on the 
        spheroid surface at geodetic latitude lat
        """
        #        formula from:
        #        http://en.wikipedia.org/wiki/Earth_radius#Notable_radii
        gradius= math.sqrt(\
            ((self.a**2 * math.cos(lat))**2 + (self.b**2 * math.sin(lat))**2) / \
            ((self.a * math.cos(lat))** 2 + (self.b * math.sin(lat))** 2))
        return gradius

    def getParallelRadius(self, lat):
        """
        radius of parallel at latitude lat
        """
        return self.getNormal(lat) * math.cos(lat)

    def getParallelArc(self, lat, fromLon, toLon):
        """
        arc of parallel between two longitude (fromLong toLon) at laitutde lat
        """
        return self.normal(lat) * math.cos(lat) * (toLon - fromLon)

    def getMeridianArcDelambre(self,lat):
        """
        arc of meridian from equator to lat
        Delambre series
        """
        e2, e4, e6, e8 = (math.pow(self.getFirstEccentricity(), i) for i in (2, 4, 6, 8))
        k = self.a * (1 - e2)
        A0 = k * (1 + 3 / 4 * e2 + 45 / 64 * e4 + \
                  175 / 256 * e6 + 11025 / 16384 * e8)
        A2 = -(k / 2) * (3 / 4 * e2 + 15 / 16 * e4 + \
                   525 / 512 * e6 + 2205 / 2048 * e8)
        A4 = k / 4 * (15 / 64 * e4 + 105 / 256 * e6 + 2205 / 4096 * e8)
        A6 = -(k / 6) * (35 / 512 * e6 + 315 / 2048 * e8)
        A8 = (k / 8) * (315 / 16348 * e8)
        return A0*lat+A2*math.sin(2*lat)+A4*math.sin(4*lat)+ \
                    A6*math.sin(6*lat)+ A8*math.sin(8*lat)

    def getMeridianArcBessel(self,lat):
        """
        arc of meridian from equator to lat
        Bessel, third flattering
        """
        n=(self.a-self.b)/(self.a+self.b)
        n2,n3,n4,n5=(math.pow(n, i) for i in range(2,6))
        D0=1+(9/4*n**2)+(255/64*n4)
        D2=(3/2*n)+(45/16*n3)+(525/128*n5)
        D4=(15/16*n2)+(105/64*n4)
        D6=(35/48*n3)+(315/256*n5)
        return self.a*(1-n)**2*(1+n)*((D0*lat)-(D2*math.sin(2*lat))+\
                    (D4*math.sin(4*lat))-(D6*math.sin(6*lat)))
    
    def getSferaLocale(self,lat):
        """raggio della sfera locale """
        return math.sqrt(self.getMeridional(lat)*self.getNormal(lat))

if __name__ == "__main__":
    '''qualche dato di prova'''
    bessel=ellipsoid(6377397.155,f=0.003342773154,name='bessel')
    wgs84 = ellipsoid(6378137.0, b=6356752.3142,name='wgs84')
    grs80=ellipsoid(6378137.0,rf=298.257222101,name='grs80')
    ell=[wgs84,bessel,grs80]
    x=bessel.getFlattering()
    lat=44
    lat_rad=math.radians(lat)
    print "Principali caratterisitche di alcuni ellissoidi"
    print "Quando viene richiesto, la latitudine è uguale a {0} gradi".format(lat)
    for e in ell:
        print "========================================================="
        print "Nome: {0}".format(e.name)
        print "Flattering \'f\': {0}".format(e.getFlattering())
        print "Inverse flatteirng \'rf\': {0}".format(e.getInverseFlattering())
        print 'Eccentricità'.rjust(30,'-')
        print "Eccentricità: {0}"\
            .format(e.getFirstEccentricity())
        print "Eccentricità al quadrato: {0}"\
            .format(e.getFirstEccentricitySquared())
        print "Eccentricità seconda: {0}"\
            .format(e.getSecondEccentricity())
        print "Eccentricità seconda al quadrato: {0}"\
            .format(e.getSecondEccentricitySquared())
        print 'Raggi'.rjust(30,'-')
        print "Raggio equatoriale \'a\': {0}".format(e.a)
        print "Raggio polare \'b\': {0}".format(e.b)
        print "Raggio geodetico: {0:.15}"\
            .format(e.getGeodeticRadius(lat_rad))
        print "Gran normale alla latitudine di {0} gradi: {1}"\
            .format(lat,e.getNormal(lat_rad))
        print "Raggio di curvatura del meridiano alla latitudine di {0} gradi: {1}"\
            .format(lat, e.getMeridional(lat_rad))
        print "Raggio di curvatura del parallelo alla latitudine di {0} gradi: {1}"\
            .format(lat, e.getParallelRadius(lat_rad))
        print "Raggio di una sfera di area equivalente: {0}"\
            .format(bessel.getAuthalicRadius())
        print "Raggio di una sfera di volume equivalente: {0}"\
            .format(bessel.getVolumetricRadius())
        print "Raggio della sfera locale alla  latitudine di {0} gradi: {1}"\
            .format(lat,bessel.getSferaLocale(lat_rad))
        print 'Archi di meridiano'.rjust(30,'-')
        print "Arco di meridiano dall'equatore alla latitudine di {0} gradi (Bessel): {1}"\
            .format(lat,e.getMeridianArcBessel(lat_rad))
        print "Arco di meridiano dall'equatore alla latitudine di {0} gradi (Delambre): {1}"\
            .format(lat,e.getMeridianArcDelambre(lat_rad))