====== Raster con GDAL e MapServer ======
===== Raster type: Float32, Int16, Byte =====
La gestione dei **raster Tiff** con GDAL e MapServer non funziona sempre alla perfezione come si vorrebbe. Si deve fare sempre attenzione al tipo di raster con cui si ha a che fare, potrebbe trattarsi di **Float32**, **Int16** oppure **Byte**.
Con **''gdalinfo''** è possibile ottenere l'informazione, ecco alcuni esempi di output su raster costituiti da una singola banda:
Band 1 Block=2160x1 Type=Float32, ColorInterp=Gray
NoData Value=-9999
Band 1 Block=2160x1 Type=Int16, ColorInterp=Gray
NoData Value=-9999
Band 1 Block=64x64 Type=Byte, ColorInterp=Gray
La conversione con **''gdal_translate''** - ad esempio dal formato **Arc/Info ASCII Grid** - potrebbe non indovinare il formato corretto, con il risultato di troncare alcuni valori.
Ecco come forzare l'utilizzo del tipo adeguato:
gdal_translate -of GTiff -ot Float32 minrainfal.asc minrainfal.tif
gdal_translate -of GTiff -ot Int16 maxrainfal.asc maxrainfal.tif
===== Classificazione di raster con MapServer =====
Per ottimizzare la classificazione di un raster **Float32** o **Int16** con MapServer, conviene aggiungere le direttive per la gestione dei raster tramite **''LAYER.PROCESSING''**, in particolare **''SCALE''** e **''SCALE_BUCKETS''**. Ecco un esempio in cui i valori compresi tra -55.0 e 45.0 vengono classificati in 24 classi:
LAYER
NAME "temperature"
PROCESSING "SCALE=-55.00,45.00"
PROCESSING "SCALE_BUCKETS=24"
CLASS
NAME "< -50"
EXPRESSION ([pixel] < -50.83)
STYLE
COLOR 0 53 255
END
END
CLASS
NAME "-50 ÷ -46"
EXPRESSION ([pixel] >= -50.83 and [pixel] < -46.67)
STYLE
COLOR 0 106 255
END
END
...
Vedere in proposito [[http://www.mapserver.org/input/raster.html#raster|MapServer Raster Data]]. Tali direttive non servono (anzi **creano problemi**) con raster a 8 bit.
Per verificare l'elaborazione del raster si può utilizzare **''shp2img''**, in questo modo:
shp2img -m /path/to/file.map -all_debug 10 -o /path/to/raster
...
msDrawRasterGDAL_16BitClassification(GEBCO):
scaling to 8913 buckets from range=-5600.5,3312.5.
...
===== Ottimizzare raster con gdal_translate =====
Alcune caratteristiche di un file TIFF possono incidere sulle performance di MapServer. Su file molto grossi ad esempio la **memorizzazione in blocchi** (256x256) piuttosto che in righe di pixel può migliorare i tempi di accesso. Anche la **compressione** può incidere negativamente sulle prestazioni.
Ecco un comando che converte un file TIFF in un'altro file TIFF, l'organizzazione sarà a blocchi di 256x256 pixel e nessuna compressione. Le informazioni geografiche saranno memorizzate in un world file (.tfw) esterno:
gdal_translate -co TILED=YES -co TFW=YES -co COMPRESS=NONE src.tif dst.tif
Con **''gdalinfo''** si potrà vedere il blocksize:
Band 1 Block=256x256 Type=Byte, ColorInterp=Palette
===== Georiferire un file PNG =====
Quei mattacchioni della NASA fanno [[http://visibleearth.nasa.gov/view_set.php?categoryID=2355|scaricare il BlueMarble]] in formato PNG, non georiferito!
Bisognerebbe creare un //world file// secondo la struttura documentata da [[http://mapserver.org/input/raster.html#georeference-with-world-files|Georeference with World Files]]. Il problema è sapere i parametri! Per fortuna abbiamo un GeoTiff del TrueMarble alla stessa risoluzione e nello stesso sistema di riferimento, con **''gdalinfo''** vediamo i parametri del GeoTiff:
gdalinfo TrueMarble.2km.21600x10800.tif
...
Origin = (-180.000000000000000,90.000000000000000)
Pixel Size = (0.016666666666667,-0.016666666666667)
...
Un trucco per creare il world file a partire da un GeoTiff usando i tool gdal è il seguente:
gdal_translate -co "TFW=YES" TrueMarble.2km.21600x10800.tif temp.tif
mv temp.tfw TrueMarble.2km.21600x10800.tfw
rm temp.tif
Creiamo il world file con lo stesso nome del file ''.png'', ma estensione **''.wld''** oppure **''.tfw''** (un world file con tale estensione viene trovato sia da MapServer che dalla libreria GDAL). Le coordinate del centro del pixel in alto a sinistra sono le cordinate dell'origine riportate da ''gdalinfo'', meno la metà di un pixel:
0.0166666667
0.0000000000
0.0000000000
-0.0166666667
-179.9916666667
89.9916666667
Peccato che la velocità di accesso non sia accettabile, il manuale di GDAL recita: //PNG files are linearly compressed, so random reading of large PNG files can be very inefficient//.
Una **soluzione non ottimale** è quella di convertire in Tiff non compresso con ImageMagick (il world file da usare rimane esattamente lo stesso):
convert world_topo_bathy.png -compress none world_topo_bathy.tif
La **soluzione migliore** è creare un **GeoTiff completo**, che non ha bisogno del world file perché incorpora tutte le informazioni geografiche necessarie. A tale scopo si può usare il comando:
gdal_translate -a_srs EPSG:4326 world_topo_bathy.png world_topo_bathy.tif
Nel GeoTiff risultante vengono incorporate sie le informazioni relative al sistema di riferimento che quelle contenute nel world file.
===== Overview Images (piramidi) =====
Se il file è troppo grande conviene creare le overview, dette anche //piramidi//, altrimenti uno zoom alla massima estensione obbliga a leggere e ridimensionare al volo tutto il raster.
In genere si generano delle copie dell'immagine con risoluzione 1/2, 1/4, 1/8, ... dell'originale. Le overview possono essere contenute all'interno dello stesso GeoTiff:
gdaladdo TrueMarble.8km.5400x2700.tif 2 4 8
Se invece vogliamo tenere le overview in un file separato (che avrà estensione **''.tif.ovr''**) ed utilizzare un algoritmo di ricampionamento migliore (//gauss//, ma anche //average// è migliore del //nearest// predefinito):
gdaladdo -ro -r gauss TrueMarble.8km.5400x2700.tif 2 4 8
Con le versioni di GDAL < 1.6.0 si deve usare una sintassi diversa, il file esterno con le overview avrà formato //Erdas Image// ed estensione **''.aux''**:
gdaladdo --config USE_RRD YES TrueMarble.8km.5400x2700.tif 2 4 8
MapServer e la libreria Gdal accedono automaticamente alle overview anche se si trovano nel file separato **''.aux''**.
===== Creare una overview di una copertura tiled =====
Anche se nelle tile sono state aggiunte le overview, ai fattori di scala maggiori (zoom minore) è necessario accedere a tutti i file (tile) della copertura per generare la mappa. Questo fa decadere di molto le prestazioni, è quindi opportuno creare una **versione ridotta della copertura**, contenuta in un **singolo file**.
Tale file conviene che sia TIFF memorizzato a blocchi ed eventualmente deve comprendere le overview opportune. Ecco come usare **''gdal_merge''** per ottenere il risultato:
gdal_merge.py -o output.tif -of GTiff -co TILED=YES -co TFW=YES -ps 0.9765625 0.9765625 2001/*.jpg
La dimensione del pixel (in metri) viene calcolata dividento la larghezza della copertura per il numero di pixel desiderati. Ad esempio se si ha una superficie larga 16 km e si vuole creare un'immagine larga 16384 pixel:
(670510 - 654510) / 16384 = 0.9765625