Come ritagliare porzioni di mappe...

S

sun75

Guest
GUIDA: Come ritagliare porzioni di mappe...

Ecco qui una breve guida su come poter ritagliare delle porzioni di mappe gia' georeferenziate, senza esportarsi una jpeg dell'area interessata e ri-georeferenziare nuovamente, perdendo in precisione!

Il procedimento che qui illustro, lo uso da anni, per ritagliarmi dalle regioni_IGM il pezzetto di cartografia che mi serve nelle mie escursioni, senza dover caricare l'intera regione. Ha il vantaggio di non perdere mai la georeferenziazione originale, e quindi la precisione della carta "figlia" e' identica a quella originale!

Procediamo con ordine!

Software necessari:
  • OziExplorer o equivalente software cartografico sul quale rilevare le coordinate dell'area da ritagliare. Io uso Ozi, quindi le info qui contenute saranno riferite ad Ozi.
  • FWTools, compilation free di software GIS a linea di comando per manipolare la cartografia elettronica. Si scaricano da qui: FWTools: Open Source GIS/RS Binary Kit
Una volta che sia OziExplorer che gli FWTools sono installati correttamente sul PC, inizia la procedura.

Ammettiamo, nel nostro esempio, di voler ritagliare tutta l'area del Velino-Cafornia-Bicchero-Costone-Duchessa-Valle del Salto nel Parco Regionale Velino-Sirente.

Assunti:
  • La mia directory di lavoro e': T:\Cartografia\Data
  • La mia directory degli originali cartografici e': T:\Cartografia\IGM
1) Apriamo OziExplorer e l'FWtools Shell. Nella Shell ci spostiamo subito nella directory di lavoro!
2) Si carica Abruzzo_IGM_25K.map in OziExplorer
3) Va abilitato il reticolo su Ozi in questo modo:
  • Menu' "Mappa" -> Impostazioni Linea Reticolo
  • Altro Reticolo
  • Flag su "Reticolo Attivo"
  • Intervallo Linea: "Nessun reticolo" (preferibile, oppure "1Km" se volete vederlo coincidere con quello reale della carta)
  • Colori Linea (entrambi): "blu"
  • Flag su "Scala Automatica"
  • Intervallo Etichette: "1Km"
  • Colore etichetta: nero
  • Colore sfondo: giallo
  • Dimensione font: 15
  • Visualiz. Etichette: "su tutti gli schermi"
  • Numeri: "niente metri"
  • Flag su "mostra linea margini"
  • Quindi "Salva" -> Default -> Chiudi
Ora vedrete il reticolo e le coord. UTM stampate a lato
Se il reticolo blu non coincide con quello disegnato sulla carta, dovete impostare come Datum "ED50" (European 1950) e non WGS84 che e' il default di Ozi. Sistemato il Datum, i due reticoli _devono_ coincidere!

4) Ora si ricavano le coordinate dei 4 lati dell'area di nostro interesse: a noi interessano solo l'angono superiore sinistro e quello inferiore destro.
Nel nostro caso le coordinate sono:
  • Sup. sx: 357000 4673000
  • Inf. dx: 371000 4663000
5) Ora, sulla shell FWTools, ritagliamo la mappa in due passaggi. Il nostro obiettivo e' di avere un nuovo Ecw da quello iniziale, e pertanto dobbiamo passare per un formato intermedio, in quanto non e' possibile "convertire" da ecw a ecw direttamente: il convertitore si aspetta sempre un formato finale diverso da quello iniziale. Quindi scegliamo il GeoTiff come formato intermedio, visto che e' un formato senza compressione, ovvero non perdiamo definizione tra i due salti di conversione. Infatti quando ripasseremo da GeoTiff a Ecw, quest'ultimo compresso, si perdera' un minimo di definizione grafica, anche se impercettibile a occhio!

5a) Da ecw a tif con ritaglio e georeferenza applicata. Il comando di conversione e' il seguente:

Codice:
gdal_translate.exe -of GTiff -co LARGE_OK=YES -co TILED=YES -co TFW=YES -a_srs EPSG:32633 -projwin {up sx} {down dx} {src} {dest}

