Il formato Shapefile è ovviamente superiore al GPX, però questa conversione può tornare utile ad esempio per caricare in JOSM uno shapefile da usare come base per il disegno.

Questa ricetta funziona solo per la conversione di uno shapefile di tipo line, per quelli di tipo polygon ci vuole un po' di pre-processing ad esempio in GRASS.

La cosa più importante da scoprire è il sistema di riferimento utilizzato nello shape, per la cartografia italiana è abbastanza consueto trovare il Gauss-Boaga Roma40. Per avere questa informazione si deve ricorrere alla documentazione che dovrebbe accompagnare lo shapefile.

Esiste il programma gpx2shp, ma contrariamente a quelle che si legge sulla home page non consente la trasformazione inversa (shp2gpx).

Nell'esempio che segue si converte uno shapefile Roma40 in un file GPX WGS84. Lo shapefile originale è composto da tre file: VX.dbf, VX.shp e VX.shx.

Il primo programma da utilizzare è ogr2ogr. La home page del programma è www.gdal.org, mentre Debian lo fornisce nel pacchetto gdal-bin.

Purtroppo per la conversione da Roma40 a WGS84 è necessario dare ad ogr2ogr un parametro un po' complicato; si tratta del parametro di correzione per la libreria PROJ, usata da ogr2ogr. Per comodità lo scriviamo in un file di nome roma40.srs, mettendo tutto su una riga sola:

+proj=tmerc +lat_0=0 +lon_0=9 +k=0.999600 +x_0=1500000 +y_0=0 +ellps=intl +units=m
    +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68

Possiamo eseguire la riproiezione con il comando:

ogr2ogr -f "ESRI Shapefile" -s_srs "roma40.srs" -t_srs "EPSG:4326" wgs84_shp VX.shp

Nella directory wgs84_shp troveremo il nuovo shapefile nel sistema di riferimento WGS84 (identificato dal codice EPSG 4326).

Adesso si usa il programma GPSBabel per convertire da shapefile a GPX. Il comportamento predefinito di GPSBabel sarebbe quello di convertire le linee dello shape in rotte GPX, noi invece vogliamo delle tracce, ecco il motivo dei due parametri transform:

cd wgs84_shp
gpsbabel -i shape -f VX -x transform,trk=rte -x nuketypes,routes -o gpx -F VX.gpx

Procedura decisamente più complessa, utile però a capire un po' la logica dell'accoppiata QGIS+GRASS.

1. Aprire lo shapefile in QGIS e attribuire il sistema di riferimento giusto (Roma40).
2. Creare un mapset GRASS Roma40 e importare lo shape da QGIS.
3. Creare un mapset GRASS WGS84 e importare (riproiettando) il layer GRASS Roma40.
4. Esportare il layer GRASS in formato shapefile.
5. Usare GPSBabel per convertire da shapefile a GPX.

Controllare e risistemare!

1. QGIS, menu File, Nuovo progetto.
2. Dal menu Impostazioni, Proiezione personalizzata definire una proiezione di nome Gauss-Boaga Roma40 con i parametri che abbiamo visto nel paragrafo precedente (file roma40.srs).
3. Riavviare QGIS per far rileggere l'elenco delle proiezioni.
4. Menu Plugin, GRASS, Nuovo mapset.
   • Il database di solito è $HOME/grass.
   • Creare una nuova location: osm_arezzo_roma40
   • Selezionare una prioezione: User defined Gauss-Boaga Roma40
   • Selezionare una regione: Italy
   • Dare un nome al mapset: arezzo
   • Il mapset appena creato dovrebbe diventare quello in uso, altrimenti usare: Plugin, GRASS, Apri mapset.
5. Aggiungere in QGIS lo shapefile originale col pulsante Aggiungi un layer vettoriale. In Proprietà del layer (click destro), Generale, Sistema di Riferimento Spaziale scegliere User defined Gauss-Boaga Roma40.
6. Importare lo shapefile in GRASS; dal pulsante Apri strumenti GRASS avviare il modulo File, Import, Import vector, v.in.ogr (Import OGR/PostGIS vector layer).
   • Indicare come origine il layer appena aggiunto in QGIS.
   • Indicare il nome che questo layer avrà in GRASS, es: VX_GRASS_Roma40.
7. Creare un altro mapset: menu Plugin, GRASS, Nuovo mapset.
   • Location: osm_arezzo_wgs84
   • Proiezione WGS84 (EPSG 4326)
   • Selezionare una regione: Italy
   • Dare un nome al mapset: arezzo
   • Il mapset appena creato dovrebbe diventare quello in uso.
8. Convertire il layer vettoriale GRASS da una location all'altra; dal pulsante Apri strumenti GRASS avviare il modulo Vector, Develop map, Geometry management, Reproject vector from another location, v.proj:
   • Name of input vector map: VX_GRASS_Roma40
   • Path to GRASS database of input location: $HOME/grass/
   • Location containing input vector map: osm_arezzo_roma40
   • Mapset containing input vector map: arezzo
   • Name of output vector map: VX_GRASS_WGS84
9. Esportare il layer GRASS in formato shapefile:

v.out.ogr input=AO2k_line type=line dsn=AO2k_line_shp olayer=AO2k_line_shp layer=1 format=ESRI_Shapefile

Per una semplice conversione di coordinate si può utilizzare il comando cs2cs fornito dal pacchetto proj.

Questo esempio converte una coppia di coordinate dal sistema EPSG:3003 (Monte Mario, Italy zone 1 - Roma40) al sistema EPSG:4326 (WGS84). La cosa è resa complicata perché la definizione di EPSG:3003 fornita da Proj non include i parametri di trasformazione verso WGS84, che vanno quindi aggiunti a mano:

echo "1674257 4846294" | cs2cs \
    +init=epsg:3003 \
    +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68 \
    +to +init=epsg:4326 -f "%.8f"
11.16398206     43.74906364 45.83517126

Il comando restituisce una terna di numeri (le coordinate x, y e z), cioè anche l'altezza.