VSOP2013 e nuove routine di precisione

Ho scoperto che il gruppo francese dell'Observatoire de Paris ha continuato a elaborare le teorie semianalitiche per il calcolo della posizione dei pianeti, fino ad una edizione del 2013, di elevatissima precisione, di cui gli autori parlano in questo articolo.

Le VSOP87, di cui ho presentato una piccola libreria personale in linguaggio Python, sono ancora molto diffuse. Il lavoro degli studiosi francesi ha pero' prodotto delle novità nel corso degli anni, producendo una versione 2010 e quella di cui parlo oggi, aggiornata al 2013.

Tra le novità più salienti, l'inclusione nella teoria del pianeta Plutone, il calcolo del baricentro Terra_Luna, e una teoria apposita per i pianeti esterni.

Inoltre la precisione di calcolo si estende nel tempo ad un intervallo che va dal 4000 a.C. all'8000 d.C. con mantenimento di elevata accuratezza.

Tutto il materiale necessario per il calcolo è fornito nel sito ftp: ftp://ftp.imcce.fre.

Per gli esempi contenuti in questo post vi consiglio di scaricare l'intera directory ftp://ftp.imcce.fr/pub/ephem/planets/vsop2013/solution/ che contiene anche le tabelle dati per il calcolo, molto voluminose per l'inclusione di decine di migliaia di termini correttivi che tengono conto di tutte le perturbazioni al moto prodotte dai corpi celesti del sistema solare.

Il programma che presento di seguito è, di fatto, una traduzione in python del programma in linguaggio FORTRAN allegato dagli autori, che serve a produrre una tabella di risultati per 11 date a partire dal 26 giugno 1890 ore 12 e intervallate fra loro di 4000 giorni Per garantire la precisione del risultato ho incluso la libreria mpmath che consente l'esecuzione di calcoli in elevata precisione. Se confrontate la tabella prodotta a video dal mio programma con il file VSOP2013.ctl, constaterete un'elevata corrispondenza, che arriva fino alla decima cifra decimale se si imposta una precisione pari ad almeno 1e-12. Con il crescere della precisione crescono esponenzialmente i tempi di calcolo, per l'inclusione di un sempre maggiore numero di termini. Operativamente, per le nostre necessità, una precisione di 1e-7 puo' essere ritenuta sufficiente ed è impostata di default nella classe che vi presento di seguito.

import mpmath as mpm
from mpmath import mp

import struct

mp.dps = 50
mp.pretty = True

def cal2jul(year, month, day, hour=0, minute =0, second=0):
    month2 = month
    year2 = year
    if month2 <= 2:
        year2 -= 1
        month2 += 12
    if (year*10000 + month*100 + day) > 15821015:
        a = int(year2/100)
        b = 2 - a + int(a/4)
    else:
        a = 0
        b = 0
    if year < 0:
        c = int((365.25 * year2)-0.75)
    else:
        c = int(365.25 * year2)
    d = int(30.6001 *(month2 + 1))
    return b + c + d + day + hour / 24.0 + minute / 1440.0 + second / 86400.0 + 1720994.5        