Dove:
  • -of GTiff e' il formato di destinazione
  • -co LARGE_OK=YES serve a superare il limite dei 500Mb nella trattazione dei file cartografici (si intende applicato al file "raw" da trattare e non a quello compresso!! Quindi Abruzzo.ecw che sul disco e' 477Mb, Raw sono 10,01GB!!)
  • -co TILED=YES per la generazione di immagini tiff omogenee
  • -co TFW=YES per il salvataggio dell'ESRI World File (dati di georeferenziazione)
  • -a_srs EPSG:32633 se manca questo parametro, di default, il file di destinazione viene "strippato" della georeferenza. La chiave e' tutta qui! Per capire il parametro: a_srs indica a gdal di fare l'override dei dati di datum e proiezione e di ricalcolarli e riferenziare il ritaglio effettuato. E' un discorso molto ampio, ma il modo piu' semplice di farlo e' referenziando la nuova immagine usando il modello EPSG, i cui codici sono tutti scritti nel file pcs.csv nella directory "data" sotto l'installazione di FWTools in C:\Programmi o dove l'avete installato. Se aprite quel file, noterete che, in fondo, 32632 e 32633 corrispondono ai codici EPSG per il Datum WGS84, Zone N-UTM32 e N-UTM33 (per l'Italia) dove N sta' per emisfero Nord e non per la fascia (S e T in Italia). Perche' WGS84? Perche' nel nostro caso stiamo lavorando con cartografia IGM 1:25K georeferenziata con Datum WGS84 (lo potete vedere con il comando "gdalinfo {file}", prime righe in alto: "PROGCS" e "GEOCS"). Il fatto che il reticolo sia ED50 e' un dettaglio irrilevante ai fini della traslazione dei dati di referenza: quelli sono e rimangono WGS84!
Codice:
Driver: ECW/ERMapper Compressed Wavelets
Files: t:\cartografia\igm\Abruzzo_IGM_25K.ecw
Size is 64000, 56000
Coordinate System is:
PROJCS["NUTM33",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9108"]],
        AXIS["Lat",NORTH],
        AXIS["Long",EAST],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",15],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
  • -projwin {up sx} {down dx} niente altro che le coppie di coord. est e nord degli angoli superiore sx e inferiore dx del rettangolo che a noi interessa
  • {src} {dest} sono i file di origine e destinazione del processo
Quindi il nostro "fantastico" comando diviene:

Codice:
gdal_translate.exe -of GTiff -co LARGE_OK=YES -co TILED=YES -co TFW=YES -a_srs EPSG:32633 -projwin 357000 4673000 371000 4663000 T:\Cartografia\IGM\Abruzzo_IGM_25K.ecw velino.tif

Volendo si puo' verificare che il nuovo file sia correttamente georeferenziato con il seguente comando:

Codice:
gdalinfo velino.tif

Controllando la presenza dei tag "PROGCS" e "GEOCS"

5b) Ora dobbiamo convertire la nostra nuova carta in Ecw, e questo e' semplice semplice quanto il seguente comando:

Codice:
gdal_translate.exe -of ECW -co LARGE_OK=YES {src} {dest}

Dove:
  • -of ECW e' il formato di destinazione
  • -co LARGE_OK=YES serve a superare il limite dei 500Mb nella trattazione dei file cartografici (in questo caso potrebbe non servire se trattiamo file ridotti....ma metterlo non costa nulla!!)
  • {src} {dest} sono i file di origine e destinazione del processo
Nota: come potete vedere non ho inserito il "fondamentale" parametro "-a_srs EPSG:32633" in questo secondo comando...perche'? Semplice: nel primo caso trattavasi di un ritaglio con ricalcolo della georeferenza applicata alla nuova immagine e in quanto tale andava istruito gdal su come e cosa calcolare e come codificarlo; in questo caso, invece l'immagine non viene toccata, viene solo convertito il formato della mappa, e in quanto tale i dati di georeferenza rimangono assolutamente invariati e quindi trasferiti al nuovo formato "as-is"!

Quindi il nostro "fantastico" comando diviene:

Codice:
gdal_translate.exe -of ECW -co LARGE_OK=YES velino.tif velino.ecw

Una controllatina veloce con gdalinfo tanto per avere la certezza che i dati di georeferenza ci siano:

Codice:
gdalinfo velino.ecw

Controllando la presenza dei tag "PROGCS" e "GEOCS"

Bene, abbiamo la nostra nuova e fantastica carta ritagliata e georeferenziata "fotocopia" dell'originale! Ora .... usiamola!

