giovedì 11 dicembre 2008

Piani Forestali, ogr2ogr e conversione di coordinate

Da qualche tempo lavoro con i Piani Forestali della Regione Piemonte che possono essere scaricati da qui previa registrazione.
È possibile scaricare per singoli Comuni oppure per aree forestali.
Lavorare con un gran numero di shapefiles è abbastanza semplice, quando si capisce la possibilità di iterare all'interno di una directory contenente gli shapefiles medesimi.
Ecco i primi passaggi che faccio nell'ipotesi di aver scaricato un'area forestale.

L'iterazione ha sempre la stessa struttura base: for files in directory; do qualcosa; done

Ad esempio:

laptop@pcgis:~/pft$ for f in *.*; do echo $f; done

004002.ZIP

004006.ZIP

004025.ZIP

004079.ZIP

004093.ZIP

... etc ...



Appena scaricati, sul mio pc Linux-Ubuntu i files zippati non hanno permessi di scrittura-lettura


laptop@pcgis:~/pft$ for f in *.*; do echo $f; chmod ugoa+wr $f; done

004002.ZIP

004006.ZIP

004025.ZIP

004079.ZIP

004093.ZIP

... etc ...



A questo punto scompatto gli zip:

laptop@pcgis:~/pft$ for f in *.*; do echo $f;unzip $f; done004002.ZIP

Archive: 004002.ZIP

inflating: 004002.TXT

inflating: 004002.DBF

inflating: 004002.SHX

inflating: 004002.SHP

004006.ZIP

Archive: 004006.ZIP

inflating: 004006.TXT

inflating: 004006.DBF

inflating: 004006.SHX

inflating: 004006.SHP

004025.ZIP

Archive: 004025.ZIP

inflating: 004025.TXT

inflating: 004025.DBF

inflating: 004025.SHX

inflating: 004025.SHP

... etc ...




Gli shapefiles dei PFT della Regione Piemonte sono in ED50 UTM32 Nord che corrisponde ad un codice epsg pari a 23032 (parametri per la rototraslazione in WGS84: TOWGS84=-87.000,-98.000,-121.000); siccome la direttiva che ho ricevuto è quella di lavorare in WGS 84 UTM 32 Nord, devo fare una conversione tra datum

laptop@pcgis:~/pft$ suffix="_wgs84"



laptop@pcgis:~/pft$ for f in *.SHP; do name=${f%.*}; ogr2ogr \

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

-t_srs '+init=epsg:32632' "${name}${suffix}.shp" $f; done



Con poche righe di comando ho ottenuto nella cartella pft una lunga serie di shapefiles nominati con il vecchio nome dello shapefile in ED50 UTM 32N più un suffisso _wgs84. Tutti gli shapefiles generati da ogr2ogr hanno il file prj che descrive il sistema di riferimento (WGS84 UTM 32N)

In figura si può vedere il medesimo shapefile riferito al sistema ED50 UTM 32 N (in marrone) e il risultato della traslazione effettuata con ogr2ogr (in viola).



From Grass, gis e dintorni


Ovvio che il tutto può essere inglobato in uno script, per migliorarne la leggibilità e la riusabilità

Per saperne di più su ogr2ogr si può consultare il suo manuale qui.

I comandi impartiti richiedono una spiegazione che sarà oggetto di un prossimo post.

Nessun commento:

Posta un commento