class VSOP2013():
    
    def __init__(self, t, planet, precision=1e-7):
        # calculate millennia from J2000
        self.JD = t
        self.t = mpm.fdiv((t - cal2jul(2000,1,1,12)), 365250.0)
        # predefine powers of self.t
        self.power = []; self.power.append(mpm.mpf('1.0')); self.power.append(self.t)
        for i in xrange(2,21):
            t = self.power[-1]
            self.power.append(mpm.fmul(self.t,t))
        # choose planet file in a dict
        self.planet = planet
        self.planets = {'Mercury':'VSOP2013p1.dat',
                        'Venus'  :'VSOP2013p2.dat',
                        'EMB'    :'VSOP2013p3.dat',
                        'Mars'   :'VSOP2013p4.dat',
                        'Jupiter':'VSOP2013p5.dat',
                        'Saturn' :'VSOP2013p6.dat',
                        'Uranus' :'VSOP2013p7.dat',
                        'Neptune':'VSOP2013p8.dat',
                        'Pluto'  :'VSOP2013p9.dat'}
        # VSOP2013 routines precision
        self.precision = precision
        
        # lambda coefficients
        # l(1,13) : linear part of the mean longitudes of the planets (radian).
        # l(14): argument derived from TOP2013 and used for Pluto (radian).
        # l(15,17) : linear part of Delaunay lunar arguments D, F, l (radian).

        self.l = (
            (mpm.mpf('4.402608631669'), mpm.mpf('26087.90314068555')),
            (mpm.mpf('3.176134461576'), mpm.mpf('10213.28554743445')),
            (mpm.mpf('1.753470369433'), mpm.mpf(' 6283.075850353215')),
            (mpm.mpf('6.203500014141'), mpm.mpf(' 3340.612434145457')),
            (mpm.mpf('4.091360003050'), mpm.mpf(' 1731.170452721855')),
            (mpm.mpf('1.713740719173'), mpm.mpf(' 1704.450855027201')), 
            (mpm.mpf('5.598641292287'), mpm.mpf(' 1428.948917844273')),
            (mpm.mpf('2.805136360408'), mpm.mpf(' 1364.756513629990')),
            (mpm.mpf('2.326989734620'), mpm.mpf(' 1361.923207632842')),
            (mpm.mpf('0.599546107035'), mpm.mpf('  529.6909615623250')),
            (mpm.mpf('0.874018510107'), mpm.mpf('  213.2990861084880')),
            (mpm.mpf('5.481225395663'), mpm.mpf('   74.78165903077800')),
            (mpm.mpf('5.311897933164'), mpm.mpf('   38.13297222612500')),
            (mpm.mpf('0.000000000000'), mpm.mpf('    0.3595362285049309')),
            (mpm.mpf('5.198466400630'), mpm.mpf('77713.7714481804')),
            (mpm.mpf('1.627905136020'), mpm.mpf('84334.6615717837')),
            (mpm.mpf('2.355555638750'), mpm.mpf('83286.9142477147'))
            )

        # planetary frequencies in longitude
        self.freqpla = {'Mercury'   :   mpm.mpf('0.2608790314068555e5'),  
                        'Venus'     :   mpm.mpf('0.1021328554743445e5'),
                        'EMB'       :   mpm.mpf('0.6283075850353215e4'),
                        'Mars'      :   mpm.mpf('0.3340612434145457e4'), 
                        'Jupiter'   :   mpm.mpf('0.5296909615623250e3'),  
                        'Saturn'    :   mpm.mpf('0.2132990861084880e3'),
                        'Uranus'    :   mpm.mpf('0.7478165903077800e2'),  
                        'Neptune'   :   mpm.mpf('0.3813297222612500e2'), 
                        'Pluto'     :   mpm.mpf('0.2533566020437000e2')}  

        
            
        # target variables
        self.ax = mpm.mpf('0.0')            # major semiaxis
        self.ml = mpm.mpf('0.0')            # mean longitude
        self.kp = mpm.mpf('0.0')            # e*cos(perielium longitude)
        self.hp = mpm.mpf('0.0')            # e*sin(perielium longitude)
        self.qa = mpm.mpf('0.0')            # sin(inclination/2)*cos(ascending node longitude)
        self.pa = mpm.mpf('0.0')            # sin(inclination/2)*cos(ascending node longitude)
        
        self.tg_var = {'A':self.ax, 'L':self.ml, 'K':self.kp,
                       'H':self.hp, 'Q':self.qa, 'P':self.pa } 

        # eps  = (23.d0+26.d0/60.d0+21.41136d0/3600.d0)*dgrad
        self.eps = mpm.mpf((23.0+26.0/60.0+21.411360/3600.0)*mpm.pi/180.0)
        self.phi  = mpm.mpf(-0.051880 * mpm.pi / 180.0 / 3600.0)
        self.ceps = mpm.cos(self.eps)
        self.seps = mpm.sin(self.eps)
        self.cphi = mpm.cos(self.phi)
        self.sphi = mpm.sin(self.phi)
        
        # rotation of ecliptic -> equatorial rect coords
        self.rot = [[self.cphi, -self.sphi*self.ceps,  self.sphi*self.seps],
                    [self.sphi,  self.cphi*self.ceps, -self.cphi*self.seps],
                    [0.0,        self.seps,            self.ceps          ]]
        
        self.fmt = struct.Struct("""6s 3s 3s 3s 3s x 3s 3s 3s 3s 3s x 4s 4s 4s 4s x
                       6s x 3s 3s 3s 20s x 3s 20s x 3s x""")

        self.gmp = {
       'Mercury' : mpm.mpf('4.9125474514508118699e-11'),
       'Venus'   : mpm.mpf('7.2434524861627027000e-10'),
       'EMB'     : mpm.mpf('8.9970116036316091182e-10'),
       'Mars'    : mpm.mpf('9.5495351057792580598e-11'),
       'Jupiter' : mpm.mpf('2.8253458420837780000e-07'),
       'Saturn'  : mpm.mpf('8.4597151856806587398e-08'),
       'Uranus'  : mpm.mpf('1.2920249167819693900e-08'),
       'Neptune' : mpm.mpf('1.5243589007842762800e-08'),
       'Pluto'   : mpm.mpf('2.1886997654259696800e-12')
       }

        self.gmsol = mpm.mpf('2.9591220836841438269e-04')
        self.rgm = mpm.sqrt(self.gmp[self.planet]+self.gmsol)
        # run calculus routine
        self.calc()
        

    def __str__(self):
        
        vsop_out = "{:3.13} {:3.13} {:3.13} {:3.13} {:3.13} {:3.13}\n".format(
                                          self.tg_var['A'],
                                          self.tg_var['L'],
                                          self.tg_var['K'],
                                          self.tg_var['H'],
                                          self.tg_var['Q'],
                                          self.tg_var['P'])
        vsop_out += "{:3.13} {:3.13} {:3.13} {:3.13} {:3.13} {:3.13}\n".format(
                                          self.ecl[0],
                                          self.ecl[1],
                                          self.ecl[2],
                                          self.ecl[3],
                                          self.ecl[4],
                                          self.ecl[5])
        vsop_out += "{:3.13} {:3.13} {:3.13} {:3.13} {:3.13} {:3.13}\n".format(
                                          self.equat[0],
                                          self.equat[1],
                                          self.equat[2],
                                          self.equat[3],
                                          self.equat[4],
                                          self.equat[5])
                        

        return vsop_out



    def calc(self):
        with open(self.planets[self.planet]) as file_in:
            terms = []
            b = '*'
            while b != '':
                b = file_in.readline()
                if b != '':
                    if b[:5] == ' VSOP':
                        header = b.split()
                        #print header[3], header[7], header[8], self.t**int(header[3])
                        no_terms = int(header[4])
                        for i in xrange(no_terms):
                            #6x,4i3,1x,5i3,1x,4i4,1x,i6,1x,3i3,2a24
                            terms = file_in.readline()
                            a = self.fmt.unpack(terms)
                            S = mpm.fmul(a[18],mpm.power(10,int(a[19])))
                            C = mpm.fmul(a[20],mpm.power(10,int(a[21])))
                            if mpm.sqrt(S*S+C*C) < self.precision:
                                break
                            aa = 0.0; bb = 0.0;
                            for j in xrange(1,18):
                                aa += mpm.fmul(mpm.mpf(a[j]), self.l[j-1][0])
                                bb += mpm.fmul(mpm.mpf(a[j]), self.l[j-1][1])
                            arg = aa + bb * self.t
                            power = int(header[3])
                            comp = self.power[power] * (S * mpm.sin(arg) + C * mpm.cos(arg))
                            if header[7] == 'L' and power == 1 and int(a[0]) == 1:
                                pass
                            else:
                                self.tg_var[header[7]] += comp
        self.tg_var['L'] = self.tg_var['L'] + self.t * self.freqpla[self.planet]
        self.tg_var['L'] = self.tg_var['L'] % (2 * mpm.pi)
        if self.tg_var['L'] < 0:
            self.tg_var['L'] += 2*mpm.pi
        print "Julian date {}".format(self.JD)
        file_in.close()
        ##print self.tg_var
    
        ####    def ELLXYZ(self):

        xa = self.tg_var['A']
        xl = self.tg_var['L']
        xk = self.tg_var['K']
        xh = self.tg_var['H']
        xq = self.tg_var['Q']
        xp = self.tg_var['P']
        
        # Computation
  
        xfi = mpm.sqrt(1.0 -xk * xk - xh * xh)
        xki = mpm.sqrt(1.0 -xq * xq - xp * xp)
        u = 1.0/(1.0 + xfi)
        z = complex(xk, xh)
        ex = abs(z)
        ex2 = ex  * ex
        ex3 = ex2 * ex
        z1 = z.conjugate()
        #
        gl = xl % (2*mpm.pi)
        gm = gl - mpm.atan2(xh, xk)
        e  = gl + (ex - 0.1250 * ex3) * mpm.sin(gm)
        e += 0.50 * ex2 * mpm.sin(2.0 * gm)
        e += 0.3750 * ex3 * mpm.sin(3.0 * gm)
        #
        while True:
            z2 = complex(0.0, e)
            zteta = mpm.exp(z2)
            z3 = z1 * zteta
            dl = gl - e + z3.imag
            rsa = 1.0 - z3.real
            e = e + dl / rsa
            if abs(dl) < 1e-15:
                break
        #
        z1  =  u * z * z3.imag
        z2  =  complex(z1.imag, -z1.real)
        zto = (-z+zteta+z2)/rsa
        xcw = zto.real
        xsw = zto.imag
        xm  = xp * xcw - xq * xsw
        xr  = xa * rsa
        #
        self.ecl = []; self.equ = {}
        self.ecl.append(xr * (xcw -2.0 *xp * xm))
        self.ecl.append(xr * (xsw +2.0 *xq * xm))
        self.ecl.append(-2.0 * xr * xki * xm)
        #
        xms = xa *(xh + xsw) / xfi
        xmc = xa *(xk + xcw) / xfi
        xn  = self.rgm / xa ** (1.50)
        #
        self.ecl.append( xn *((2.0 * xp * xp - 1.0) * xms + 2.0 * xp * xq * xmc))
        self.ecl.append( xn *((1.0 -2.0 * xq * xq) * xmc -2.0 * xp * xq * xms))
        self.ecl.append( 2.0 * xn * xki * (xp * xms + xq * xmc))
        
        # Equatorial rectangular coordinates and velocity

        #         
        #
        # --- Computation ------------------------------------------------------
        #
        self.equat = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]    
        for i in xrange(3):
           for j in xrange(3):
               self.equat[i]   = self.equat[i]   + self.rot[i][j] * self.ecl[j]
               self.equat[i+3] = self.equat[i+3] + self.rot[i][j] * self.ecl[j+3]