6) Importazione e verifica in OziEplorer. Chiudiamo la mappa dell'Abruzzo, e importiamo la nuova:
  • Menu' File -> Chiudi Mappa
  • Menu' File -> Importa Mappa -> Mappa DRG Singola
  • Premere OK al Passo 1
  • Selezionare il file .ecw appena creato, quindi velino.ecw
  • Premere OK al Passo 2
  • Specificare la directory dove _creare_ il file .map di OziExplorer: generalmente nella stessa directory dov'e' l'.ecw...quindi nel mio caso "T:\Cartografia\Data"
  • Datum Mappa: WGS84
  • Proiezione Mappa: UTM
  • Zona Reticolo: 33
  • Emisfero: N (Nord)
  • Confemare premendo OK.
Dopo il messaggio di conferma, apparira' la nuova porzione di mappa perfettamente georeferenziata come l'originale!
Per verificarlo, basta andare su:
  • Menu' Mappa -> Impostazioni Linea Reticolo
  • Altro Reticolo -> Spunta su "Reticolo Attivo" (se prima avevate salvato con "Default" ora dovreste avere tutti gli altri campi gia' popolati)
  • Quindi "Chiudi"
Il Datum deve, ovviamente, essere sempre ED50 (European 1950) perche' e' vero che la georeferenza e' WGS84, ma il reticolo "disegnato" sulla cartina e' sempre ED50!!! Con il Datum giusto, dovreste vedere i due reticoli (blu' e originale) assolutamente sovrapposti.

7) avete finito! Prendete il vostro fantastico nuovo file .ecw e caricatelo sul vostro GPS cartografico! Buon cammino!!:cool::si:
 
Ultima modifica di un moderatore:
Come diceva un mio amico di infanzia, a proposito di qualsiasi cosa trovassimo e che lui si metteva in tasca... "Può far comodo" :D

Grazie delle info che ho provveduto a salvare opportunamente.
 
Molto utile ed interessante. Grazie Sun.

Esiste un equivalente per mac di FWtools? Vedo che nel link che hai postato c'è la versione per linux. Dovrei provare ad usare quella.
 
S

sun75

Guest
Ultima modifica di un moderatore:
Che poi non vorrei dir fesserie ma il gis che ho installato sul mio mac mi sa che le ritaglia già georefernziate le carte... stasera controllo.
 
S

sun75

Guest
Che poi non vorrei dir fesserie ma il gis che ho installato sul mio mac mi sa che le ritaglia già georefernziate le carte... stasera controllo.

Se usi Grass, e' possibile!!:roll:

Me lo sto giusto scaricando: se riesco a fare tutto quello che devo fare su Grass, che e' open-source e gira ovunque, ho un paio di progettini :) su alcuni miei server Linux....;)
 
Ultima modifica di un moderatore:
Ho detto una fesseria. Con il programma che ho sul mac, Macgpspro, non posso fare proprio niente. Devo utilizzare la procedura descritta da te. Grass provai a scaricarlo una volta, solo che è pesantissimo da quello che ricordo.
 
S

sun75

Guest
Ho detto una fesseria. Con il programma che ho sul mac, Macgpspro, non posso fare proprio niente. Devo utilizzare la procedura descritta da te.

Che poi, alla fine, sembra complicata, ma sono solo due comandi effettivi...una volta presa "la mano"...e' roba di 5 minuti!;)

Grass provai a scaricarlo una volta, solo che è pesantissimo da quello che ricordo.

Si, vedo....e mi sembra anche alquanto ostico, tanto che una coppia di scrittori ci ha scritto un libro di 450pg!!:oops:
 
Me lo sto giusto scaricando: se riesco a fare tutto quello che devo fare su Grass, che e' open-source e gira ovunque, ho un paio di progettini :) su alcuni miei server Linux....;)
anch'io ho grass, installato su ubuntu che uso come sistema operativo di default

a dire la verità non mi sembra pesante...anche se ammetto di non aver mai provato ad impararlo :roll:

daniele, stasera provo a seguire la procedura che hai indicato ;)

p.s. la procedura che ho fatto e che ti ho descritto in quell'altro topic era corretta? se sì, perchè facendo in quel modo si perde in precisione?
 
S

sun75

Guest
p.s. la procedura che ho fatto e che ti ho descritto in quell'altro topic era corretta? se sì, perchè facendo in quel modo si perde in precisione?

Te l'ho gia' scritto di la': tu hai georeferenziato una mappa usando come datum ED50, prendendo le coordinate dal PCN, che usa come datum WGS84!!

