giovedì 22 settembre 2011

Condividere progetti QGis tra Linux e Windows (memo)

I motivi per cui si deve (purtroppo) mantenere una partizione windows in una postazione gis sono molti e, per quanto mi riguarda, sono imposti da altri...
Capita quindi di dover usare in Windows progetti di Quantum Gis creati in Linux; per far questo vado a correggere i percorsi che puntano ai dati con una semplice linea di comando (che scrivo qui per ricordarmela):

cat nomeProgettoCreatoInLinux.qgs | sed sed 's/\/media\/windows/c:/g' > nuovoNomeProgettoInWin.qgs

Ovvio che i dati devono stare in directories che entrambi i sistemi operativi possano leggere e scrivere; inoltre il percorso di mount della partizione windows (nell'esempio: /media/windows) e la lettera ('c:') indicante il disco possono cambiare da pc a pc.
In windows  sed e cat  possono essere forniti  dalla shell di grass (nel plugin di grass in QGis) o da MinGW  (Minimalist GNU for Windows) installato insieme a grass per Windows .

mercoledì 21 settembre 2011

Riproiettare un'intera cartella di shapefiles da linea di comando (memo)

Giusto per ricordare a me stesso e, forse, ad altri: riproiettare al volo da linea di comando un'intera cartella di shapefiles:

for f in inputFolder/*shp; do 
     shp=$(basename $f)
      ogr2ogr  -s_srs '+init=epsg:32632' -t_srs '+init=epsg:4326' \
      outputFolder/${shp} $f
done

dove:
inputFolder è la cartella con gli shapefiles da riproiettare;
outputFolder quella ove verranno memorizzati gli sphapefiles risultato.

Nell'esempio si riproietta da  WGS84 UTM32N a WGS84 (gradi decimali).

Ne ho già viste decine e decine di versioni, ma quando serve non  mi ricordo mai...

martedì 12 luglio 2011

Ancora su Garmin e GPSBabel

In una'installazione fresca di Ubuntu ho ancora avuto problemi di permessi nel leggere la porta USB alla quale viene collegato l'apparecchio garmin.

pcgis@pcgis-netbook:~$ lsb_release -a
No LSB modules are available.
Distributor ID: Ubuntu
Description:    Ubuntu 11.04
Release:        11.04
Codename:       natty

La soluzione è sempre quella indicata nel precedente post:

pcgis@pcgis-netbook:~$ cat /etc/udev/rules.d/71-garmin.rules 
# # http://wiki.openstreetmap.org/wiki/USB_Garmin_on_GNU/Linux
# # gps garmin con porta usb senza permessi
ATTRS{idVendor}=="091e", ATTRS{idProduct}=="0003", MODE="0660", GROUP="plugdev"

Notare che al nome del file si è dato un numero iniziale più alto (71) di quelli già presenti nella cartella.

Leggere il file README presente nella cartelle per provare a capire il perchè.

pcgis@pcgis-netbook:~$ dir /etc/udev/rules.d/
70-persistent-cd.rules  70-persistent-net.rules  71-garmin.rules  README

Il file gps_garmin risulta già in blacklist:

pcgis@pcgis-netbook:~$ cat /etc/modprobe.d/blacklist.conf | grep gps
blacklist garmin_gps



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