if __name__ == '__main__':
    for planet in ('Mercury', 'Venus', 'EMB', 'Mars', 'Jupiter',
          'Saturn', 'Uranus', 'Neptune', 'Pluto'):
        print "PLANETARY EPHEMERIS VSOP2013  "+ planet + "(TDB)\n"+""" 
  1/ Elliptic   Elements:    a (au), lambda (radian), k, h, q, p        - Dynamical Frame J2000
  2/ Ecliptic   Heliocentric Coordinates:  X,Y,Z (au)  X',Y',Z' (au/d)  - Dynamical Frame J2000
  3/ Equatorial Heliocentric Coordinates:  X,Y,Z (au)  X',Y',Z' (au/d)  - ICRS Frame J2000
"""
        init_date = cal2jul(1890,6,26,12)
        set_date = init_date

        while set_date < init_date + 41000:
            v = VSOP2013(set_date, planet, 1e-12)
            print v
            set_date += 4000
          
        

Il programma accetta in input tre parametri:

1. la data giuliana comprensiva dell'ora in frazione decimale di giorno
2. il nome del pianeta richiesto (in inglese con l'iniziale maiuscola, vedere la lista
3. la precisione richiesta, se non è inserita scatta quella di default (1e-7)

Il programma è in ogni caso lento nell'esecuzione. Nel prossimo post cerchero' di sfruttare qualche tecnica per migliorare la velocità esecutiva, ma dubito di poter lavorare in Python alla velocità di Fortran. Al prossimo post.

Ancora un utile widget, per gestire informazioni geografiche - parte 2

Ora che il database sqlite è stato realizzato, possiamo dedicarci al nuovo widget. Per questo piccolo progetto utilizzo ancora la libreria Tkinter, disponibile nell'installazione di base di Python (nulla vieta di usare altre librerie, che sicuramente offrono un'interfaccia meno scarna, come le wx o le qt, ma Tkinter è facile da usare e pronto da subito, per cui lo preferisco).

