Grass, gis e dintorni

Condividere progetti QGis tra Linux e Windows (memo)

2011-09-22T21:00:00.000+02:00

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 .

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

2011-09-21T18:00:00.003+02:00

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

Ancora su Garmin e GPSBabel

2011-07-12T12:45:00.001+02:00

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

La mia prima volta... con python

2011-02-23T13:00:00.004+01:00

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

QLandkarte GT e lo screenshot capovolto

2010-05-19T23:31:00.001+02:00

Da qualche tempo uso l'ottimo programma QLandKarte GT per gestire il GPS Garmin.
Una funzionalità che può servire è quella della cattura della schermata del dispositivo.
Non so per quale strana ragione, ma da quando ho aggiornato il sistema ad Ubuntu Lucid lo screenshot ottenuto è ruotato e rispecchiato orizzontalmente (ma forse l'aggiornamento è del tutto ininfluente).
È comunque facile, con ImageMagick convertire gli screenshots ottenuti:
convert -rotate 180 -flop screen_bad.png screen_ok.png


Screenshot capovolto	...corretto con ImageMagick

Da Shapefile a KML

2010-01-27T19:02:00.002+01:00

Alcune brevi note sulla conversione di uno shapefile in un file kml per Google Earth.
Si voglia, ad esempio, visualizzare le aree protette del Piemonte su Google Earth, partendo da uno shapefile nel sistema UTM ED50 (codice epsg 23032).

Scaricare il file parchi_ed50.zip da questo indirizzo; scompattare il file con unzip:
unzip parchi_ed50.zip

Un solo comando con ogr2ogr è sufficiente per ottenere il file kml delle aree protette:
ogr2ogr -f "KML" -s_srs '+init=epsg:23032 +towgs84=-87.000,-98.000,-121.000' \
    -t_srs '+init=epsg:4326' parchi_1.kml parchi_ed50.shp

Una lettura del manuale in linea di ogr2ogr potrà aiutare la comprensione.
Si può scindere l'operazione in due parti; la prima è la conversione dal sistema (proiettato) UTM con datum ED50 al sistema geografico WGS84 (codice epsg 4326):
ogr2ogr -s_srs '+init=epsg:23032 +towgs84=-87.000,-98.000,-121.000' \
    -t_srs '+init=epsg:4326' parchi_wgs84.shp parchi_ed50.shp

da notare i parametri towgs84=-87.000,-98.000,-121.000 necessari per una corretta trasformazione (vedere questo precedente post).

La seconda operazione è la trasformazione dello shapefile in kml :
ogr2ogr -f "KML" parchi_2.kml parchi_wgs84.shp

ogr2ogr offre la possibilità di impostare interrogazioni SQL; ad esempio, se si vuole un file kml con il solo poligono riguardante l'area protetta denominata 'confluenza del Varaita' impartire il comando:
ogr2ogr -f "KML" -sql "SELECT * FROM parchi_wgs84 WHERE LOCALITA='confluenza del Varaita'" \
    -dsco NameField=LOCALITA conf_varaita.kml parchi_wgs84.shp

Le specifiche kml definiscono un campo denominato 'name' che contiene il nome dell'elemento (nel nostro caso un poligono); l'opzione -dsco imposta 'name' al contenuto del campo 'LOCALITA'.

From Grass, gis e dintorni

Un analogo risultato si può ottenere con il plugin OGR Layer converter di Quantum Gis :

From Grass, gis e dintorni

Qgis 1.4: interfaccia alle utilities GDAL.

2010-01-06T18:28:00.001+01:00

Qgis 1.4 presenta numerose nuove funzionalità. La prima ad essere apprezzata è stata l'interfaccia a molte utilities della libreria gdal delle quali non ricordo quasi mai la sintassi: in pratica un compositore che facilita enormemente l'immissione corretta di righe di comando.

Ecco un esempio:

From Grass, gis e dintorni

Il comando così composto può essere eseguito subito, oppure copiato ed incollato in uno script.

Qgis e la sorgente dati non valida.

2009-12-30T22:59:00.003+01:00

Raramente uso Qgis per l'elaborazione di raster memorizzati in mapset di GRASS; ho quindi, per lungo tempo, sopportato passivamente l'errore che ricevevo quando tentavo di caricare un raster di GRASS in Qgis: sorgente dati non valida.
L'ultima versione di Qgis ("Enceladus") è però molto ricca e l'uso della shell di GRASS nel plugin per Qgis pare uno strumento efficace; ho fatto allora qualche tentativo (positivo) per risolvere il problema dei raster.

Ecco gli appunti (non è... oro colato) di una soluzione poco ortodossa e temporanea:

Linux Ubuntu Karmic; nel file sources.list ho, tra gli altri, questi repositories:

deb http://ppa.launchpad.net/ubuntugis/ubuntugis-unstable/ubuntu karmic main
deb-src http://ppa.launchpad.net/ubuntugis/ubuntugis-unstable/ubuntu karmic main

Per il gioco delle dipendenze, nel sistema sono installati questi files:

/usr/lib/gdal15plugins/gdal_GRASS.so
/usr/lib/gdal15plugins/ogr_GRASS.so

/usr/lib/gdal16plugins/ogr_GRASS.so
/usr/lib/gdal16plugins/gdal_GRASS.so

Sospetto che i due pacchetti gdal vadano in conflitto; rinomino il pacchetto meno aggiornato e creo un link:

pcgis@laptop:/usr/lib/gdal15plugins$ sudo mv gdal_GRASS.so copia_di_gdal_GRASS.so
pcgis@laptop:/usr/lib/gdal15plugins$ sudo ln -s /usr/lib/gdal16plugins/gdal_GRASS.so gdal_GRASS.so

Qgis 1.4 ora carica senza problemi i raster di GRASS:

From Grass, gis e dintorni

GPS Garmin in Windows: i due (+1) indispensabili

2009-12-21T14:54:00.005+01:00

Dare gli strumenti indispensabili, facili da installare e liberi, per lavorare in Windows con apparecchi GPS in... un'ottica GIS...: brevi istruzioni per l'uso rivolte a nuovi utenti (entry level).

Map Source

Non essendo libero e facile da installare non dovrebbe essere menzionato in questo post; un accenno, però, occorre farlo. Considero Map Source utiile nel solo caso si debbano gestire mappe acquistate da Garmin: pare che questo software si preoccupi in modo quasi maniacale di tutto ciò che riguarda autenticazione, sblocco e registrazione delle mappe acquistate tralasciando di fornire all'utente un prodotto... utile e fruibile in tutte le occasioni di normale utilizzo (download e upload di waypoitns, tracce, rotte). Per quanto l'esperienza di chi scrive, l'utilità di Map Source inizia e finisce nel momento in cui, dopo una lunga e cervellotica fase di autenticazione nella quale vengono richiesti 3-4 codici di attivazione, le mappe acquistate sono sbloccate e possono essere caricate nel GPS.

Quantum Gis

Qgis può essere scaricato da qui: consiglio di scegliere il file denominato QGIS-1.3.0-3-No-GrassSetup.exe (o un analogo file con un numero più recente di versione). Scaricarlo ed installarlo accettando la licenza ed ogni altra preimpostazione (non è necessario scaricare i set di dati di esempio proposti durante l'installazione).

Una guida introduttiva a Qgis può essere trovata qui.

Il plugin GPS di Qgis sfrutta gpsbabel, un altro programma che deve essere presente (installato) nel sistema.

Gpsbabel

A giudizio di chi scrive, gpsbabel è uno tra i migliori software per lavorare con i GPS.
Scaricare da qui il file denominato gpsbabel-1.3.6.zip (o un analogo file con un numero più recente di versione). Per l'installazione:

creare una cartella sull'Hard disk (ad esempio: C:\Programmi\gpsbabel)
scomapttare il file gpsbabel-1.3.6.zip nella cartella appena creata; alla fine di questa operazione il contenuto di quest'ultima dovrebbe presentarsi così:

From gps
gpsbabel ha due eseguibili: GPSBabelGUI.exe che è un'interaccia grafica ; gpsbabel.exe che è il vero e proprio motore, eseguibile all'interno della finestra MS-DOS di Windows ;
gpsbabel , pur essendo installato, può essere - al momento - eseguito solo dalla cartella ove sono collocati i files (nell'esempio: C:\Programmi\gpsbabel). Per rendere disponibili le sue funzionalità all'esterno di questa cartella e per tutti gli altri programmi che lo utilizzano (ad es. qgis), occorre inserire il percorso degli eseguibili di gpsbabel nella variabile di sistema path; per far ciò (in Windows XP):
Impostazioni->pannello dai controllo->sistema->avanzate->variabili di ambiente: nella casella Variabili di sistema selezionare Path e modificarlo, aggiungendo alla linea il percorso ove sono stati installati gli eseguibili di gpsbabel (nell'esempio: C:\Programmi\gpsbabel; ricordarsi di separare i diversi percorsi con ";")
a questo punto gpsbabel dovrebbe funzionare...

attraverso la riga di comando della finestra MS-DOS di Windows :
attraverso l'interfaccia grafica di gpsbabel (doppio click su GPSBabelGUI.exe nella cartella C:\Programmi\gpsbabel oppure creare un link sul desktop)

From gps
attraverso il plugin GPS di qgis:

GRASS, gpsd, l'MTK e il bluetooth

2009-12-19T16:59:00.001+01:00

Nel precedente post si è visto come connetere l'MTK 747 Transystem a gpsd via bluetooth in Ubuntu 9.10.
Ecco un semplice memo di come utilizzare la connessione in real-time in GRASS.

Lanciare GRASS in un terminale ed entrare in una location con sistema di riferimento cartografico noto; dal terminale di GRASS lanciare gpsd e quindi gpspipe:

From Grass, gis e dintorni

Ora, con semplici comandi della shell è possibile ottenere in real-time diverse elaborazioni; per esempio con:
gpspipe -w -n 30 | grep 'GGA' | awk 'BEGIN{FS=" "}{print $5,$4}' |\
m.proj -i --q | awk '{sx+=$1;sy+=$2}END{ print " posizione media: ", sx/NR, sy/NR " samples: ", NR }'

si riproiettano (m.proj) le coordinate (Latitudine e Longitudine) inviate dal ricevitore nel sistema cartografico di riferimento della location corrente (nell'esempio: epsg: 32632), viene calcolata la posizione media del ricevitore collegato al pc e si stampa a video il risultato.

Complicando ulteriormente l'esempio, è possibile editare in real-time, catturando il punto in un vettoriale:

From Grass, gis e dintorni

Per le frasi NMEA è possibile consultare questo link.

MTK, GPSD via bluetooth

2009-12-19T16:14:00.001+01:00

Brevi note per connettere l'MTK 747 Transystem a  gpsd via bluetooth in Ubuntu 9.10.

Ho una chiavetta USB SURECOM per le connessioni bluetooth che viene vista dal sistema così:
lsusb
[...]
Bus 003 Device 003: ID 0a12:0001 Cambridge Silicon Radio, Ltd Bluetooth Dongle (HCI mode)
[...]

Accendere l'MTK 747 Transystem in posizione NAV ; provare se la perifirica Bluetooth viene riconosciuta:
hcitool scan
Scanning ...
    00:1C:88:00:CC:9C    iBT-GPS

La periferica viene riconosciuta ed associata al codice 00:1C:88:00:CC:9C.

Configuro la connessione:
paolo@laptop:~$ sudo gedit /etc/bluetooth/rfcomm.conf

inserisco questi settaggi:
rfcomm0 {
#    # Automatically bind the device at startup
     bind yes;
#    # Bluetooth address of the device
     device 00:1C:88:00:CC:9C;
#    # RFCOMM channel for the connection
    channel    1;
#    # Description of the connection
    comment "iBT-GPS";
}

Salvo e chiudo. N.B.: non è detto che debba per forza essere rfcomm0; se si hanno altri dispositivi bluetooth associare il dispositivo ad un'altra connessione (rfcomm1, rfcomm2, ecc.)

Quindi:
paolo@laptop:~$ rfcomm connect 0

oppure:
rfcomm connect 0 00:1C:88:00:CC:9C

se tutto è ok si avrà la seguente risposta:
Connected /dev/rfcomm0 to 00:1C:88:00:CC:9C on channel 1
Press CTRL-C for hangup

Aprire un altro terminale e digitare:
paolo@laptop:~$ gpsd /dev/rfcomm0
paolo@laptop:~$ xgps

From Grass, gis e dintorni

La disconnessione può essere effettuata dal gestore Bluetooth di Gnome:

From Grass, gis e dintorni

oppure premendo CONTROL-C nel primo terminale.
Buona cosa, inoltre, uccidere il processo gpsd
gpsd killall

GoogleEarth e files kml

2009-12-18T19:20:00.004+01:00

Google Earth in Linux: problema di visualizzazione non corretta dei file kml importati.

Dopo aver scaricato ed installato Google Earth in Ubuntu ho notato che i files kml importati non venivano posizionati correttamente; ho risolto il problema grazie a questo post; in pratica è sufficiente modificare il file googleearth nella cartella /home/nomeutente/google-earth (nel caso di installazione in locale, nella home utente) come segue:

LD_LIBRARY_PATH=.:${GOOGLEEARTH_DATA_PATH}:${LD_LIBRARY_PATH}
export LD_LIBRARY_PATH
export LC_ALL=us_US.UTF-8
export LANG=it_IT.UTF-8

Ecco un esempio, con una traccia denominata track_15_12_2009.gpx_trk.gpx:

step 1: conversione da gpx a kml con gpsbabel:
paolo@laptop:~$ gpsbabel -t -i gpx -f track_15_12_2009_trk.gpx -o kml -F track.kml

step 2: aprire Google Earth e caricare la traccia:

From Grass, gis e dintorni

Codici epsg per l'Italia Nord-Occidentale.

2009-12-09T22:05:00.005+01:00

Codici epsg: memo.

Ecco i codici di più frequente uso nell'Italia Nord-Occidentale (Piemonte, Liguria, Valle d'Aosta).

Codice	Descrizione	wkt
4806	Roma1940, in gradi.	wkt
26591	Sistema nazionale Gauss-Boaga, fuso Ovest. Datum Roma1940. Falsa origine delle ascisse: 1.500 km. Deprecato (usare 3003).	wkt
3003	Idem c.s.	wkt
23032	European 1950 datum (ED50). Grigliato UTM, fuso 32 Nord.	wkt
4326	WGS84 datum. Lat/Long.	wkt
32632	WGS84 datum. Grigliato UTM, fuso 32 Nord.	wkt

I parametri per le trasformazioni possono essere ottenuti con il comando di GRASS g.proj.

Garmin GPSmap 60 Cx e GPSBabel

2009-11-09T21:56:00.009+01:00

Usare il Garmin 60Cx con GPSBabel. Seguire le guide di OpenStreetMap e di GPSBabel.
In pratica, creare il file 51-garmin.rules:

paolo@laptop:~$ cat /etc/udev/rules.d/51-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"

Inserire il modulo garmin_gps in blacklist:

paolo@laptop:~$ cat /etc/modprobe.d/blacklist
#prevent garmin gops loaded so generic USB can be used instead
blacklist garmin_gps

A questo punto il ricevitore GPS viene riconosciuto come periferica usb:.
Per testare il tutto caricando una traccia dal pc al ricevitore:

cancellare ogni traccia dal ricevitore;
collegare il ricevitore;

caricare una traccia dal pc al ricevitore con:

gpsbabel -t -i gpx -f miatraccia.gpx -o garmin -F usb:

la traccia viene (dovrebbe!) essere caricata, ma il ricevitore
non la nomina;
nel menù tracce salvare la traccia dandole un nome

Point in polygon (terza parte)

2009-04-25T21:45:00.005+02:00

Joins spaziali e SQL

In questa sezione vengono affrontati gli aspetti principali del point in polygon maggiormente legati all'SQL.

Preparazione dei dati e prime verifiche

I dati necessari per l'esercitazione sono quelli ottenuti nell'esercitazione precedente; riassiumiamo:

la location simpleloc con il mapset pointInPoly;
un vettoriale poligonale polygons con attributi cat (la categoria dell'area), label (l'etichetta identificativa di ciascuna area) e GRASSRGB (il colore dell'area secondo la codifica RGB);
un vettoriale di punti mypt collegato, sul layer 1, alla tabella points_code con i campi cat e code (quest'ultimo codifica la posizione del punto rispetto agli elementi dei poligoni).

Merita porre l'attenzione sull'ultimo punto: un vettoriale è stato collegato ad una tabella con un nome totalmente diverso. Spesso si è indotti a pensare che ciò non sia possibile (a causa forse dell'abitudine ad usare formati di file quali lo shapefile) ma tutti i moduli di GRASS usati in questo esercizio funzionano benissimo anche in questo caso. Vediamo un esempio:

# # # vector and table with a different name
# # # mypt connected to points_table
table="points_code"
column="polyCat"
# # # if column not exists will be created
if ! checkColumn "$table" "$column" ; then
 echo "ALTER TABLE $table ADD COLUMN $column integer" | db.execute
fi
# # upload polygons cat to table
v.what.vect vect=mypt col=polyCat \
 qvector=polygons qcolumn=cat
# # # how many rows in table?
db.select sql="SELECT count(cat) FROM $table"
# # # how many points overlap polygons?
db.select sql="SELECT count(cat) FROM $table \
 WHERE $column NOT NULL"
# # # how many points don't overlap polygons?
db.select sql="SELECT count(cat) FROM $table \
 WHERE $column IS NULL"
# # # how many points overlaps boundaries, nodes, vertices, centroids, ecc.
db.select sql="SELECT code,count(cat) FROM $table \
 WHERE $column NOT NULL GROUP BY code"

Al solito vengono lasciati i commenti.

Una nota di spiegazione meritano le funzioni checkColumn e checkTable (usata dopo). Tutto il codice riportato è contenuto in una script che viene testato ogni qualvolta si apportano modifiche; il tentativo di creare una colonna (campo) già esistente fornisce un errore e lo stesso dicasi per il tentativo di creare una tabella già presente nel database. Ripetendo quindi l'esecuzione dello script senza un preventivo controllo sull'esistenza o meno di tabelle e colonne che si volgiono creare fornirebbe numerosi errori. Per evitare questo è possibile usare queste funzioni:

checkTable()
{
# true=0
# false=1
 table=$1
 db.tables -p | grep '\<'$table'\>' 2>/dev/null
 if [ $? -eq 0 ]; then
 exists=0
 else
 g.message -w message="Table $table doesn't exist"
 exists=1
 fi
 return $exists
}
checkColumn()
{
# true=0
# false=1
 local table=$1
 local column=$2
 if checkTable $table; then
 db.describe -ct $table | grep '\<'$column'\>' 2>/dev/null
 if [ $? -eq 0 ]; then
 exists=0
 else
 g.message -w message="Column $column doesn't exist"
 exists=1
 fi
 return $exists
 else
 exit 1
 fi
}

È chiaro che chi non ha dimestichezza con lo scripting di shell può lo stesso fare questi esercizi, senza dover per forza capire il meccanismo di funzionamento delle funzioni sopra riportate! Nell'esempio di prima basta omettere le righe che richiamano checkColumn:

# # # if column not exists will be created
 echo "ALTER TABLE points_code ADD COLUMN polyCat integer" | db.execute

...ricordandosi di cancellare la colonna polyCat con db.dropcol prima di ripetere l'esercizio, per evitare errori del tipo:

GRASS 6.5.svn (simpleloc):~/script_grass > echo "ALTER TABLE points_code ADD COLUMN polyCat integer" | db.execute
DBMI-SQLite driver error:
Error in sqlite3_prepare():
duplicate column name: polyCat

ERROR: Error while executing: 'ALTER TABLE points_code ADD COLUMN polyCat
 integer
 '

Riprendendo l'esercitazione, se proprio si vuole avere vettoriale e tabella con lo stesso nome si può procedere semplicemente così, copiando il vettoriale con un nuovo nome:

 # # vector and table with the same name
g.copy vect=mypt,mypt2 --overwrite 
db.describe -tc mypt2

oppure esportando e reimportando il vettoriale con v.out.ogr e v.in.ogr:

# # # vector and table with the same name
v.out.ogr in=mypt layer=1 type=point dsn=. \
 olayer=tmp_pippo format=ESRI_Shapefile -c
v.in.ogr dsn=. output=mypt2 layer=tmp_pippo --o

Si noti che con quest'ultimo metodo si ottiene una tabella leggermente diersa rispetto all'originaria (usare db.describe... viene inserito un campo cat_; il campo polyCat viene trasformato da integer a double precision).

Per esercitarsi appieno con SQL, si crei una tabella che memorizzi la desrizione della codifica dei punti così come definita nell'esercitazione precedente

# # # create a table for points descriptions
table="points_descr"
# # # if table exists will be cancelled
if checkTable $table ; then
 db.droptable -f $table
fi
echo "CREATE TABLE $table (cat integer, txtDescr varchar(25) )" | db.execute
# # insert data in table with loop and here document
while read data
 do
 fcat=$(echo $data | cut -d"|" -f1)
 ftxtDescr=$(echo $data | cut -d"|" -f2)
 echo "INSERT INTO $table (cat,txtDescr) \
 values ("$fcat",'"$ftxtDescr"')" | db.execute
done<<EOF
1|out (completly out!)
2|in (but not on boundaries nor centroids)
3|on_centroids
4|on_nodes (but not on vertices)
5|on_vertices (but not on nodes)
6|on_boundaries (but not on vertices nor nodes)
EOF
db.select $table

v.distance semplificato: v.what.vect

Si è visto nell'esercitazione precedente che v.distance permetto il trasferimento della categoria o degli attributi dai poligoni ai punti; v.what.vect non è altro che uno script per semplificare v.distance in questo genere di operazioni.

Per esempio si può aggiornare il campo polyCat dei punti con la categoria dei poligoni o trasferire un attributo (in questo caso il campo GRASSRGB):

# # # upload polygons cat to points attribute
v.what.vect vect=mypt2 col=polyCat qvector=polygons qcolumn=cat
# # # coloro i punti con lo stesso colore del poligono
# # # nel quale sono contenuti
if ! checkColumn "mypt2" "GRASSRGB" ; then
 v.db.addcol mypt2 col="GRASSRGB VARCHAR(11)"
fi
# # # upload polygons colors to points
v.what.vect vect=mypt2 col=GRASSRGB qvector=polygons qcolumn=GRASSRGB

Queste operazioni sono del tutto analoghe a quelle già viste in precedenza con v.distance.

Punti, poligoni e SQL

Si propone qui una semplice carrellata di esempi per dare un'idea generale sull'uso dell'SQL collegato con lo spatial join e il point in polygon; si rimanda agli ottimi tutorial presenti in rete per gli approfondimenti riguardanti l'SQL.

Verifiche sul numero dei punti

# # # how many points?
db.select sql="SELECT count(cat) FROM mypt2"
# # # how many points with polyCat not empty?
# # # how many points overlap polygons?
db.select sql="SELECT count(cat) FROM mypt2 \
 WHERE polyCat NOT NULL"
# # # as above-mentioned but group by code
db.select sql="SELECT code,count(cat) FROM mypt2 \
 WHERE polyCat NOT NULL GROUP BY code"
# # # how many points polyCat=NULL? (=points outside areas)
db.select sql="SELECT count(cat) FROM mypt2 \
 WHERE polyCat IS NULL"

Inner Join (equal join)

Nell'Inner join vengono estratti solo i record di entrambe le tabelle nei quali il campo chiave coincide.

Nel caso di mypt2 e polygons il campo chiave è rispettivamente polyCat e cat.

db.select sql="SELECT mypt2.cat AS pointCat, \
 polygons.cat, polygons.label AS polyLabel, \
 polygons.GRASSRGB AS polyColor \
 FROM mypt2 INNER JOIN polygons \
 ON mypt2.polyCat=polygons.cat"

L'espressione con i campi chiave mypt2.polyCat=polygons.cat viene valutata true solo quando i campi, provenienti dalle due tabelle, sono uguali.

È facile verificare che è indifferente porre la clausula WHERE mypt.polyCat NOT NULL perchè nell'inner join vengono estratti solo i record di entrambe le tabelle nei quali il campo chiave coincide:

db.select sql="SELECT mypt2.cat AS pointCat, \
 polygons.cat, polygons.label AS polyLabel, \
 polygons.GRASSRGB AS polyColor \
 FROM mypt2 INNER JOIN polygons \
 ON mypt2.polyCat=polygons.cat \
 WHERE mypt2.cat NOT NULL"

Si può usare anche un'altra notazione, che porta a risultati identici:

db.select sql="SELECT a.cat AS pointCat, \
 b.cat, b.label AS polyLabel, \
 b.GRASSRGB AS polyColor \
 FROM mypt2 a, polygons b \
 WHERE a.polyCat=b.cat"

È possibile ottenere i punti raggruppati per codice con la descrizione contenuta nella tabella points_descr:

db.select sql="SELECT mypt2.code AS code, count(mypt2.code) AS NumberOfPoints, \
 points_descr.txtDescr AS description \
 FROM mypt2 INNER JOIN points_descr \
 ON mypt2.code=points_descr.cat GROUP BY mypt2.code"

...o l'elenco completo dettagliato di tutti i punti

db.select sql="SELECT mypt2.cat, mypt2.code AS code, \
 points_descr.txtDescr AS description \
 FROM mypt2 INNER JOIN points_descr \
 ON mypt2.code=points_descr.cat ORDER BY code"

Un problema abbastanza comune è il calcolo del numero di punti contenuti in ciascun poligono:

db.select sql="SELECT COUNT(mypt2.cat) AS totalPoint, \
 polygons.label AS polyLabel FROM mypt2 \
 INNER JOIN polygons ON \
 mypt2.polyCat=polygons.cat GROUP BY polyCat"

Questo risutato può essere memorizzato in una tabella autonoma:

table="countPointsPerPoly"
# # # if table exists will be cancelled
if checkTable $table ; then
 db.droptable -f $table
fi
echo "CREATE TABLE $table AS SELECT count(mypt2.cat) AS numpoint, \
 polygons.cat AS refPoly FROM mypt2 \
 INNER JOIN polygons \
 ON mypt2.polyCat=polygons.cat \
 GROUP BY polygons.cat" | db.execute

...oppure memorizzato in un campo nella tabella attributi dei poligoni:

if ! checkColumn "polygons" "numpoints" ; then
 v.db.addcol polygons col="numpoints integer"
fi
#
for i in $(v.db.select polygons col=cat -c); do 
echo "UPDATE polygons SET numpoints= \
 (SELECT count(mypt2.cat) FROM mypt2 \
 WHERE mypt2.polyCat=$i \
 GROUP BY polyCat) WHERE polygons.cat=$i" | db.execute
done

L'aggiornamento del campo numpoints è realizzato attraverso un loop; per ogni cat della tabella dei poligoni viene eseguita la subquery rappresentata dalle istruzioni all'interno delle parentesi tonde che restituisce un elenco di valori. La successiva clausula WHERE polygons.cat=$i filtra i valori ottenuti dalla subquery imponendo la condizione polygons.cat=$i ottenendo un solo valore per ogni ciclo che viene inserito nel campo numpoints.

È possibile concatenare l'inner join; nell'esempio si ottiene, per ogni poligono il punto contenuto con il suo codice e la descrizione:

# # # double inner join
db.select sql="SELECT polygons.label AS polyLabel, mypt2.cat AS pointCat,\
 mypt2.code as pointCode, points_descr.txtDescr AS pointDescr \
 FROM polygons \
 INNER JOIN mypt2 ON polygons.cat=mypt2.polyCat \
 INNER JOIN points_descr ON mypt2.code=points_descr.cat"

Left join (left outer join)

Il risultato del left join include tutti i records della prima tabella e solo i records della seconda che rendono vera l'espressione con i campi chiave:

db.select sql="SELECT mypt2.cat AS pointCat, \
 polygons.cat, polygons.label AS polyLabel, \
 polygons.GRASSRGB AS polyColor \
 FROM mypt2 LEFT JOIN polygons \
 ON mypt2.polyCat=polygons.cat"

Attualmente (aprile 2009) il right join (right outer join) non è supportato dal driver Sqlite di Grass.

Analogamente a quanto fatto in precedenza, è possibile estrarre i dati dalle tre tabelle usando sia l'inner join che il left join:

db.select sql="select mypt2.cat AS pointCat, mypt2.code AS pointCode, \
 points_descr.txtDescr AS pointDescr, polygons.label AS polyLabel \
 FROM mypt2 \
 INNER JOIN points_descr ON mypt2.code=points_descr.cat \
 LEFT JOIN polygons ON mypt2.polyCat=polygons.cat"

Casi particolari

Si voglia per esempio estrarre, per ciascun poligono avente area maggiore ad un certo valore, il numero di punti in esso contenuti:

if ! checkColumn "polygons" "area"; then
 v.db.addcol polygons col="area double precision"
fi
v.to.db polygons opt=area col="area"
db.select sql="SELECT count(mypt2.polyCat) AS numPoint, \
 polygons.label as polyLabel, \
 polygons.area AS polyArea FROM mypt2 \
 INNER JOIN polygons ON mypt2.polyCat=polygons.cat \
 GROUP BY polygons.cat HAVING polyArea>=7 ORDER BY polyArea DESC"

Oppure è possibile calcolare il rapporto di densità tra punti contenuti ed area del poligono contenitore:

db.select sql="SELECT count(mypt2.polyCat) AS numPoint, \
 polygons.label as polyLabel, \
 polygons.area AS polyArea, round(numpoints/area,2) AS density FROM mypt2 \
 INNER JOIN polygons ON mypt2.polyCat=polygons.cat \
 GROUP BY polygons.cat ORDER BY density DESC"

Conclusioni

In queste tre esercitazioni si è visto a grandi linee come GRASS affronta il problema del point in polygon. I tools base sono v.select e v.distance; v.what.vect è uno script che semplifica l'uso di v.distance; l'SQL è stato usato per l'estrazione dei dati e per creare i joins. Merita ricordare che quando l'esigenza è di unire (merge) due tabelle, c'è a disposizione v.db.join: il suo funzionamento interno, al di lá dei diversi controlli , ricorda le operazioni usate per l'aggiornamento del campo numpoints attraverso un loop.

La conoscenza dell'SQL pare uno strumento indispensabile per l'utente GIS, in particolare quando si lavora con il modello vettoriale.

Tutto il codice riportato in questo esercizio può essere scaricato da qui

Tre anni di Gis Open Source

2009-04-16T22:51:00.002+02:00

Per una serie di motivi e coincidenze, nel marzo 2006 decisi di passare dal Gis proprietario in Windows al Gis Open Source in Linux.

Dopo tre anni circa mi sento di poter affermare la grande soddisfazione che questa scelta ha portato.

Poche cose mi mancano veramente del Gis-Windows:

l'ambiente layout di ArcView (dalla vers. 8.2 in poi); ps.map ha grandi potenzialità ma allestimenti raffinati richiedono un ambiente dedicato WYSIWYG o WYSIWYM (What You See Is What You Mean): in ogni caso allestendo una mappa di un certo livello qualcosa (e non solo il codice) si deve vedere in tempo reale man mano che si lavora e non solo dopo compilazione. gvSIG offre un ambiente del tipo ArcView vecchia maniera (vers. 3.2) ma come sarebbe bello avere tutto in GRASS!
l'editor vettoriale di Terranova ShArc: in 3 anni GRASS ha fatto passi da gigante dal punto di vista dell'editing vettoriale, gvSIG offre un ottimo ambiente ma non topologico; purtroppo funzionalità del tipo inseguimento linee del raster o rimozione interattiva dei dangles sono ancora lontane; nell'editing massivo di grandi quantità di elementi geometrici contano molto l'interfaccia, gli acceleratori da tastiera, la velocità di ridisegno, forma e dimensione del puntatore, ecc.
la gestione delle etichette (labels) di ArcView (dalla vers. 8.2); tra le nuove caratteristiche incluse in gvSIG (vers. 1.9) vi è un motore per il labeling avanzato; speriamo in bene...

Point in polygon (seconda parte)

2009-03-26T13:20:00.007+01:00

Classificazione di punti

Prima di affrontare gli aspetti del point in polygon maggiormente legati all'SQL è bene riassumere ed approfondire l'uso di v.select e v.distance.

L'esercizio che ci si propone di risolvere è il seguente: dati una copertura poligonale e un vettoriale di punti, attribuire i seguenti codici numerici:

1, ai punti completamente esterni alle aree;
2, ai punti completamente all'interno delle aree ma non sui boundaries né sui centroidi;
3, ai punti posti sui centroidi;
4, ai punti posti sui nodi ma non sui vertici dei boundaries;
5, ai punti posti sui vertici ma non sui nodi;
6, ai punti posti sui boundaries (ma non sui vertici né sui nodi).

In pratica si vuole una classificazione dei punti in funzione della posizione rispetto agli elementi dei poligoni (centroidi, boundaries, nodi, vertici).

La soluzione che si propone non è, probabilmente, ottimale ma... chi scrive non è riuscito a trovarne di migliori. Ben vengano i suggerimenti e le correzioni!

In generale la procedura prevede:

l'attribuzione di una categoria corrispondente al codice voluto agli elementi dei poligoni (centroidi, nodi, vertici, boundaries),
il trasferimento della categoria dagli elementi dei polgioni agli attributi dei punti mediante v.distance,
l'estrazione dei punti con il codice desiderato (v.extract) in un nuovo vettoriale;
la cancellazione dei punti nel vettoriale originario: vengono cancellati solo i punti di volta in volta già estratti con v.extract

Preparazione dei dati

Al solito, viene usata la location simpleloc. Entrare nel mapset pointInPoly creato in precedenza.

db.connect driver=sqlite database='$GISDBASE/$LOCATION_NAME/$MAPSET/pointInpoly.db'
g.copy vect=points,mypt,polyTopology,polygons --overwrite
# # # add table to mypt
v.db.addtable mypt
v.db.addcol mypt col="code integer"
g.copy vect=mypt,mypt_tmp --overwrite

Tutto il codice riportato in questo esercizio può essere scaricato da qui

Punti (solo) sui centroidi

Si attribuisce con v.category una categoria ai centroidi corrispondente al codice voluto (3); si trasferisce la categoria dai centroidi al campo code di mypt_tmp; si estraggono da mypt_tmp tutti i punti con il codice pari a 3 creando il vettoriale on_centroids; si cancellano in mypt_tmp tutti i punti con codice 3.

v.category in=polygons out=polycode \
 type=centroid layer=2 cat=3 step=0 --overwrite
v.distance from=mypt_tmp to=polycode \
 to_type=centroid to_layer=2 upload=cat \
 col=code dmax=0
v.extract in=mypt_tmp out=on_centroids where="code=3" --overwrite
v.edit map=mypt_tmp tool=delete where="code=3"
g.remove vect=polycode

Punti sui nodi (ma non sui vertici)

Si assume la definizione di nodo fornita in questa esercitazione. Occorre estrarre i nodi con v.to.points; in questo modo però, il punto finale di un boundaries può coincidere con il punto iniziale del boundaries successivo. Si rimuovo allora i nodi duplicati con v.clean. Si attribuisce poi la categoria 4 ai nodi estratti e, così come fatto per i centroidi, si trasferisce la categoria dai nodi al campo code di mypt_tmp. Si estraggono i punti in un nuovo file on_nodes e si rimuovono nel vettoriale my_tmp.

v.to.points -n in=polygons out=multiNode \
 type=boundary --overwrite
# # # remove duplicate nodes
v.clean in=multiNode out=nodeSimple \
 type=point tool=rmdupl --overwrite
# # # get a unique cat for each node
v.category nodeSimple out=node \
 layer=3 opt=add cat=4 step=0 --overwrite
# # # upload distance fron nearest node
v.distance from=mypt_tmp to=node to_layer=3 \
 to_type=point upload=cat \
 column=code dmax=0
v.extract in=mypt_tmp out=on_nodes where="code=4" --overwrite
v.edit map=mypt_tmp tool=delete where="code=4"
g.remove vect=multiNode,nodeSimple,node


Punti (solo) sui centroidi	Punti sui nodi

Punti sui vertici (ma non sui nodi)

Le operazioni sono del tutto simili a quelle già svolte per i nodi.

v.to.points -v in=polygons out=multiVertex \
 type=boundary --overwrite
# # # remove duplicate nodes
v.clean in=multiVertex out=vertexSimple \
 type=point tool=rmdupl --overwrite
# # # set unique cat for each vertx
v.category vertexSimple out=vertex \
 layer=3 opt=add cat=5 step=0 --overwrite
# # # upload distance fron nearest node
v.distance from=mypt_tmp to=vertex to_layer=3 \
 to_type=point upload=cat \
 column=code dmax=0
v.extract in=mypt_tmp out=on_vertices where="code=5" --overwrite
v.edit map=mypt_tmp tool=delete where="code=5"
g.remove vect=multiVertex,vertexSimple,vertex

Punti sui boundaries (ma non sui nodi né sui vertici)

A questo punto è doveroso spiegare il perchè si sono cancellati i punti da mypt_tmp ogni volta che si sono estratti in un nuovo file: se non lo si facesse sarebbe impossibile riuscire ad estrarre i punti posti sui boundaries (o meglio, lungo i boundaries) senza estrarre anche i punti posti sui vertici e i sui nodi dei boundaries medesimi (provare per credere). Come si è detto questa procedura è lungi dall'essere ottimale (sig!)...

v.category in=polygons type=boundary \
 out=arcs layer=2 opt=add cat=6 step=0 --overwrite
v.distance from=mypt_tmp to=arcs to_layer=2 \
 to_type=boundary to_layer=2 upload=cat \
 col=code dmax=0
v.extract in=mypt_tmp out=on_boundaries where="code=6" --overwrite
v.edit map=mypt_tmp tool=delete where="code=6"
g.remove vect=arcs


Punti (solo) sui vertici	Punti sui boundaries

Punti esterni alle aree

Con v.select è, a questo punto, banale estrarre i punti esterni alle aree (codice 1):

v.select -r ainput=mypt_tmp \
 binput=polygons \
 op=overlap out=out_area --overwrite
v.db.update out_area col=code value=1

Punti nelle aree (ma non sui centroidi né sui boundaries)

Si usa v.select con una certa sicurezza: avendo cancellato via via i punti da mypt_tmp non si dovrebbe (ancora una volta: la procedura qui proposta è lungi dall'essere ottimale) correre il rischio di estrarre punti non corretti.

v.select ainput=mypt_tmp \
 binput=polygons \
 op=overlap out=in_area --overwrite
v.db.update in_area col=code value=2


Punti esterni alle aree	Punti nelle aree (ma non sui centroidi, nodi,vertici, boundaries)

Creazione della tabella dei punti

Verrebbe da usare, a questo punto, v.patch:

v.patch -e in=in_area,out_area,on_boundaries,\
 on_nodes,on_vertices,on_centroids \
 out=points_code_vpatch --overwrite
v.db.select points_code_vpatch

ma ecco il risultato nella tabella attributi:

cat|code
12|2
[...]
18|2
36|2
38|1
[...]
47|1
81|6
[...]
87|6
110|4
[...]
152|5
171|3
[...]
174|3

Le categorie originarie di mypt vengono perse; se si vogliono preservare le categorie, occorre creare una tabella unica partendo dagli attributi dei diversi files di volta in volta estratti (on_centroids, on_nodes, ecc.) e connetterla (v.db.connect) sul layer 1 sovrascrivendo la connessione esistente mypt:

# # # test if table exists 
# # # from v.db.addtable script
table="points_code"
database=`db.connect -p | grep '^database' | cut -d':' -f2`
driver=`db.connect -p | grep '^driver' | cut -d':' -f2`
db.tables database="$database" driver="$driver" 2> /dev/null | grep "^$table$"
# # # if table exists it will be deleted
if [ $? -eq 0 ]; then
 db.droptable -f $table
fi

echo "CREATE TABLE $table AS \
 SELECT * FROM on_boundaries \
 UNION \
 SELECT * FROM on_nodes \
 UNION \
 SELECT * FROM on_vertices \
 UNION \
 SELECT * FROM in_area \
 UNION \
 SELECT * FROM out_area \
 UNION \
 SELECT * FROM on_centroids" | db.execute
v.db.connect map=mypt table=points_code key=cat layer=1 -o
v.db.select mypt layer=1

In questo modo la tabella points_code viene connessa al vettoriale mypt sul layer 1 (non è detto che il nome del vettoriale e il nome della tabella ad esso collegata debbano proprio essere uguali!).

La prima parte del codice serve solamente per testare se la tabella points_code esiste già; in caso positivo viene sovrascritta con la nuova creata dal codice SQL riportato.

Come si vede dal risultato di v.db.select la tabella sul layer 2 di mypt contiene due campi: la cat che corrisponde ovviamente a quella del vettoriale originario (points) e il code contenente la codifica basata sulla localizzazione del punto rispetto agli elementi della copertura poligonale.

Un'ultima semplice e grossolana verifica sul numero dei punti nei diversi files intermedi rispetto al vettoriale contenente tutti i punti (mypt):

or vectorpoint in "on_vertices" \
 "on_boundaries" \
 "on_nodes" \
 "on_centroids" \
 "out_area" \
 "in_area" \
 "mypt"; do
echo "$vectorpoint $(v.info -t $vectorpoint | grep points | cut -d"=" -f2)"
done

on_vertices 7
on_boundaries 5
on_nodes 5
on_centroids 4
out_area 10
in_area 8
mypt 39

La somma dei punti nei diversi vettoriali (on_centroids, on_nodes, ecc.) corrisponde al numero dei punti contenuti in mypt.

Point in polygon (prima parte)

2009-03-01T22:28:00.006+01:00

Obbiettivi

Trovare se una determinata area contiene/non contiene un dato punto è tra i più classici dei problemi GIS.

GRASS risolve il problema con alcuni tools:

v.select, v.distance e v.what.vect.

Questa esercitazione introduttiva affronta il problema del point in polygon, ed il relativo passaggio degli attributi dagli elementi poligonali ai punti contenuti (o a quelii più vicini, sovrapposti, ecc): lo spatial join.

Set di dati per l'esercitazione

Al solito l'esercitazione sarà condotta con la location simpleloc, la location creata apposta per fini didattici.

È necessario creare un nuovo mapset pointInPoly ed accedervi.

Per poter usufruire appieno dell'SQL (per la seconda parte di questa esercitazione) occorre connettere il mapset ad un database SQLite e copiare i due vettoriali dal mapset PERMANENT al mapset pointInPoly

db.connect driver=sqlite \
 database='$GISDBASE/$LOCATION_NAME/$MAPSET/pointInpoly.db'
g.copy vect=points,mypt,polyTopology,polygons --overwrite
v.db.addtable mypt

Le tabelle dei due vettoriali sono ora memorizzate nel database SQLite denominato pointInpoly.db.

Il vettoriale di punti usato è stato creato apposta per consentire di avere tutte le combinazioni possibili (punto in area, fuori area, su taluni centroidi, su alcuni nodi, su alcuni vertici, sui boundaries, molto vicino ai centroidi, molto vicino ai boundaries, ecc.).

Selezione semplice di punti (in,out,on...)

Dati due vettoriali, v.select consente di selezionare ed estrarre in un nuovo vettoriale gli elementi del primo (ainput) sovrapposti al secondo (binput).

Ecco in sequenza le diverse operazioni; vengono lasciati i commenti che servono da spiegazione:

# # # select all points into polygons
v.select ainput=mypt \
 binput=polygons btype=area \
 op=overlap out=inarea --overwrite

# # # select all points outside polygons
v.select -r ainput=mypt \
 binput=polygons \
 op=overlap out=outarea --overwrite

# # # select all points on centroid
v.select ainput=mypt \
 binput=polygons btype=centroid\
 op=overlap out=oncentroids --overwrite
 
# # #select all points on boundary
v.select ainput=mypt \
 binput=polygons btype=boundary \
 op=overlap out=onboundaries --overwrite

Al momento overlap è l'unico operatore di v.select ma nella roabmap di GRASS si legge anche di operatori di contenimento e adiacenza.

Le mappe che seguono illustrano i risultati ottenuti con v.select e sono ottenute in modo automatico tramite lo script che può essere scaricato da questo link


Punti all'interno delle aree	Punti sui boundaries

Punti sui centroidi	Punti non sovrapposti alle aree

Selezione di punti con passaggio di attributi

Il modulo v.distance spaventa un po': la sua sintassi è davvero complicata (almeno per me!) ma è ampiamente giustificata dall'elevato numero di analisi possibili. Il manuale riporta che v.distance trova, per ogni elemento del vettoriale from, il più vicino elemento nel vettoriale to: decisamente diminutivo! In realtà v.distance è un formidabile tool per realizzare joins spaziali.

Una nota generale: Si impara ad usare v.distance usandolo..., come per moltissimi altri comandi di GRASS che a prima vista spaventano per l'elevato numero di opzioni.

Upload della cat del più vicino poligono.

Il primo e più semplice dei problemi è il seguente: dato un vettoriale costituito da punti e un altro vettoriale costituito da poligoni, trovare per ciascun punto il più vicino poligono e caricare la categoria di quest'ultimo nella tabella degli attributi dei punti.

Si crea quindi una sorta di legame spaziale (=join) tra ciascun punto e il poligono più vicino.

# # # upload category of the NEAREST polygons
v.db.addcol mypt col="polyCatNearest integer"
# # # dmax=-1 (default, no limit): all points
v.distance from=mypt to=polygons \
 upload=cat col=polyCatNearest dmax=-1

Settando il parametro dmax al valore di -1 non si pongono limiti alla ricerca del più vicino poligono; pertanto, anche gli attributi dei punti esterni (non sovrapposti) alle aree vengono aggiornati con la cat del poligono più vicino; se si setta il valore pari a 0, saranno trovati solo i punti sovrapposti ai poligoni.

v.db.addcol mypt col="polyCatOverlap integer"
# # # dmax=0: only points that OVERLAP polygons
v.distance from=mypt to=polygons \
 upload=cat col=polyCatOverlap dmax=0

Si osservi ora la tabella degli attributi dei punti:

v.db.select mypt col=cat,polyCatNearest,polyCatOverlap

I punti 6 e 7 non sembrano essere stati presi in considerazione: la cat del poligono più vicino, per questi due punti, è nulla. I punti 6 e 7 sono in un'isola (area senza centroide e con categoria nulla): in questo caso il valore caricato da v.distance sugli attributi dei punti è nullo.

Se si vuole realizzare un upload di tutte le cat dei centroidi dei poligoni sugli attributi dei punti occorre settare il parametro to_type a centroid:

v.db.addcol mypt col="centroidCat integer"
v.distance from=mypt to=polygons to_type=centroid \
 upload=cat col=centroidCat dmax=-1 \
 out=distFromCentroid --overwrite

Si osservi che ora anche ai punti 6 e 7 è stata caricata la cat del centroide più vicino:

v.db.select mypt col=cat,polyCatNearest,polyCatOverlap,centroidCat \
 where="cat=6 OR cat=7"

Upload della cat del centroide posto ad una distanza minima data.

Sino ad ora si è usato v.distance per trovare gli oggetti (poligoni o centroidi) più vicini; è pure possibile aggiornare la tabella attributi dei punti con la cat dei centroidi trovati ad una distanza minima data:

v.db.addcol mypt col="centroidCatDmin integer"
v.distance from=mypt to=polygons to_type=centroid \
 upload=cat col=centroidCatDmin \
 dmin=3 out=distFormCentroidDmin --overwrite

In questo caso v.distance cerca i centroidi distanti almeno 3 unità di mappa da ciascun punto.


Distanza dai centroidi	Centroidi distanti almeno 3 unità di mappa dai punti

Upload dell'area dei poligoni sugli attributi dei punti sovrapposti.

Per creare una relazione spaziale tra punti e poligoni, è di per sè sufficiente caricare la cat dei poligoni sugli attributi dei punti con i metodi appena visti; come si vedrà in seguito, in tale maniera è possibile interrogare la tabella dei poligoni in funzione della posizione dei punti. Non sempre però è questo il risultato desiderato: si potrebbe ad esempio richiedere un passaggio diretto di taluni attributi dei poligono sulla tabella dei punti.

Quale esempio di upload di attributi si può calcolare l'area dei poligoni e caricarla nella tabella dei punti. In questo esempio verranno presi in considerazione solo i punti sovrapposti ai poligoni:

# # # upload all polygons (areas) attributes
# # # required fields (columns) in "from" vector
# # # must be present before run v.distance
v.db.addcol mypt col="polyArea double precision"
# # # add column for area polygons
v.db.addcol polygons col="area double precision"
# # # calculate polygons area size
v.to.db polygons opt=area col=area
# # # upload polygons attributes to OVERLAP points
v.distance from=mypt to=polygons \
 upload=to_attr col=polyArea \
 to_col=area dmax=0
v.db.select mypt

Upload della distanza fra punti e centroidi

v.distance è un calcolatore di distanze fra oggetti; ecco due esempi dell'upload delle distanze tra:

i punti e i centroidi dei poligoni;
i punti e i boundaries dei poligoni

# # # upload distance from the nearest centroid
v.db.addcol mypt col="centroidDist double precision"
v.distance from=mypt to=polygons to_type=centroid \
 upload=dist col=centroidDist \
 dmax=-1 
# # # upload distance from the nearest boundary (ALL points)
v.db.addcol mypt col="boundaryDist double precision"
v.distance from=mypt to=polygons to_type=boundary \
 upload=dist col=boundaryDist \
 dmax=-1 out=distFromBondary --overwrite

v.db.select mypt col=cat,centroidDist,boundaryDist

Come si sarà certamente notato è sempre possibile materializzare la distanza fra gli oggetti, settando il parametro opzionale out con il nome del file vettoriale desiderato.


Distanze tra punti e boundaries	Distanze tra punti e nodi

Casi particolari

Gli esercizi precedenti hanno consentito di apprendere le operazioni basilari con v.distance.

Ecco alcuni casi meno usuali.

Upload della distanza, della cat univoca e delle coordinate del più vicino nodo

Si assume la definizione di nodo fornita in questa esercitazione. Data quindi una copertura poligonale, estrarre i nodi rimuovendo ogni sovrapposizione, identificarli univocamente e applicare v.distance in modo da ottenere l'upload, sulla tabella attributi dei punti, di:

distanza fra nodo e più vicino punto;
cat univoca del nodo;
coordinate del nodo.

v.db.addcol mypt col="nodeDist double precision, \
 nodeCat integer, \
 node_x double precision, \
 node_y double precision"
# # # extract nodes from polygons
v.to.points -n in=polyTopology out=multiNode \
 type=boundary --overwrite
# # # remove duplicate nodes
v.clean in=multiNode out=nodeSimple \
 type=point tool=rmdupl --overwrite
# # # get a unique cat for each node
v.category nodeSimple out=node \
 layer=3 opt=add --overwrite
# # # upload distance fron nearest node
v.distance from=mypt to=node to_layer=3 \
 to_type=point upload=dist \
 col=nodeDist dmax=-1
# # # upload cat from nearest node
v.distance from=mypt to=node to_layer=3 \
 to_type=point upload=cat \
 col=nodeCat dmax=-1
# # # upload coordinates from nearest node
v.distance from=mypt to=node to_layer=3 \
 to_type=point upload=to_x,to_y \
 col=node_x,node_y \
 dmax=-1 out=distFromNode --overwrite
v.db.select mypt col=cat,nodeDist,nodeCat,node_x,node_y

Upload della cat univoca dei boundaries sugli attributi dei punti posti sui boundaries medesimi

Nel modello vettoriale implementato in GRASS i boundaries costituiscono i contorni delle aree e, normalmente, non sono identificati univocamente con una category.

v.db.addcol mypt col="boundCatOverlap"
# # # get unique cat (on layer 2)
v.category in=polyTopology type=boundary \
 out=arcs layer=2 opt=add --overwrite
v.distance from=mypt to=arcs to_layer=2 \
 to_type=boundary upload=cat \
 col=boundCatOverlap dmax=0
v.db.select mypt col=cat,boundCatOverlap

Upload della cat univoca dei boundaries sugli attributi dei punti posti entro due valori di distanza dati

È una veriante di quanto appena visto: i punti sui quali caricare la cat dei boundaries sono posti entro una distanza minima e massima dai boundaries medesimi.

v.db.addcol mypt col="boundCatZone double precision"
v.distance from=mypt to=arcs to_layer=2 \
 to_type=boundary upload=cat \
 col=boundCatZone out=distFromZone \
 dmin=0.2 dmax=2 --overwrite
v.db.select mypt col=cat,boundCatOverlap,boundCatZone

Punti compresi entro due valori di distanza (0.2 e 2) dai boundaries

Random Network Analisys in GRASS

2009-02-09T22:10:00.004+01:00

Questa è un'esercitazione introduttiva sulla Network Analisys che verrà condotta su un reticolo di linee e punti creati in modo casuale; su detto reticolo verranno applicati i diversi tools: creazione del grafo delle isodistanze dai punti, percorso minimo, problema del commesso viaggiatore, ecc.

Si farà uso del monitor per la rappresentazione dei risultati: (d.mon); alcune volte le operazioni usate potrebbero risultare di difficile comprensione a chi non avesse dimestichezza con la shell; si consiglia, in tal caso, di visualizzare i risultati nel Map Display nel modo usuale. L'intento è evidente: fare qualche semplice studio sui layouts per riutilizzarli in altri lavori.

Set di dati per l'esercitazione

Verrà usata, come al solito, location simpleloc; essendo la creazione del network completamente casuale, qualsiasi location può andare bene.

La location simpleloc è stata creata appositamente per fini didattici; per avere maggiori notizie vedere questo precedente blog.

La location simpleloc può essere scaricata da qui.

Creazione del network

Un network è formato da archi e nodi tra loro connessi a formare un reticolo.

Il procedimento per creare il reticolo in modo casuale è il seguente:

si estraggono a sorte un certo numero di punti (v.random);
da ogni punto si ricavano le coordinate (v.out.ascii) e si memorizzano in un file di testo;
si leggono le coordinate dal file e si crea un vettoriale di linee (v.in.ascii);
si pulisce il vettoriale di linee (v.clean)
si estraggono nuovamente a sorte un certo numero di punti che andranno a costituire i nodi del network;
si crea la connessione (v.net) tra nodi e reticolo di linee (archi).

È consigliabile scrivere tutti i comandi in un file di testo ed eseguirlo come script; sarà interessante ripetere l'esercitazione variando il numero di punti casuali ed alcuni altri parametri.

Si riportano nel seguito i comandi impartiti, lasciando i relativi commenti.

# # # set region to default values
g.region -d
# # # set number of random points
numpoint="120"
# # # create a random points
v.random out=rndpoints n=$numpoint --overwrite
# # # create a temporary file to store points coordiantes
tmpascii="$(g.tempfile pid=$$)"
# # # first row in standard mode: (L)ine number of points
echo "L $numpoint" > ${tmpascii}
# # # append coordinates to ascii file; formatting by awk
v.out.ascii in=rndpoints fs=" "| awk '{print $1 " " $2}' >> ${tmpascii}
# # # inport ascii file; bad => before cleaning
v.in.ascii -n in=${tmpascii} out=rndlines_bad format=standard --overwrite

Il numero di punti è settato a 120. È particolarmente utile osservare come si possa creare il reticolo di linee importando le coordinate dal file ASCII generato dai punti: quale applicazione pratica si pensi alla creazione di linee partendo da punti rilevati con strumenti GPS.

L'output di v.out.ascii viene incanalato (pipe, |) attraverso awk per formattare l'output delle coordinate.

Si passa ora alla pulizia del vettoriale:

# # # v.clean: order of tools is significant:
# # # try to change it (also you have to change order of
# # # thresholds) and see how vector shape changes
v.clean input=rndlines_bad out=rndlines_nocat type=line \
 tool=break,snap,rmline,rmsa,rmdupl \
 thresh=0,0.4,0,0,0 --overwrite
v.build.polylines in=rndlines_nocat out=polylines --overwrite
# # # add category to vector lines
v.category in=polylines out=rndlines opt=add --overwrite
# # # clean up vectors
g.remove vect=rndpoints,rndlines_bad,rndlines_nocat,polylines
rm $tmpascii

v.clean è tra i comandi più difficili di GRASS: provare ad invertire l'ordine dei tools per capire come agiscono...

Il passaggio attraverso la creazione di polylines (v.build.polylines) è necessario per rimuovere gli pseudonodi.

Si noti che il vettoriale di linee non ha categoria e pertanto occorre aggiungerla con v.category.

Il reticolo di linee è pronto; si passa ai nodi, randomizzando nuovamente:

g.region -d
numpoint="5"
v.random out=netpoints n=$numpoint --overwrite

Ottenuto il reticolato di linee (=archi) e i nodi, si crea il network connettendoli. Per prima cosa aggiungere una tabella a netpoints:

# # # add table to points
v.db.addtable map=netpoints col="cat integer, dist double precision"

Il campo dist viene popolato con il comando v.distance che calcola, per ciscun elemento puntuale, la distanza alla più vicina linea.

# # # upload distances from points to the nearest line
v.distance from=netpoints to=rndlines upload=dist col=dist

L'opzione connect di v.net connette i punti al reticolo di linee entro una certa soglia (threshold); per determinare tale soglia si sceglie il valore massimo che assume il campo dist:

# # # get the maximun distnce from points to lines
max_dist=$(db.select -c sql="SELECT max(dist) FROM netpoints" )

...e lo si aumenta di una certa quantità per avere la sicurezza che tutti i punti vengano connessi: viene usato awk per eseguire il semplice conto.

max_dist=$(echo "$max_dist" | awk '{print $1+0.01}')

L'ultimo passo è la creazione del network vero e proprio con la connessione dei punti algi archi:

# # # create network by connecting points and arcs at given threshold
v.net in=rndlines points=netpoints \
 out=network op=connect \
 alayer=1 thresh=$max_dist --overwrite

L'output di v.net è un vettoriale con due layer: sul layer 1 vi sono le linee (archi); sul layer 2 vi sono i punti (nodi).

Lo script sl_draw.sh può essere scaricato da questo indirizzo; deve essere ovviamente scomapttato, reso eseguibile (chmod ugoa+x sl_draw.sh) e posto in una directory in modo che l'interprete dei comandi possa trovarlo (controllare con echo $PATH).

./sl_draw.sh
d.vect network layer=1 display=shape color=grey
d.vect network layer=2 display=shape icon=basic/circle \
 size=8 color=red fcolor=blue
d.vect network layer=2 display=cat llayer=2 \
 lsize=10 lcolor=blue

Il random-network creato

Network Analysis

La network analisys si basa sul concetto di costo necessario per percorrere un determinato arco o per attraversare un determinato nodo. In questa sede si assumono i valori di default:

il costo di viaggio lungo ciascun arco è pari alla lunghezza dell'arco medesimo;
il costo di attraversamento di ciascun nodo è pari a 0

Isodistanze (isocosti)

Il comando v.net.iso peremtte di definire delle zone (costs) di uguale costo di viaggio a partire dai nodi inidcati (ccats):

v.net.iso in=network out=isonet alayer=1 nlayer=2 \
 ccats=1-100 costs=2,3,5 --overwrite
v.db.addtable isonet
./sl_draw.sh
for catcolor in "1 red 2" "2 green 2" \
 "3 blue 2" "4 grey 1"; do
 # # # parsing parameters in group with set
 # # # $1=category $2=color $3=width
 # # # see http://tldp.org/LDP/abs/html/loops1.html
 set -- $catcolor
 d.vect isonet width=$3 display=shape \
 color=$2 where="cat=$1" 
done
d.vect network layer=2 display=cat \
 icon=basic/circle size=6 llayer=2 \
 bgcolor=white bcolor=orange \
 fcolor=black lcolor=black xref=center yref=center

Il grafo dello isodistanze

Percorso minimo

Con v.net.path vengono calcolati i percorsi minimi dal nodo di categoria 1 a tutti gli altri nodi.

Con g.tempfile viene creato un file di testo temporaneo che memorizza, per ogni riga, i seguenti elementi:

la categoria univoca del percorso che si va a cercare;
la categoria del punto di partenza (nell'esempio si cercano le distanze dal punto con cat=1);
la categoria del punto di arrivo (nell'esempio tutti gli altri punti creati).

Il file temporaneo viene usato come imput per v.net.path.

 # # distance from point with cat=1 to all others points
points_net=`g.tempfile pid=$$`
# # cat_short_path from_node_cat to_node_cat
echo "1 1 2" > $points_net
echo "2 1 3" >> $points_net
echo "3 1 4" >> $points_net
echo "4 1 5" >> $points_net
v.net.path input=network output=shortpath \
 file=$points_net --overwrite
rm $points_net
# # # display 
./sl_draw.sh
d.erase
# # # split monitor 
d.frame -c frame=1 at=51,99,1,49
d.frame -c frame=2 at=51,99,51,99
d.frame -c frame=3 at=1,49,1,49
d.frame -c frame=4 at=1,49,51,99
# # # loop; catpath is the shortest category path
# # # and the frame name (1,2,3,4)
for catpath in 1 2 3 4; do
 d.frame frame=$catpath -s
 d.vect network display=shape color=grey
 d.vect map=shortpath display=shape,cat \
 width=2 color=red cats=$catpath \
 bgcolor=white bcolor=blue \
 lcolor=blue
 d.vect network layer=2 display=shape,cat \
 icon=basic/circle size=6 llayer=2 \
 bgcolor=white bcolor=orange \
 fcolor=black lcolor=black xref=center yref=center
 string=$(v.db.select shortpath col=cost -c where="cat=$catpath")
 d.text at=2,94 text="cost=$string" color=blue
done

Merita notare le istruzioni impartite per leggere, per ogni percorso, la distanza tra i due punti; v.net.path crea un vettoriale con attributi:

v.db.select shortpath col=cost -c where="cat=$catpath"

I percorsi minimi dal nodo 1 a tutti gli altri nodi.

Il problema del commesso viaggiatore e l'albero di Steiner

Connettere una serie di punti attraverso un percorso più corto possibile toccando (attraversando) ciascun punto una sola volta: questo è il problema del commesso viaggiatore (v.net.salesman).

L'albero di Steiner (v.net.steiner) realizza la connessione tra i nodi dati con il percorso a minor costo (più corto possibile).

Per ricordare: se voglio pianificare una gita o una vacanza scelgo l'algoritmo del travel salesman; se voglio posare un tubo per portare acqua a diversi irrigatori uso l'albero di Steiner.

La sintassi dei due comandi è abbastanza simile:

v.net.salesman in=network \
 output=salesmanpath alayer=1 \
 nlayer=2 ccats=1,2,3,4,5 --overwrite
# # #
v.net.steiner input=network \
 output=steiner alayer=1 \
 nlayer=2 tcats=1-5 --overwrite

Si passa alla visualizzazione dei risultati:

# # # dispay
./sl_draw.sh; d.erase
# # # split monitor
d.frame -c frame=1 at=51,99,1,49
d.frame -c frame=2 at=51,99,51,99
d.frame -c frame=3 at=1,49,1,49
d.frame -c frame=4 at=1,49,51,99
# # # network in frame 1
d.frame -s frame=1
d.vect network display=shape color=grey
d.vect network layer=2 display=shape \
 icon=basic/circle size=6 fcolor=red
d.vect network layer=2 display=cat \
 llayer=2 lcolor=red lsize=10 \
 xref=left yref=bottom
d.text at=2,2 text="Network" color=black 
# # # salesmanpath in frame 2
d.frame -s frame=2
d.vect salesmanpath display=shape \
 color=magenta width=2
d.text at=2,2 text="Salesman path" color=black 
d.vect network layer=2 display=cat \
 bgcolor=white bcolor=orange llayer=2 \
 fcolor=black lcolor=black xref=center yref=center
# # # steiner tree in frame 3
d.frame -s frame=3
d.vect steiner display=shape \
 color=violet width=2
d.text at=2,2 text="Steiner tree" color=black 
d.vect network layer=2 display=cat \
 bgcolor=white bcolor=orange llayer=2 \
 fcolor=black lcolor=black xref=center yref=center
# # # which is the shortest ?
# # # steiner or salesman ?
# # # text in frame 4
d.frame -s frame=4
top=90
for vector in "salesmanpath" "steiner"; do
 # # # add table to vector
 v.db.addtable $vector col="length double precision"
 # # # upload vector length
 v.to.db $vector opt=length col=length
 # # # summarize length 
 length=$(db.select -c sql="SELECT sum(length) FROM $vector")
 # # # print 
 d.text.freetype -s text="$vector = $length" \
 size=10 at=2,$top color=black font=Vera
 # # # decrease text position in frame
 top=$(($top-8))
done

Il salesmann path e lo Steiner tree.

L'uso dei frames rende possibile la visualizzazione delle differenze; la stessa area geografica viene ripetuta in tre frames, ogni volta rappresentando qualcosa di diverso ma mantenendo costanti alcuni riferimenti (nel nostro caso: i nodi del network): non si vuole far vedere dove esattamente l'albero di Steiner e il Salesman path passano sul network; si vuole semplicemente rtrasmettere il messaggio che i due vettoriali realizzano connessioni con proprietà diverse.

Allocazione delle risorse

v.net.alloc crea una ripartizione del network tra i diversi nodi forniti in input:

v.net.alloc input=network output=alloc \
 alayer=1 nlayer=2 ccats=1,2,3,4,5 --overwrite

Il vettoriale di output non ha tabella attributi;se ne aggiunge una calcolando la lunghezza complessiva dei diversi subnets (uno per ciascuna riga della tabella) e quella totale di tutto il network:

# # # add table 
v.db.addtable alloc layer=1 col="length double precision"
# # # upload lenght
v.to.db alloc layer=1 col=length opt=length
# # # summarize total netwrok lenght
tot_length=$(db.select -c \
 sql="SELECT sum(length) FROM alloc")

Si procede ora alla visualizzazione:

# # # colorize
v.colors alloc range=1,5 color=bcyr column=cat

./sl_draw.sh; d.erase;
# # # split monitor into 4 frames
d.frame -c at=1,99,1,65 frame=1
d.frame -c at=67,99,66,99 frame=2
d.frame -c at=34,66,66,99 frame=3
d.frame -c at=1,33,66,99 frame=4

# # # overview in frame 1
d.frame -s frame=1
d.vect -a alloc
d.vect network layer=2 display=cat \
 bgcolor=white bcolor=orange llayer=2 \
 fcolor=black lcolor=black xref=center yref=center
 d.text.freetype -s at=2,2 text="Total lenght: ${tot_length}" \
 size=8 color=black
# # # first number = frame
# # # second number = subnet
for subnet in "2 2" "3 3" "4 5"; do
 set -- $subnet
 d.frame -s frame=$1
 d.vect alloc color=grey
 d.vect -a alloc where="cat=$2" width=2
 subnet_length=$(db.select -c \
 sql="SELECT sum(length) FROM alloc WHERE cat=$2")
 d.text.freetype -s at=2,2 text="subnet n. ${2}: ${subnet_length}" \
 size=8 color=black 
done

Il grafo dell'allocazione delle risorse.

Conclusioni

Non si è affrontato la Network Analisys in modo approfondito: giusto una presentazione dei principali tools che GRASS mette a disposzione, senza considerare i versi di percorrenza degli archi ed i relativi costi.

Il codice di questa esercitazione si presta a molteplici varianti: creazione di vettoriali random, effetto dei tools di v.clean nel conferire un particolare aspetto finale al vettoriale, allestimento dei layouts con i soli comandi d.mon e d.vect.

Gdal e Grass from source in Ubuntu Intrepid Ibex - 2

2009-02-03T22:01:00.004+01:00

Faccio un seguito a questo precedente blog: installare GDAL e GRASS da sorgente è un tormentone, almeno per me.

Ecco gli appunti che uso, sotto forma di script; spero possano essere utili a qualcuno.

Non lanciare alla cieca lo script: quasi sicuramente non funzionarà; leggerlo per capire le diverse operazioni ed eventualmente modificarlo. Gran parte del contenuto è tratto dai link riportati.

#! /usr/bin/env bash
# # # january 2009
# # # this isn't a script! It's a reminder only
# # # SO THINK TWICE BEFORE RUNNING (AT YOUR OWN RISK)

# # # you should read these links before try:
# # # http://grass.osgeo.org/wiki/Compile_and_Install
# # # http://grass.osgeo.org/grass64/source/REQUIREMENTS.html
# # # http://casoilresource.lawr.ucdavis.edu/drupal/node/123
# # # http://mpa.itc.it/markus/useful/conf_install_gdal_ogr_grass_plugin.sh
# # # http://trac.osgeo.org/gdal/wiki/GRASS
# # # http://trac.osgeo.org/gdal/wiki/BuildingOnUnix

PROGRAM=`basename $0`
CURRENTDIR="`pwd`"

# # # change 

# # # gdal source directory
DIR_SRC_GDAL="/usr/local/src/gdal"
# # # gdal install directory
DIR_INST_GDAL="/usr/local"
# # # grass source directory
DIR_SRC_GRASS="/usr/local/src/grass6_devel"
# # # grass install directory
DIR_INST_GRASS="/usr/local/grass-6.5.svn"
# # # plugin grass-gdal source directory
DIR_SRC_PLUGIN="/usr/local/src/gdal-grass-1.4.3"
# # # plugin grass-gdal install directory
DIR_INST_PLUGIN="/usr/local/lib/gdalplugins"


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
ask_yes_no()
{
 local yn=
 while [ "$yn" = "" ]; do
 echo "$1"
 read yn
 case $yn in
 y|Y) yn=0 ;;
 n|N) yn=1 ;;
 e|E)
# cd ${CURRENTDIR}
 exit 1 ;;
 *) yn=
 echo "Invalid response - please answer y or n"
 ;;
 esac
 done
 return $yn
}

if ask_yes_no "Update/upgrade system [Y/N/(E)xit] ?"; then
 sudo apt-get update
 sudo apt-get upgrade
fi

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
# # #
# # # GDAL - Geospatial Data Abstraction Library
# # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 

# # # cleaning
if ask_yes_no "Clean GDAL source [Y/N/(E)xit] ?"; then
 cd ${DIR_SRC_GDAL}
 make distclean
fi

# # # upgrade from svn
if ask_yes_no "Upgrade GDAL source from svn [Y/N/(E)xit] ?"; then
# # #
# # # FIRST INSTALL:
# # # sudo mkdir /usr/local/src
# # # sudo chown my_user_name:my_user_group /usr/local/src/
# # # cd /usr/local/src
# # # current development version:
# # # svn checkout https://svn.osgeo.org/gdal/trunk/gdal gdal 
# # # current stable version:
# # # svn checkout https://svn.osgeo.org/gdal/branches/1.5/gdal gdal
# # #
# # # AFTER FIRST INSTALL:
 cd ${DIR_SRC_GDAL}
 svn up
 cd ${CURRENTDIR}
fi

# # # configure ***WITHOUT GRASS***
if ask_yes_no "Configure GDAL [Y/N/(E)xit] ?"; then
 cd ${DIR_SRC_GDAL}
 ./configure \
 --with-pg=/usr/bin/pg_config \
 --prefix=${DIR_INST_GDAL} \
 --with-python \
 --with-ogdi \
 --with-sqlite \
 --with-mysql=/usr/bin/mysql_config \
 --with-libtiff=internal \
 --without-grass \
 -with-hide-internal-symbols
 cd ${CURRENTDIR}
fi

# # # make
if ask_yes_no "Make GDAL [Y/N/(E)xit] ?"; then
 cd ${DIR_SRC_GDAL}
 make
 cd ${CURRENTDIR}
fi

# # # make install
if ask_yes_no "Install GDAL [Y/N/(E)xit] ?"; then
 cd ${DIR_SRC_GDAL}
 sudo make install
 cd ${CURRENTDIR}
fi

# # # Is directory, where binary have been installed, in PATH variable?
if ask_yes_no "Check if PATH variable contains path to GDAL binaries [Y/N/(E)xit] ?"; then
 echo "GDAL binary are here ...${DIR_INST_GDAL}/bin"
 echo ""
 cd "${DIR_INST_GDAL}/bin"
 ls -l
 echo "PATH variable:"
 echo "$PATH"
fi

if ask_yes_no "You can change PATH variable now...[Y/N/(E)xit] ?"; then
 echo "There many ways to do this; perhaps you can modify .bashrc..."
 gedit $HOME/.bashrc
fi

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
# # #
# # # GRASS GIS (Geographic Resources Analysis Support System) 
# # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 

# # # cleaning
if ask_yes_no "Clean GRASS source[Y/N/(E)xit] ?"; then
 cd ${DIR_SRC_GRASS}
 make distclean
 cd ${CURRENTDIR}
fi

# # # upgrade from svn
if ask_yes_no "Upgrade GRASS from svn [Y/N/(E)xit] ?"; then
# # #
# # # FIRST INSTALL: (6.5.0 development branch)
# # # Note: development here is discouraged and
# # # should take place in GRASS 7.
# # # 
# # # svn checkout https://svn.osgeo.org/grass/grass/branches/develbranch_6 grass6_devel
# # # To switch from SVN 'trunk' (now GRASS 7 development)
# # # to 'develbranch_6' (now GRASS 6.4 development) use
# # # cd /path/to/your/local/copy/trunk
# # # svn switch https://svn.osgeo.org/grass/grass/branches/develbranch_6 .
# # #
# # # AFTER FIRST INSTALL:
 cd ${DIR_SRC_GRASS}
 svn up
 cd ${CURRENTDIR}
fi

# # # configure GRASS
if ask_yes_no "Configure GRASS [Y/N/(E)xit] ?"; then
 cd ${DIR_SRC_GRASS}
# # # --with-gdal="${DIR_INST_GDAL}/bin/gdal-config" \
 # # # CFLAGS="-O2 -march=pentium-m -Wall" LDFLAGS="-s"
 # # # CFLAGS="-ggdb -Wall -Werror-implicit-function-declaration"
 CFLAGS="-g -Wall" LDFLAGS="-s" ./configure \
 --with-gdal="${DIR_INST_GDAL}/bin/gdal-config" \
 --with-tcltk-includes=/usr/include/tcl8.5 \
 --with-postgres=yes \
 --with-postgres-includes=/usr/include/postgresql \
 --with-sqlite \
 --with-cxx \
 --with-blas \
 --with-lapack \
 --with-cairo \
 --with-fftw \
 --with-freetype=yes \
 --with-freetype-includes=/usr/include/freetype2 \
 --with-readline \
 --with-opengl-includes=/usr/include/GL \
 --with-mysql \
 --with-mysql-includes=/usr/include/mysql \
 --with-python=/usr/bin/python2.5-config \
 --with-wxwidgets=/usr/bin/wx-config \
 --enable-largefile=yes \
 --with-proj-includes=/usr/include \
 --with-proj-share=/usr/share/proj \
 --with-proj-libs=/usr/lib
 cd ${CURRENTDIR}
fi

# # # make
if ask_yes_no "Make GRASS [Y/N/(E)xit] ?"; then
 cd ${DIR_SRC_GRASS}
 make
 cd ${CURRENTDIR}
fi

# # # make vdigit nviz
if ask_yes_no "Make vdigit and nviz [Y/N/(E)xit] ?"; then
# # # symbolic link (overwrite if exists)
 sudo ln -s -f /usr/lib/python2.5/site-packages/wx-2.8-gtk2-unicode/wx/_gdi_.so \
 /usr/local/lib/libgdi.so 
 cd "${DIR_SRC_GRASS}/gui/wxpython/vdigit"
 make
 cd "${DIR_SRC_GRASS}/gui/wxpython/nviz"
 make
 cd ${CURRENTDIR}
fi

# # # make install
if ask_yes_no "Install GRASS [Y/N/(E)xit] ?"; then
 cd ${DIR_SRC_GRASS}
 sudo make install
 cd ${CURRENTDIR}
fi

# # # check/set/modify shared library
if ask_yes_no "Check ld.so.conf [Y/N/(E)xit] ?"; then
 echo ""
 echo "REMEMBER: ${DIR_INST_GRASS}/lib must be shared!"
 echo ""
 cat /etc/ld.so.conf
fi
# # # /etc/ld.so.conf
if ask_yes_no "Modify ld.so.conf [Y/N/(E)xit] ?"; then
 sudo gedit /etc/ld.so.conf
 sudo ldconfig
fi

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
# # #
# # # PLUGIN GDAL - GRASS
# # # http://download.osgeo.org/gdal/gdal-grass-1.4.3.tar.gz
# # # 
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 

if ask_yes_no "Compile gdal-grass plugin [Y/N/(E)xit] ?"; then
# # # FIRST INSTALL
# # # cd /usr/local/src
# # # tar xfvz gdal-grass-1.4.3.tar.gz
# # # AFTER FIRST INSTALL:
 cd ${DIR_INST_GDAL}/lib
 # # # link GRASS libs to the lib/ folder
 sudo ln -sf ${DIR_INST_GRASS}/lib/*.so .
 cd ${DIR_SRC_PLUGIN}
 make distclean

 ./configure --with-grass=${DIR_INST_GRASS} \
 --with-gdal=${DIR_INST_GDAL}/bin/gdal-config \
 --with-autoload=${DIR_INST_PLUGIN} \
 --with-ld-shared="g++ -shared"
 
 make
# # # sudo mkdir /usr/local/lib/gdalplugins # # # first install
 sudo cp *.so ${DIR_INST_PLUGIN}
fi

Full planar topology in GRASS (Seconda parte)

2009-01-06T23:26:00.006+01:00

Node Topology

Occorre ottenere un vettoriale di nodi; per fare ciò si può usare il tool v.to.points:

v.to.points -n input=arcs out=multiNode \
 llayer=2 type=boundary --overwrite

I nodi vengono estratti dal layer 2 di arcs; volendo mantenere una relazione tra arco e nodo si usa l'opzione llayer che indica quale layer utilizzare (il ) per memorizzare la categoria dell'arco (=boundary) negli attributi dei nodi. v.to.points genera due layer: il primo con gli archi originali, il secondo con i nodi:

 v.db.select multiNode layer=2
cat|lcat|along
1|1|0
2|1|6.24264068711928
3|2|0
4|2|3
5|3|0
6|3|6.70820393249937
[...]
23|12|0
24|12|2.41421356237309
25|13|0
26|13|7.07106781186548

Il campo lcat memorizza la categoria dell'arco (layer 2 di arc) sul quale veine creato il nodo.

Sul layer 2 di multiNode si sono ottenuti i punti (=nodi) iniziali e finali di ciascun arco: va da sé che i nodi (finale di un arco ed iniziale dell'arco successivo) sono sovrapposti uno sull'altro.

Questo esercizio è una simulazione della Full Planar Topology: la tabella Node Topology necessita che i nodi siano unici ed univocamente individuati. Si procede allora a una pulizia dei nodi sovrapposti con v.clean tool=rmdupl:

v.clean in=multiNode out=nodeSimple \
 type=point tool=rmdupl --overwrite

Al vettoriale risultato si attribuisce un valore di categoria univoco:

v.category nodeSimple out=node \
layer=3 opt=add --overwrite

v.clean tool=rmdupl ha rimosso i nodi sovrapposti da multiNode; il meccanismo di memorizzazione dei vettori implementato in GRASS rende però possibile il recupero dei valori di categoria per ciascun layer grazie al mantenimento di una relazione uno a molti: interrogando dal layer 3 (parte uno) di nodeSimpel il layer 2 (parte molti), è possibile conoscere quanti e quali valori di categoria sono presenti in quest'ultimo layer. Non solo: con lo stesso meccanismo si possono popolare i campi del layer 3 con altri attributi del layer 2.

Per prima cosa occorre aggiungere alcuni campi al layer 3:

v.db.addtable node layer=3 \
 columns="cat integer, countnodes integer, \
 countarcs integer, nodeslist char(12), \
 arcslist char(12), arcslistnames char(20)"

Il tool che consente di realizzare quanto sopra descritto è v.to.db opt=query : per ogni geometria esegue la funzione contenuta in qcolumn interrogando il layer qlayer e popolando il campo definito in column. I commenti lasciati fra le righe di comando servono da spiegazione:

# # # numero di nodi sovrapposti nel layer 2
v.to.db node option=query layer=3 \
 qlayer=2 column=countnodes qcolumn="count(cat)"
# # # numero di archi collegati ai nodi 
v.to.db node option=query layer=3 \
 qlayer=2 column=countarcs qcolumn="count(lcat)"
# # # elenco delle categorie dei nodi sovrapposti sul layer 2
v.to.db node option=query layer=3 \
 qlayer=2 column=nodeslist qcolumn="group_concat(cat)"
# # # elenco delle categorie delgi archi collegati
v.to.db node option=query layer=3 \
 qlayer=2 column=arcslist qcolumn="group_concat(lcat)"
# # # elenco dei nomi degli archi collegati 
v.to.db node option=query layer=3 \
 qlayer=2 column=arcslistnames qcolumn="group_concat(arcLabel)"

Il risultato è il seguente:

v.db.select node layer=3
cat|countnodes|countarcs|nodeslist|arcslist|arcslistnames
1|2|2|1,2|1,1|a1,a1
2|3|3|6,9,17|3,5,9|a3,a5,a9
3|3|3|15,18,19|8,9,10|a8,a9,a10
4|3|3|4,5,20|2,3,10|a2,a3,a10
5|4|4|7,10,11,22|4,5,6,11|a4,a5,a6,a11
6|4|4|8,12,14,23|4,6,7,12|a4,a6,a7,a12
7|3|3|21,24,25|11,12,13|a11,a12,a13
8|4|4|3,13,16,26|2,7,8,13|a2,a7,a8,a13

Si può verificare la situazione delle categorie nei diversi layers:

v.category node opt=report
Layer/table: 1/node_1
type count min max
point 25 1 13
line 0 0 0
boundary 0 0 0
centroid 0 0 0
area 0 0 0
all 25 1 13
Layer/table: 2/node_2
type count min max
point 26 1 26
line 0 0 0
boundary 0 0 0
centroid 0 0 0
area 0 0 0
all 26 1 26
Layer/table: 3/node_3
type count min max
point 8 1 8
line 0 0 0
boundary 0 0 0
centroid 0 0 0
area 0 0 0
all 8 1 8

La tabella Node Topology è semplicemente:

 v.db.select node layer=3 col=cat,arcslistnames
cat|arcslistnames
1|a1,a1
2|a3,a5,a9
3|a8,a9,a10
4|a2,a3,a10
5|a4,a5,a6,a11
6|a4,a6,a7,a12
7|a11,a12,a13
8|a2,a7,a8,a13

L'arco a1 si chiude su se stesso nel nodo di categoria 1. La tabella prodotta presenta in corrispondenza del nodo 1 la ripetizione del nome dell'arco: ciò costituisce una differenza fra quanto, sulla base del modello Full Planar Topology, dovrebbe essere.

Si crea ora una visualizzazione del vettoriale node :

./sl_draw.sh
d.font.freetype font=timesbd
d.vect arcs layer=2 display=shape,attr attrcol=arcLabel \
 type=boundary llayer=2 \
 xref=center yref=center color=grey \
 lsize=10 bgcolor=white bcolor=black
d.vect node layer=3 display=shape,cat llayer=3 \
 size=5 icon=basic/triangle lcolor=aqua \
 xref=center yref=center \
 lsize=10 bgcolor=white bcolor=aqua

Relazioni tra i vettoriali node e arcs.

Polygon Topology

La terza ed ultima tabella Polygon Topology viene ottenuta combinando gli attributi sui due layers del vettoriale arcs.

Nel campo arcList sarà memorizzato l'elenco degli archi che delimitano ciascuna area:

 v.db.addcol arcs col="arcList char(15)"

Si provi dapprima questa istruzione SQL:

echo "select a.cat as cat,a.label as label, \
 group_concat(b.arcLabel) as arcList \
 from arcs a, arcs_2 b \
 WHERE (a.cat==b.catPolyLeft) OR \
 (a.cat==b.catPolyRight) \
 GROUP BY a.cat" | db.select

Ancora una volta si fa uso della funzione di aggregazione comando di aggregazione group_concat supportato nelle recenti versioni di SQLite; inoltre le diverse tabelle sono indicate da semplici lettere.

Si crea una tabella temporanea per memorizzare il risultato dell'istruzione SELECT e subito dopo la si interroga per aggiornare il valore del campo arcList:

echo "CREATE TEMPORARY TABLE foo AS \
 select arcs.cat as cat, group_concat(arcs_2.arcLabel) as List \
 from arcs, arcs_2 \
 WHERE (arcs.cat==arcs_2.catPolyLeft) OR \
 (arcs.cat==arcs_2.catPolyRight) \
 GROUP BY arcs.cat;
 UPDATE arcs SET arcList=(select List FROM \
 foo WHERE arcs.cat=foo.cat)" | db.execute

Notare che:

si è usata la notazione consueta: nomeTabella.nomeCampo anzichè le lettere per indicare le tabelle;
si sono concatenate le due istruzioni SQL con il ';';
viene creata una tabella temporanea; la sintassi completa del comando CREATE TABLE è riportato manuale di SQLite;
la tabella temporanea verrà distrutta appena la connessione al database cesserà.

Conclusioni

Arrivati al fondo di questa tortuosa esercitazione merita puntualizzare alcune cose:

l'esercitazione condotta è stata una simulazione: GRASS implementa il modello topologico in modo ben più efficiente!
la procedura si presta a numerose varianti: si provi ad esempio ad usare l'opzione -v (anziché -n) del tool v.to.points; la via seguita è quella che si è ritenuto più didattica...
mentre le tabelle Node Topology e Arc Topology sono state ottenute interrogando le geometrie, la terza tabella, Polygon Topology è stata ottenuta ricombinando altre tabelle; si provi ad ottenere Polygon Topology direttamente dalle geometrie...
per i debuttanti di GRASS: non scoraggiarsi; alcune cose del tipo
```
 v.category nodeSimple out=node layer=3 opt=add --overwrite
```
non si capiscono al volo (a meno di persone particolarmente portate). Questo è normale: rifare alcune volte un esercizio non è una perdita di tempo.

Tutti i comandi usati in questa esercitazione sono contenuti in uno script scaricabile da da questo link.

Full planar topology in GRASS (Prima parte)

2009-01-06T23:23:00.005+01:00

Un breve richiamo: la Full Planar Topology

La Full Planar Topology è un modello topologico poligonale che soggiace ad alcune regole fondamentali:

Connettività:: gli archi sono tra loro connessi con un nodo in comune; ogni arco inizia e termina con un nodo (ogni arco ha due nodi); tutte le intersezioni tra archi sono nodi.
Contiguità:: poligoni adiacenti hanno un solo arco in comune; ogni arco divide due aree.
Contenimento:: un poligono chiuso da archi ha un'area misurabile; serie di archi definiscono perimetri e aree chiuse. Ogni area è circondata da archi e nodi. Ogni nodo è contenuto in aree.

Poste queste regole, il vettoriale polyTopology della location simpleloc rappresentato qui di seguito, dovrebbe avere le tabelle riportate sotto:

Le entità geometriche fondamentali: nodi, archi, poligoni. Il verso degli archi si riferisce al verso di digitazione.

Id (cat) nodo	Archi
1	a1
2	a3,a5,a9
3	a8,a9,a10
4	a2,a3,a10
5	a4,a5,a6,a11
6	a4,a6,a7,a12
7	a11,a12,a13
8	a2,a7,a8,a13

Node topology.

Id arco	L_Poly	R_Poly
a1	polyA	polyC
a2	polyA	polyB
a3	polyA	polyD
a4	universe	polyE
a5	polyA	universe
a6	polyA	universe
a7	polyF	polyA
a8	universe	polyB
a9	universe	polyD
a10	polyB	polyD
a11	universe	polyE
a12	polyF	polyE
a13	polyF	universe

Arc Topology

Id (cat) poligono	label	Archi
1	polyA	a1,a2,a3,a5,a6,a7
2	polyB	a2,a8,a10
3	polyC	a1
4	polyD	a3,a9,a10
5	polyE	a4,a11,a12
6	polyF	a7,a12,a13

Polygon Topology

Obbiettivi

GRASS è un GIS topologico; le 3 tabelle fondamentali viste sopra devono quindi essere facilmente ottenibili.

Il comando

 v.build map=polyTopology opt=dump

consente la visualizzazione di tutti gli elementi della topologia del vettoriale polyTopology.

Lo scopo di questo esercizio è didattico: nella costruzione delle tre tabelle topologiche, si condurrà una simulazione senza ricorrere a v.build; questo consentirà di esplorare al meglio alcuni tools che interagiscono con le entità geometriche fondamentali (nodi,archi,poligoni).

Set di dati per l'esercitazione

La copertura poligonale scelta per l'esercitazione è polyTopology, contenuta nella location simpleloc.

Questa location è stata creata appositamente per fini didattici; per avere maggiori notizie vedere questo precedente blog.

La location simpleloc può essere scaricata da qui.

Procedura

Si creerà per prima la tabella Arc Topology, quindi la Node Topology e da ultimo la Polygon Topology.

Le fasi logiche sono le seguenti:

creazione di un layer di archi (=boundaries) numerato univocamente;
agli attributi degli archi vengono aggiunti gli identificativi (=cat) dei poligoni posti a destra e a sinistra di ciascun elemento;
creazione di un layer di nodi (=punti) numerati univocamente;
agli attributi dei nodi vengono aggiunti gli identificativi degli archi ai quali sono uniti;
dalla combinazione delle tabelle dei nodi e degli archi si ottiene la Polygon Topology.

Arc Topology

I boundaries per default non hanno una categoria. Si crea un nuovo layer aggiungendo un valore di categoria univoco:

v.category in=polyTopology type=boundary \
 layer=2 out=arcs option=add cat=1 --overwrite

Agli attributi del layer appena creato si aggiunge una serie di campi: alcuni conterrano solo testo e serviranno per la visualizzazione. Inoltre alcuni campi testuali offorno l'occasione per fare esercizio anche con l'SQL.

v.db.addtable arcs layer=2 columns="cat integer, \
 arcLabel char(3), \
 catPolyLeft integer, \
 polyLeftLabel char(5), \
 catPolyRight int, \
 polyRightLabel char(5), \
 arcLength double precision"

Si inizia a popolare il campo arcLabel:

v.db.update map=arcs layer=2 column=arcLabel qcolumn="'a'||cat"

L'operatore SQL || concatena le stringhe; in questo caso concatena la stringa a con il contenuto del campo cat.

Con il tool v.to.db opt=sides si ottengono le categorie dei poligoni rispettivamente a sinistra e a destra di ciascun arco (boundary):

v.to.db arcs layer=2 option=sides col=catPolyLeft,catPolyRight

Si calcola inoltre la lunghezza di ciascun arco:

v.to.db arcs layer=2 option=length col=arcLength

Viene ora il momento di usare l'SQL: sulla base delle categorie delle aree ottenute con v.to.db opt=sides si popolano i due campi polyLeftLabel e polyRightLabel:

echo "UPDATE arcs_2 SET polyLeftLabel=IFNULL( \
 (SELECT arcs.label FROM arcs \
 WHERE arcs.cat=arcs_2.catPolyLeft),'universe')" | db.execute

echo "UPDATE arcs_2 SET polyRightLabel=IFNULL( \
 (SELECT arcs.label FROM arcs \
 WHERE arcs.cat=arcs_2.catPolyRight),'universe')" | db.execute

L'operazione effettuata merita qualche spiegazione; sul layer 1 del vettoriale arcs vi sono gli attributi delle aree:

v.db.select arcs layer=1
cat|label|GRASSRGB
1|polyA|141:211:199
3|polyC|190:186:218
2|polyB|255:255:179
4|polyD|251:128:114
5|polyE|128:177:211
6|polyF|253:180:98

Nei campi catPolyLeft e catPolyRight del layer2 sono state memorizzate le categorie delle aree del layer 1: esiste quindi una relazione biunivoca fra i due campi del layer 2 e il campo cat del layer1. Grazie a questa relazione si può stabilire un join fra le due tabelle imponendo la condizione WHERE arcs.cat=arcs_2.catPolyLeft e WHERE arcs.cat=arcs_2.catPolyRight.

Se tali condizioni non vengono soddisfatte significa che il valore della cat del layer 1 non ha corrispondenze nel layer 2: è il caso del poligono universe (valore -1), l'area che, nel modello topologico Full Planar Topology avvolge idealmente ogni altro oggetto.

Dopo queste operazioni gli attributi sul layer 2 sono:

v.db.select arcs layer=2
cat|arcLabel|catPolyLeft|polyLeftLabel|catPolyRight|polyRightLabel|arcLength
1|a1|3|polyC|1|polyA|6.242641
2|a2|2|polyB|1|polyA|3
3|a3|4|polyD|1|polyA|6.708204
4|a4|5|polyE|-1|universe|3.414214
5|a5|-1|universe|1|polyA|5.019765
6|a6|-1|universe|1|polyA|2.828427
7|a7|1|polyA|6|polyF|4.162278
8|a8|2|polyB|-1|universe|5.162278
9|a9|4|polyD|-1|universe|10.261297
10|a10|4|polyD|2|polyB|3
11|a11|5|polyE|-1|universe|5.828427
12|a12|5|polyE|6|polyF|2.414214
13|a13|-1|universe|6|polyF|7.071068

Si verifichi anche la situazione delle categorie sui due layers:

v.category arcs opt=report
Layer/table: 2/arcs_2
type count min max
point 0 0 0
line 0 0 0
boundary 13 1 13
centroid 0 0 0
area 0 0 0
all 13 1 13
Layer/table: 1/arcs
type count min max
point 0 0 0
line 0 0 0
boundary 0 0 0
centroid 6 1 6
area 0 0 0
all 6 1 6

A questo punto la costruzione della tabella Arc Topology è conclusa:

v.db.select arcs layer=2 \
 col=arcLabel,polyRightLabel,polyLeftLabel
arcLabel|polyRightLabel|polyLeftLabel
a1|polyA|polyC
a2|polyA|polyB
a3|polyA|polyD
a4|universe|polyE
a5|polyA|universe
a6|polyA|universe
a7|polyF|polyA
a8|universe|polyB
a9|universe|polyD
a10|polyB|polyD
a11|universe|polyE
a12|polyF|polyE
a13|polyF|universe

Si procede ora alla visualizzazione :

./sl_draw.sh
d.font.freetype font=timesbd
d.vect -a arcs layer=1 display=shape 
d.vect arcs layer=1 display=attr attrcol=Label llayer=1 \
 color=black lsize=10 xref=center yref=bottom
d.vect arcs layer=2 display=shape,attr,dir attrcol=arcLabel \
 type=boundary llayer=2 \
 xref=center yref=center color=grey \
 lsize=10 bgcolor=white bcolor=black

Il vettoriale arcs: sul layer 1 mantiene le aree; sul layer 2 gli archi (=boundaries) con gli attributi calcolati.

Trovare i contorni duplicati con GRASS

2008-12-28T00:09:00.006+01:00

Obbiettivi

Con questa esercitazione si vuole dimostrare come GRASS possa individuare quali e quanti sono, in una copertura poligonale non topologica, i contorni (boundaries) duplicati.

Si farà ricorso a moduli di alto livello, senza dover accedere ai vettoriali a basso livello, come nel caso del'impiego di v.build (alto e basso livello vengono qui usati per riferirsi al livello di interfacciabilità verso l'utente).

Set di dati per l'esercitazione

La copertura poligonale scelta per l'esercitazione è polyNoTopology, contenuta nella location simpleloc.

Questa location è stata creata appositamente per fini didattici; per avere maggiori notizie vedere questo precedente blog.

La location simpleloc può essere scaricata da qui.

polyNoTopology è un vettoriale non topologico (modello CAD o modello spaghetti) con le seguenti caratteristiche:

i contorni dei poligoni adiacenti sono duplicati;
i contorni, se si intersecano, non generano nodi.

Prime verifiche

Avviare GRASS ed entrare nella location simpleloc; si dovrebbero trovare i seguenti vettoriali:

GRASS 6.4.svn (simpleloc):~ > g.list vect
----------------------------------------------
vector files available in mapset <user1>:
node_grid polyBreakCat polyNoTopology punti xylabels
polyBoundCat polyClean polyTopology vgrid
polyBreak polyFromArc poly_arc vrastergrid

Si può verificare che vi sono lati duplicati controllando le coordinate dei vertici dei contorni semplicemente con v.out.ascii:

v.out.ascii in=polyNoTopology \
 format=standard | sed -n '/^B/,/^C/p' | sed -n '/^C/!p'

Leggendo le coordinate in output si può verificare che i lati in comune sono digitati due volte. Il pipe ('|') verso sed non è strettamente necessario e serve solo per filtrare i dati non utili in questo contesto.

Cercare i contorni duplicati

I contorni di polyNoTopology, oltre ad essere duplicati in corrispondenza delle adiacenze tra poligoni, non sono spezzati (break).

La procedura per trovare i duplicati può essere riassunta nei seguenti passi:

spezzare in boundaries (= contorni dei poligoni), inserendo un nodo nei punti di incrocio;
numerare in modo univoco i boundaries;
rimuovere i boundaries duplicati;
attribuire nuovamente una numerazione univoca ai boundaries senza duplicati
interrogare, dallo strato dei boundaries senza duplicati, le geometrie dello strato con i duplicati: per ciascun elemento geometrico del primo ottenere quanti e quali elementi del secondo esistono: se si ottiene un valore pari ad 1 il boundaries in quel tratto non è duplicato; se si ottiene un valore maggiore di 1 il boundaries è duplicato.

Per spezzare i boundaries si usa il modulo v.clean con il tool break:

v.clean in= vectInput out=polyBreak tool=break type=boundary --overwrite

Si può facilmente verificare che i boundaries non hanno categoria; nell'output di v.category opt=report, infatti, vengono riportate le categorie per i soli centroidi:

v.category in=polyBreak opt=report
Layer/table: 1/polyBreak
type count min max
point 0 0 0
line 0 0 0
boundary 0 0 0
centroid 6 1 6
area 0 0 0
all 6 1 6

Si aggiunge allora la categoria ai boundaries spezzati (break) creando il layer 2 ; considerato che i centroidi hanno categoria da 1 a 6 , per non confondere i valori, si decide di numerare i boundaries da 100 in poi:

v.category in=polyBreak type=boundary out=polyBreakCat \
 layer=2 opt=add cat=100 --overwrite

La situazione dei layers e delle categorie può essere evidenziata ancora una volta con v.category opt=report:

v.category in=polyBreakCat opt=report
Layer: 2
type count min max
point 0 0 0
line 0 0 0
boundary 18 100 117
centroid 0 0 0
area 0 0 0
all 18 100 117
Layer/table: 1/polyBreakCat
type count min max
point 0 0 0
line 0 0 0
boundary 0 0 0
centroid 6 1 6
area 0 0 0
all 6 1 6

La situazione è quindi la seguente: sul layer 1 di polyBreakCat vi sono i centroidi con categorie da 1 a 6 ; sul layer 2 vi sono i buondaries numerati, in modo univoco, da 100 a 117 .

Occorre ora rimuovere i duplicati da polyBreakCat:

v.clean in=polyBreakCat out=polyClean \
 type=boundary tool=rmdupl --overwrite

GRASS ha rimosso i pezzi di boundaries duplicati; il vettoriale ottenuto (polyClean) conserva sul layer 1 le categorie dei centroidi e sul layer 2 le categorie di tutti i boundaries, compresi quelli duplicati; ed è proprio quest'ultima caratterisitca il segreto di questo eserzio.

Se infatti ora si rinumerano i boundaries, verrano rinumerati solo ed esclusivamente i boundaries senza i duplicati ma esisterà una relazione uno a molti tra il layer 2 e il layer 3 .

v.category in=polyClean out=polyBoundCat \
 type=boundary layer=3 opt=add cat=100 --overwrite
v.category in=polyBoundCat opt=report
Layer/table: 2/polyBoundCat_2
type count min max
point 0 0 0
line 0 0 0
boundary 18 100 117
centroid 0 0 0
area 0 0 0
all 18 100 117
Layer/table: 3/polyBoundCat_3
type count min max
point 0 0 0
line 0 0 0
boundary 13 100 112
centroid 0 0 0
area 0 0 0
all 13 100 112
Layer/table: 1/polyBoundCat
type count min max
point 0 0 0
line 0 0 0
boundary 0 0 0
centroid 6 1 6
area 0 0 0
all 6 1 6

Riassumendo:

sul layer 1 : i centroidi (categorie da 1 a 6 );
sul layer 2 : i boundaries spezzati (break) ma con i duplicati (categorie da 100 a 117 );
sul layer 3 : i boundaries spezzati (break) e senza dupicati (categorie da 100 a 112 ).

Si aggiungono ora le tabelle al layer 2 e 3 di polyBoundCat:

v.db.addtable polyBoundCat layer=2

v.db.addtable polyBoundCat layer=3 \
 col="cat integer, count integer,listarcs varchar(12)"

GRASS offre un potente modulo per interrogare le geometrie e caricare i dati dell'interrogazione nella tabella degli attributi: v.to.db.

Nel modello vettoriale di GRASS, per ogni primitiva geometrica possono venir memorizzati i numeri dei layers e, per ciascun numero di layer, il numero della categoria (difficilmente si capirà il meccanismo di questo esercizio se non si ha chiaro questo concetto; vedere il manuale di GRASS).

Viene ora il momento di interrogare le geometrie del layer 2 , popolando i due campi count e listarcs:

v.to.db map=polyBoundCat option=query type=boundary \
 layer=3 qlayer=2 column=count qcolumn="count(*)"
v.to.db map=polyBoundCat option=query type=boundary \
 layer=3 qlayer=2 column=listarcs qcolumn="group_concat(cat)"

Nella seconda operazione si usa il comando di aggregazione group_concat supportato nelle recenti versioni di SQLite: questo comando concatena i valori trovati, creando così l'elenco delle categorie dei boundaries.

Si passa ora alla visualizzazione della tabella:

v.db.select polyBoundCat layer=3
cat|count|listarcs
100|1|100
101|1|105
102|1|106
103|1|107
104|1|111
105|1|112
106|1|115
107|1|117
108|2|101,103
109|2|102,104
110|2|108,109
111|2|110,113
112|2|114,116

GRASS ha correttamente indviduato quali e quanti boundaries erano duplicati:

nella colonna count, per ogni record (=per ogni geometria) del layer 3 sono stati contati quanti elementi sono presenti nel layer 2 ;
nella colonna listarcs, per ogni record (=per ogni geometria) del layer 3 sono stati elencati quali elementi sono presenti nel layer 2 ; in questo caso sono state elencate le categorie corrispondenti ai boundaries contenuti nel layer 2 .

Una visualizzazione grafica renderà più chiare le idee; si può usare la linea di comando:

# # # disegno
./sl_draw.sh
d.vect polyBoundCat layer=3 display=shape,cat color=grey llayer=3 \
 xref=center yref=center lcolor=aqua \
 lsize=10 bgcolor=white bcolor=black
# # # output del file png 
d.out.file out=polyBoundCat_l3 format=png --overwrite 
./sl_draw.sh
d.vect polyBoundCat layer=2 display=shape,cat color=grey llayer=2 \
 xref=center yref=center lcolor=red \
 lsize=10 bgcolor=white bcolor=black
# # # output del file png 
d.out.file out=polyBoundCat_l2 format=png --overwrite

Tutto il codice utilizzato in questo esercizio può esere scaricato da qui

Il risultato delle ultime righe di comando dovrebbe essere la visualizzazione dei due layer di polyBoundCat in due separati monitors:

È pure possibile usare la comoda interfaccia grafica wxpython di GRASS:

Il caso dell'import di shapefiles

Tutta la procedura vista può essere utilmente applicata al caso dell'import di shapefiles in GRASS.

Lo shapefile è una struttura dati di tipo non topologico. Si provi allora ad esportare in shapefile la copertura poligonale (topologicamente corretta) della location simpleloc denominata polyTopology con v.out.ogr:

v.out.ogr -c input=polyTopology type=area dsn=~/tmp

L'opzione -c fa in modo che nell'export siano comprese solo le features geometriche con categorie (nel caso di polyTopology tutte le aree). Si può ora reimportare in GRASS lo shapefile:

v.in ogr dsn=~/tmp layer=polyTopology out=polyFromShp snap=-1 --o

L'esercizio visto può essere ripetuto con il vettoriale ottenuto dall'import (polyFromShp) addivenendo ai medesimi risultati.

Il modulo v.clean offre dei tools appositi per la pulizia dei vettoriali importati in GRASS da strutture dati (come gli shapefiles) non topologici: rmdupl e bpol; il primo rimuove i duplicati; il secondo spezza (=break) in modo topologico il contorno dei poligoni.

Conclusioni

Questa esercitazione permette di capire quali attenzioni devono essere poste quando si importa in un GIS topologico come GRASS una struttura dati non topologica.

Il modulo v.clean rappresenta il passaggio obbligato attraverso il quale pulire i vettoriali importati in GRASS: si è vista la possibilità di rimozione dei duplicati (rmdupl), l'inserimento di nodi nei punti di intersezione (break); si è pure accennato all'apposito tool bpol (per i polgioni).

Il modulo v.to.db consente di interrogare le geometrie tra layers diversi dello stesso vettoriale con l'opzione query.

Sono stati trattati solo due dei diversi tipi di errori di topologia. Per altre anomalie, quali pseudonodi, dangles, weird polygons, sliver polygons, GRASS offre dei rimedi più o meno efficaci, che si esamineranno in dettaglio in prossime esercitazioni.

Piani Forestali e wget

2008-12-13T16:03:00.001+01:00

Visto che per me è ormai diventata ormai una fissazione, mi appunto il comando per scaricare, tutti insieme e in un colpo solo, i Piani Forestali della Regione Piemonte (tutti i comuni e tutte le aree forestali) :

wget -r -nd -nH --level=3 --accept=zip \ http://www.regione.piemonte.it/montagna/foreste/pianifor/forestale.htm

Ricordarsi di registrarsi... prima.

ogr2ogr, codici epsg e Piani Forestali del Piemonte

2008-12-13T15:47:00.008+01:00

In un precedente post ho usato ogr2ogr per convertire gli shapefiles dei Piani Forestali della Regione Piemonte da ED50 UTM 32N a WGS84 UTM 32N.
Quando si parla di conversione tra datum diversi divento diffidente: forse a causa della scarsa conoscenza della materia e delle ripetute fregature avute in passato...e nel presente!
In ogni caso voglio semplicemente provare ad usare ogr2ogr con diversi settaggi, tanto per chiarirmi le idee.
La sintassi (semplificata) di ogr2ogr per la conversione tra sistemi di riferimento diversi (di shapefiles) è la seguente:
ogr2ogr -s_srs (descrizione del sistema di riferimento dell'input)-t_srs (descrizione del sistema di riferimento dell'output) (nome shapefile di output) (nome shapefile di input)

Per quanto riguarda la descrizione del sistema di riferimento dell'imput (nel mio caso: ED50 UTM 32N) le possibilità che ho preso in considerazione sono:

una variabile che immagazzina una lunga stringa wkt, come quella riportata qui:

wkt_ed50='PROJCS["UTM Zone 32, Northern Hemisphere", GEOGCS["international", DATUM["European_Datum_1950", SPHEROID["International_1924",6378388,297], TOWGS84[-87.000,-98.000,-121.000]], PRIMEM["Greenwich",0], UNIT["degree",0.0174532925199433]], PROJECTION["Transverse_Mercator"], PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",9], PARAMETER["scale_factor",0.9996], PARAMETER["false_easting",500000], PARAMETER["false_northing",0],UNIT["meter",1]]'
Per ottenere questo wkt ho semplicemente creato una Location di GRASSS in ED50 UTM32N e da shell ho digitato g.proj -w. Lo stesso risultato potrebbe essere ottenuto da qui (fare attenzione ai parametri di traslazione TOWGS84...);
il codice epsg:
'+init=epsg:23032'
il codice epsg (23032) e i parametri di traslazione: '+init=epsg:23032 +towgs84=-87.000,-98.000,-121.000'
le stesse informazioni di cui sopra ma in formato proj:
'+proj=utm +zone=32 +ellps=intl +units=m +towgs84=-87.000,-98.000,-121.000'

Per quanto riguarda l'output in WGS84 UTM 32N ho semplicemente settato al codice epsg 32632.
Ho quindi scritto questo semplice script per testare il tutto su uno shapefile dei Piani Forestali della Regione Piemonte:





#!/bin/sh


dir_work="$HOME/pft"


shpfile="004017"


suffix="_wkt"


ogr2ogr -s_srs "$wkt_ed50" \ 

    -t_srs '+init=epsg:32632' \


    "$dir_work/${shpfile}${suffix}.shp" \ 

    "$dir_work/${shpfile}.SHP"
suffix="_epsg"




ogr2ogr -s_srs '+init=epsg:23032'\

    -t_srs '+init=epsg:32632' \


    "$dir_work/${shpfile}${suffix}.shp" \ 

    "$dir_work/${shpfile}.SHP"
suffix="_epsgtowg"




ogr2ogr -s_srs '+init=epsg:23032 +towgs84=-87.000,-98.000,-121.000' \


    -t_srs '+init=epsg:32632' \ 

    "$dir_work/${shpfile}${suffix}.shp" \


    "$dir_work/${shpfile}.SHP"
suffix="_defepsg"




ogr2ogr -s_srs '+proj=utm +zone=32 +ellps=intl +units=m
+towgs84=-87.000,-98.000,-121.000' \ 

    -t_srs '+init=epsg:32632' \


    "$dir_work/${shpfile}${suffix}.shp" \ 

    "$dir_work/${shpfile}.SHP"
suffix="traspunto"




wine /usr/local/traspunto/traspunto.exe shp pia pia \

    "${dir_work}/${shpfile}" \

    "${dir_work}/${shpfile}${suffix}" ED50 WGS84
32 32

dove nella variabile $wkt_ed50 è contenuto il file wkt visto sopra.

Da notare che uso anche traspunto, disponibile qui, con wine.

Questa è una prova (empirica) nella quale apprezzo visivamente i risultati; non ho qui l'interesse di fare alcuna analisi statistica del tipo: esiste una differenza statistica (e se si, quale?)nel posizionamento spaziale di una serie di punti trasformati con traspunto rispetto a ogr2ogr?

Una volta eseguito lo script visualizzo il tutto in OpenJUMP (ma, ovviamente, vanno bene anche Quantum Gis, gvSIG, udig, ecc).

Il primo confronto interessante è quello tra traspunto (linee in rosso) e ogr2ogr con '+init=epsg:23032 +towgs84=-87.000,-98.000,-121.000' (linee in verde).

Per apprezzare la differenza, nella figura ho settato una griglia di sfondo con una maglia di 10 metri.

Il secondo confronto riguarda i due settagggi di ogr2ogr: il primo con il solo codice epsg (in blu) (-s_srs '+init=epsg:23032'); il secondo (in verde) con il codice e i parametri di traslazione ('+init=epsg:23032 +towgs84=-87.000,-98.000,-121.000'). Ho mantenuto la griglia a 10 m.

Com'era prevedibile non esistono differenze tra il risultato della conversione effettuata con il wkt, con '+init=epsg:23032 +towgs84=-87.000,-98.000,-121.000' e con
-s_srs '+proj=utm +zone=32 +ellps=intl +units=m +towgs84=-87.000,-98.000,-121.000'.

Conclusioni

Una prova rigorosa su modalità diverse di conversione può solo essere condotta con una serie di punti dei quali si conoscono a priori con esattezza le coordinate nei diversi sistemi di interesse.
Questa è una semplice prova per fissare alcune cose:

non fermarsi al codice epsg: cercare, se vi sono, i parametri di rototraslazione relativi alla zona di interesse;
nella definizione di una nuova Location porsi il problema di settare i parametri di rototraslazione (vedere g.proj);
l'uso di traspunto può essere un'alternativa;
l'uso del software Verto unitamente ai grigliati (in vendita presso l'Istituto Geografico Militare ) potrebbe costituire la migliore tra le alternative (?) ...