rigacci.org

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== 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:

<code>
Band 1 Block=2160x1 Type=Float32, ColorInterp=Gray
  NoData Value=-9999
</code>

<code>
Band 1 Block=2160x1 Type=Int16, ColorInterp=Gray
  NoData Value=-9999
</code>

<code>
Band 1 Block=64x64 Type=Byte, ColorInterp=Gray
</code>

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:

<code>
gdal_translate -of GTiff -ot Float32 minrainfal.asc minrainfal.tif
gdal_translate -of GTiff -ot Int16   maxrainfal.asc maxrainfal.tif
</code>

===== 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:

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

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:

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

===== 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:

<code>
gdal_translate -co TILED=YES -co TFW=YES -co COMPRESS=NONE src.tif dst.tif
</code>

Con **''gdalinfo''** si potrà vedere il blocksize:

<code>
Band 1 Block=256x256 Type=Byte, ColorInterp=Palette
</code>

===== 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:
<code>
gdalinfo TrueMarble.2km.21600x10800.tif
...
Origin = (-180.000000000000000,90.000000000000000)
Pixel Size = (0.016666666666667,-0.016666666666667)
...
</code>

Un trucco per creare il world file a partire da un GeoTiff usando i tool gdal è il seguente:

<code>
gdal_translate -co "TFW=YES" TrueMarble.2km.21600x10800.tif temp.tif
mv temp.tfw TrueMarble.2km.21600x10800.tfw
rm temp.tif
</code>

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:

<file>
0.0166666667
0.0000000000
0.0000000000
-0.0166666667
-179.9916666667
89.9916666667
</file>

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

<code>
convert world_topo_bathy.png -compress none world_topo_bathy.tif
</code>

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:

<code>
gdal_translate -a_srs EPSG:4326 world_topo_bathy.png world_topo_bathy.tif
</code>

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:

<code>
gdaladdo TrueMarble.8km.5400x2700.tif 2 4 8
</code>

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

<code>
gdaladdo -ro -r gauss TrueMarble.8km.5400x2700.tif 2 4 8
</code>

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''**:

<code>
gdaladdo --config USE_RRD YES TrueMarble.8km.5400x2700.tif 2 4 8
</code>

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:

<code>
gdal_merge.py -o output.tif -of GTiff -co TILED=YES -co TFW=YES -ps 0.9765625 0.9765625 2001/*.jpg
</code>

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:

<code>
(670510 - 654510) / 16384 = 0.9765625
</code>