L'utilità di questo widget in se è, come vedremo, limitata, ma si presta ottimamente all'uso nell'ambito di un progetto più ampio. Ne verificheremo in seguito una possibile applicazione.

Scopo del widget è l'estrazione, a richiesta dell'utilizzatore, delle coordinate geografiche e della timezone di un qualsiasi luogo abitato della superficie terrestre, presente nel database localita.db, unitamente alle informazioni ulteriori che possono essere utili per risolvere omonimie e fornire dettagli utili.

L'immagine a sinistra puo' dare un'idea del funzionamento del widget, ospitato, insieme al widget calendario, di cui ho presentato il progetto in un precedente post, in un frame più grande.

Il funzionamento è molto semplice: è sufficiente inserire il nome completo della località o alcune lettere consecutive, anche interne al nome, nella casella indicata dal pulsante <-- Substr Search. in tempi molto rapidi la combobox sottostante sarà popolata di tutti i nomi delle località che soddisfano il criterio di ricerca. Di default viene presentata la prima della lista, ma si puo', ovviamente, selezionarne una diversa tra quelle che si presentano quando si clicca sulla freccetta a destra. Nella sottostante casella arancione vengono elencati, in ordine, i campi corrispondenti alla selezione. (I colori sono, naturalmente, configurabili).