Dovevi georeferenziare la mappa usando come datum WGS84, quindi una volta refereziata e visualizzata, cambiare datum di visualizzazione in ED50 (per allineare i reticoli).

Ciao ciao!;)
 
Dovevi georeferenziare la mappa usando come datum WGS84, quindi una volta refereziata e visualizzata, cambiare datum di visualizzazione in ED50 (per allineare i reticoli)

ok, tutto chiaro tranne l'ultimo punto, cioè come cambiare successivamente il datum per allineare i reticoli

ecco ora la mia situazione:
 

Allegati

  • Schermata-windows xp [In esecuzione] - Oracle VM VirtualBox.jpg
    Schermata-windows xp [In esecuzione] - Oracle VM VirtualBox.jpg
    94,4 KB · Visite: 2.290
Ultima modifica:
ok...ora i reticoli coincidono

immaginiamo ora di voler fornire le coordinate chilometriche di p.ta di paipo

devo darle in questo formato : UTM{zona}{fascia}{coord. est}E{coord. nord}N

giusto?

ed è ai valori segnati ora sulle linee del reticolo che devo riferirmi?

il quadrato in cui si trova il punto segnato sarebbe dato (in senso orario e dal vertice alto a sinistra) da [(4498000,459000), (4498000,460000), (4497000,459000), (4497000, 460000)]

ovviamente dovrei usare le opportune proporzioni
 

Allegati

  • praiano12.jpg
    praiano12.jpg
    92,2 KB · Visite: 2.326
S

sun75

Guest
ok...ora i reticoli coincidono

immaginiamo ora di voler fornire le coordinate chilometriche di p.ta di paipo

devo darle in questo formato : UTM{zona}{fascia}{coord. est}E{coord. nord}N

giusto?

ed è ai valori segnati ora sulle linee del reticolo che devo riferirmi?

il quadrato in cui si trova il punto segnato sarebbe dato (in senso orario e dal vertice alto a sinistra) da [(4498000,459000), (4498000,460000), (4497000,459000), (4497000, 460000)]

ovviamente dovrei usare le opportune proporzioni

Quando dai delle coord. le devi dare nel formato scritto aggiungendo sempre il datum utilizzato.

Quindi il punto da te richiesto, a occhio, sara'

UTM33T459600E4497400N - ED50

Le coppie di coord. vanno sempre date scrivendo prima la Est (asse X) e poi la Nord (asse Y)!!

Considerando che ogni quadrato ha il lato di 1Km e che le coord. non sono altro che la misura in metri su quella proiezione, viene da se, che se tu hai le coord. dell'angolo inferiore sx rispetto al punto che ti interessa, per avere le coord. del punto devi sommarci la distanza x,y in metri rispetto al punto "angolo" di cui sopra. Chiaro?:D

Ciao!;)
 
perfettamente ;)

ultima domanda: visto che a me non servono mappe per gps preferisco seguire (finchè non comprerò un gps) la prima parte della procedura da te descritta, cioè fin quando si arriva a visualizzare il reticolo con le etichette e farlo coincidere con quello della mappa igm, poi mi esporto il pezzo di mappa in jpg e me lo stampo su un A4

ora però vorrei chiederti: come faccio a mantenere la scala 1:25000 quando andrò a stampare?
 
S

sun75

Guest
perfettamente ;)

ultima domanda: visto che a me non servono mappe per gps preferisco seguire (finchè non comprerò un gps) la prima parte della procedura da te descritta, cioè fin quando si arriva a visualizzare il reticolo con le etichette e farlo coincidere con quello della mappa igm, poi mi esporto il pezzo di mappa in jpg e me lo stampo su un A4

ora però vorrei chiederti: come faccio a mantenere la scala 1:25000 quando andrò a stampare?

Puoi stampare anche direttamente da OziExplorer: ha una "fantastica" funzione di "stampa in scala" e tu selezioni il rapporto, ovvero 1:1 in questo caso (infatti la scala e' riferita a cio' che visualizzi e non "alla realta'"!). Prepara un plotter A0....pero'!:D

Difficilmente riuscirai a stampare nella scala 1:1 della mappa, a meno che non siano porzioni veramente piccole! Ti conviene stampare come viene...e poi le coordinate le prendi "a occhio", dividendo sempre "x 10" ogni quadrato del reticolo...cosi' da avere almeno le centinaia di metri, che e' gia' molto!:cool:
 
Alto Basso