Passiamo all'esame del listato, ho chiamato lo script python tkAtlas.py (lo trovate, come sempre, in github, sotto la directory atlas, che contiene tutti gli script per creare il database. Come al solito, riporto il codice in frammenti per commentarlo adeguatamente. Il primo frammento comprende l'inizio e la fine dello script, gli altri sono inseriti in sequenza dopo le prime righe.

from Tkinter import *
import tkFont
import ttk
import sqlite3
import os

class Cities:
    """
    This class is a new widget, aimed to deal with world cities from geoname site
    (see http://download.geonames.org/export/dump/) for details.
    A lot of useful informations are taken from a main archive (allCountries.txt)
    and linked to textual geographical informations from the same site
    (admin1CodesASCII.txt and admin2Codes.txt), then reversed in a sqlite3 db,
    which is read from inside the present program
    """

    def __init__(self, master):
        """
        All the widgets are created in __init__, all significant variables
        are instance variables and can be read from outside
        """
        self.master = master

----------


if __name__ == '__main__':
    root = Tk()
    app = Cities(root)
    root.mainloop()

Lo schema dovrebbe esservi ormai familiare: importazione di moduli, dichiarazione di classe, docstring e metodo __init__ per l'inizializzazione dell'istanza di classe, richiesta, in questo caso, dalle ultime righe dello script, eseguite solo se lo script viene lanciato direttamente anzichè essere importato come modulo da un altro script.

        root_dir = os.path.dirname(os.getcwd())
        for root, dirs, files in os.walk(root_dir, topdown=True):
            for file in files:
                if file.endswith('localita.db'):
                    self.db = os.path.join(root, file)

In queste righe si effettua la ricerca del database. Devo dire che la libreria standard di Python non rende molto facile l'esplorazione dei path, per cui, non volendo usare librerie di terze parti, ho dovuto lavorarci un po'.

Se il database è nella stessa directory dello script principale o in una qualsiasi sottodirectory, viene recuperato il path completo del file e viene messo a disposizione della successiva routine di query sql.

        self.conn = sqlite3.connect(self.db)
        self.c = self.conn.cursor()

Con queste righe si gestisce la connessione ad database sqlite. Ricordo che i metodi per la gestione dei database sqlite sono già incorporati nella libreria standard di Python, è sufficiente importare il modulo relativo con l'espressione import sqlite3, che abbiamo collocato all'inizio dello script.

        self.frame1 = Frame(self.master, background='#006f6f', padx=5, pady=5)
        self.frame1.pack(fill=X, expand=1)
        self.label1 = Label(self.frame1, text='Location')
        self.label1.pack(side=LEFT, anchor='w', fill=X, expand=1)
        self.search = StringVar(self.master)
        self.search.set('')
        self.result = ''
        self.entry1 = Entry(self.frame1, textvariable=self.search)
        self.entry1.pack(side=LEFT, fill=X, expand=1)
        self.button1 = Button(self.frame1, text='<-- Substr. Search', command=self.select)
        self.button1.pack(side=LEFT, fill=X, expand=1)
        self.frame2 = Frame(self.master)
        self.frame2.pack()
        self.combo1 = ttk.Combobox(self.frame2, text=self.result, width=60)
        self.combo1.pack(anchor='w')
        self.combo1.bind('<<ComboboxSelected>>', self.update_label)
        self.label2 = Label(self.frame2, text='\n\n\n\n\n\n\n\n\n', justify=LEFT, anchor='w', width=60, relief=GROOVE, background='orange')
        self.label2.pack(side=LEFT)

In queste righe di codice si creano vari oggetti (widget) per la visualizzazione grafica delle operazioni di query del database. Ho previsto due frame, il primo contiene una label, una entry e un button, il secondo una combobox e una label. La casella di entry è quella che utilizziamo per la stringa di ricerca che, come ho detto sopra, puo' essere il nome completo della località o una sottostringa. Maiuscole e minuscole non sono importanti: SQL è un linguaggio case-insensitive. La variabile wrapper self.search non è usata per il binding di un metodo, ma solo come contenitore della sottostringa di ricerca.

Il button invece, quando viene cliccato, lancia il metodo select che vediamo successivamente. La combobox non è tra i widget standard di Tkinter ma è contenuta nel modulo ttk, che ne è un'estensione, già comunque presente nell'installazione Python di base. Più interessante è il binding che faccio tra la scelta di un item della combobox e il lancio del metodo update_label. Il binding fatto in questo modo puo' essere scelto anche per tutti i casi di un elemento widget non responsivo al click del mouse o alla pressione della tastiera, per esempio il passaggio del mouse su una label o su un'immagine.

Infine ho previsto per la seconda label (quella in arancione) una disposizione incolonnata delle informazioni relative alla località, per ragioni di spazio e di leggibilità. All'inizio è composta di sole righe vuote, viene aggiornata con il metodo che vediamo successivamente.

    def select(self, *args):
        """
        This is the search engine in the database. As a result, the combobox
        and the bottom label texts are updated
        """
        search = self.search.get()
        p_query = ('%' + search + '%',)
        query = 'select nome_ascii, latitudine, longitudine, admin1, country, ascad1, ascad2, popolazione, elevazione, timezone                 from localita where nome_ascii like ? order by nome_ascii'
        self.c.execute(query, p_query)
        self.cities_list = self.c.fetchall()
        self.combo1.config(values='')
        values = [ [x[0],
         x[4],
         x[5],
         x[6]] for x in self.cities_list ]
        for i in range(len(values)):
            for j in range(4):
                values[i][j] = str(values[i][j])

            values[i] = tuple(values[i])

        self.combo1['values'] = tuple(values)
        self.combo1.set(values[0][0] + ' ' + values[0][1] + ' ' + values[0][2])
        index = self.combo1.index(0)
        self.update_label(index)

Questo metodo della classe Cities crea la query in base all'input utente, la esegue e raccoglie in risultati in una lista, che puo' essere esplorata con i metodi consueti. Nella combobox sono riportati solo alcuni valori estratti, cioè il nome in caratteri ascii della località. il paese e la regione (o la principale divisione del territorio secondo le regole del paese. La label invece viene popolata con tutti campi, come si vede nel codice seguente. Per ovviare alla possibile mancanza di informazioni (risultato None della query) ho preventivamente convertito in stringa tutti i risultati, None compreso, per poi rituplizzare il tutto, visto che la combobox accetta solo tuple.

    def update_label(self, *args):
        """
        This routine updates the index value to get back the selected row
        in the query results and subsequently updates the top down label
        """
        if args[0] == 0:
            index = 0
        else:
            index = self.combo1.current()
        _list = list(self.cities_list[index])
        s = [ '{}'.format(x) for x in _list ]
        long_list = 'Name: {0}\nLat: {1}\nLong: {2}\nState: {3} ({4})\n'
        long_list += 'Region: {5}\nProvince: {6}\nPopulation: {7}\n'
        long_list += 'Elevation: {8}\nTimezone: {9}'
        self.label2['text'] = long_list.format(s[0], s[1], s[2], s[3], s[4], s[5], s[6], s[7], s[8], s[9])
        self.parameters = self.cities_list[index]

L'ultimo frammento di codice è relativo alla routine (metodo) di aggiornamento in base ai risultati della query. Come si puo' vedere i parametri di ingresso possono essere in numero variabile, a me interessava che, nel momento in cui scelgo l'indice della combobox, come al termine del metodo select, la routine possa gestire l'informazione. In questo caso l'indice al termine della query è sempre fissato a 0, ma cambia in base alla scelta dell'utente sulla lista presente nella combobox. Nel primo caso l'assegnazione è automatica nel corso dell'esecuzione del metodo select, nel secondo consegue alla scelta fatta dall'utente con il successivo lancio della funzione di callback update_label. Questo metodo è pensato essenzialmente per riempire la label con le informazioni disponibili dopo l'esecuzione della query. Ho pensato, inoltre, di mettere a disposizione questi risultati inserendoli in una lista che potrà essere letta al di fuori della classe, come variabile di istanza. In sostanza i valori contenuti nella riga scelta confluiscono in self.parameters e sono leggibili come variabile di istanza Cities().parameters o, come in questo caso, come app.parameters. Da qui in avanti iniziamo a scrivere un vero e proprio software astrologico, in cui il widget calendario e il widget atlante saranno strumenti fondamentali per la gestione dell'ora e del luogo degli eventi astrologici.

Ancora un utile widget, per gestire informazioni geografiche - parte 1

Su questo blog ho già parlato della creazione di un database sqlite per gestire le informazioni geografiche che sono utili al lavoro dell'astrologo. Riprendo il tema, ma stavolta associandolo alla creazione di un altro widget complesso, che possiamo utilizzare in vari contesti, sia autonomo sia all'interno di un frame complesso.

Partiamo dal sito geonames che ci fornisce il materiale di base. I file contenuti veongono aggiornati spesso, il che ci aiuta a cogliere eventuali differenze demografiche e di geografia politica. Nello specifico ho selezionato i seguenti:

1. allCountries.zip
2. admin1CodesASCII.txt
3. admin2Codes.txt

Il primo file, molto voluminoso (più di 10 milioni di righe), contiene le coordinate geografiche di ogni luogo abitato del mondo, nonché elevazione, popolazione, timezone e codici relativi alle suddivisioni amministrative di appartenenza, i cui nomi per esteso sono contenuti nei due file successivi.

Per l'uso pratico del database delle località ho scelto di selezionare solo quei luoghi per cui la popolazione censita è di almeno 1 abitante, il che riduce considerevolmente il numero di righe che andranno a costituire il database sqlite finale. Nel seguito allego i file python che consentono di generare delle tabelle testuali per selezioanre e organizzare le informazioni esistenti.

Iniziamo dalla prima tabella testuale (file crea_tab1.py):

#1/usr/bin/env python
file_in = open('allCountries.txt','r')
file_out = open('all_gen.txt','w')
elenco=file_in.readlines()
file_in.close()
count = 1
for i in elenco:
    a=i.split('\t')
    stringa ="%s;%s;%s;%s;%s;%s;%s;%s;%s;%s;%s;%s;%s"
    # print a[1], a[2], a[4], a[5], a[8], a[10], a[14], a[15], a[16]
    if int(a[14])>100:
        # for j in (1, 2, 4, 5, 8, 10, 14, 15, 17, 18):
        file_out.write(stringa % (a[1], a[2], a[4], a[5], a[8], a[10],
                                  a[11], a[12], a[13], a[14],
                                  a[16], a[17], a[18]))
file_out.close()
file_in.close()

##00 geonameid         : integer id of record in geonames database
##01 name              : name of geographical point (utf8) varchar(200)
##02 asciiname         : name of geographical point in plain ascii characters, varchar(200)
##03 alternatenames    : alternatenames, comma separated, ascii names automatically transliterated, convenience attribute from alternatename table, varchar(10000)
##04 latitude          : latitude in decimal degrees (wgs84)
##05 longitude         : longitude in decimal degrees (wgs84)
##06 feature class     : see http://www.geonames.org/export/codes.html, char(1)
##07 feature code      : see http://www.geonames.org/export/codes.html, varchar(10)
##08 country code      : ISO-3166 2-letter country code, 2 characters
##09 cc2               : alternate country codes, comma separated, ISO-3166 2-letter country code, 60 characters
##10 admin1 code       : fipscode (subject to change to iso code), see exceptions below, see file admin1Codes.txt for display names of this code; varchar(20)
##11 admin2 code       : code for the second administrative division, a county in the US, see file admin2Codes.txt; varchar(80) 
##12 admin3 code       : code for third level administrative division, varchar(20)
##13 admin4 code       : code for fourth level administrative division, varchar(20)
##14 population        : bigint (8 byte int) 
##15 elevation         : in meters, integer
##16 dem               : digital elevation model, srtm3 or gtopo30, average elevation of 3''x3'' (ca 90mx90m) or 30''x30'' (ca 900mx900m) area in meters, integer. srtm processed by cgiar/ciat.
##17 timezone          : the timezone id (see file timeZone.txt) varchar(40)
##18 modification date : date of last modification in yyyy-MM-dd format

il file allNews.txt è un file testuale in cui ogni riga è relativa ad una località, e contiene il nome in caratteri unicode, in caratteri ascii, vari codici relativi ai paesi e alle maggiori suddivisioni amministrative, popolazione, elevazione e timezone.

Il secondo script python (crea_tab2.py) genera nomi e codici delle principali suddivisioni amministrative, rilevate dal file admin1CodesASCII.txt:

#1/usr/bin/env python
file_in = open('admin1CodesASCII.txt','r')
file_out = open('all_admin1.txt','w')
elenco=file_in.readlines()
file_in.close()

for i in elenco:
    a=i.split('\t')
    b = a[0].split('.')
    stringa = "{};{};{};{};{}"
    file_out.write(stringa.format(b[0],b[1],a[1],a[2],a[3]))

file_out.close()
file_in.close()

Analogamente il terzo script (crea_tab3.py) per le ulteriori suddivisioni amministrative:

#1/usr/bin/env python
file_in = open('admin2Codes.txt','r')
file_out = open('all_admin2.txt','w')
elenco=file_in.readlines()
file_in.close()

for i in elenco:
    a=i.split('\t')
    b = a[0].split('.')
    stringa = "{};{};{};{};{};{}"
    file_out.write(stringa.format(b[0],b[1],b[2],a[1],a[2],a[3]))

file_out.close()
file_in.close()

Infine il quarto script (crea_tab4.py) genera una tabella con le sigle dei paesi e relativo nome per esteso, secondo le codifiche ISO 3166:

file_in = open("ISO 3166 Codes (Countries).csv")
                 
cnt = file_in.readlines()
lista = []
for i in cnt:
    lista.append(i.split(';'))
file_out = open('all_countries.txt', 'w')
for i in lista:
    file_out.write("{};{}\n".format(i[0].strip(),i[1].strip()))
file_out.close()
file_in.close()

L'elenco delle denominazioni ISO 3166 è facilmente reperibile sul web. Per ottenere un file sorgente pulito ho ricombinato il tutto in un file csv.

Con le quattro tabelle testuali appena create posso costruire finalmente il database sqlite localita.db. Allo scopo ho creato uno script che va lanciato all'atto della creazione del database con il comando in linea: sqlite3 localita.db <crea_db.txt. Il programma eseguibile sqlite3 non è preinstallato nella mia distribuzione ubuntu, ma si trova nel repository. Quindi, come di consueto, sarà sufficiente il comando sudo apt-get install sqlite3.

Allego il contenuto dello script crea_db.txt:

SELECT datetime('now')||'Luoghi';
DROP TABLE IF EXISTS luoghi;
CREATE TABLE luoghi (
   nome_unicode text, 
   nome_ascii text, 
   latitudine real, 
   longitudine real, 
   admin1 text, 
   admin2 text, 
   admin3 text, 
   admin4 text, 
   regione text, 
   popolazione int,
   elevazione int,
   timezone text,
   mod_date text);
SELECT datetime('now')||'a1';
DROP TABLE IF EXISTS a1;
CREATE TABLE a1 (
   admin1 text, 
   admin2 text,
   uniad1 text,
   ascad1 text,
   codice1 int);
SELECT datetime('now')||'a2';
DROP TABLE IF EXISTS a2; 
CREATE TABLE a2 (
   admin1 text,
   admin2 text,
   admin3 text,
    uniad2 text,
   ascad2 text,
   codice2 int);
SELECT datetime('now')||'a3';
DROP TABLE IF EXISTS a3; 
CREATE TABLE a3 (
   country text,
   admin1 text);
.separator ";"
.import all_gen.txt luoghi
.import all_admin1.txt a1
.import all_admin2.txt a2
.import all_countries.txt a3

SELECT datetime('now')||'loc1';
DROP TABLE IF EXISTS loc1;
CREATE TABLE loc1 AS SELECT * from luoghi as L
     LEFT JOIN a1 as A on L.admin1 = A.admin1 and  L.admin2 = A.admin2;

SELECT datetime('now')||'localita';
DROP TABLE IF EXISTS localita;
CREATE TABLE localita AS SELECT nome_ascii, nome_unicode, latitudine, longitudine, L.admin1, B.country, L.admin2, popolazione, elevazione, timezone, mod_date, uniad1, ascad1, uniad2, ascad2 from loc1 as L
   LEFT JOIN a2 as A on L.admin1 = A.admin1 and L.admin2 = A.admin2 and L.admin3 = A.admin3
   LEFT JOIN a3 as B on L.admin1 = B.admin1;
SELECT datetime('now');
   DROP TABLE luoghi;
   DROP TABLE a1;
   DROP TABLE a2;
   DROP TABLE a3;
   DROP TABLE loc1;

VACUUM;

Lo script di creazione del database procede attraverso la generazione dello schema di quattro tabelle, a partire dai file testuali precedentemente compilati, quindi vengono create in sequenza due tabelle ulteriori 'loc1' e 'localita' per estendere con procedure di JOIN la tabella delle località con i nomi per esteso delle suddivisioni amministrative. Al termine le tabelle di appoggio vengono cancellate e si effettua la compattazione del database sqlite (VACUUM). Nel seguito allego il trace della creazione del database:

2015-03-10 12:26:35Luoghi
2015-03-10 12:26:35a1
2015-03-10 12:26:35a2
2015-03-10 12:26:35a3
2015-03-10 12:26:38loc1
2015-03-10 12:26:40localita
2015-03-10 12:26:42

Come si vede, l'operazione è abbastanza veloce (7 secondi in tutto). Il database finale ha una dimensione di 53,7 megabyte, piuttosto voluminoso. La sua lettura da python risulta comunque molto rapida, come vedremo nel prossimo post. Per chi è interessato, trasferisco i file nel repository github nella sottocartella astrometry/atlas