Dati da un server WMS: scaricarli, riproiettarli, comprimerli, tassellarli e piramidarli da riga di comando
Come faremmo senza GDAL?
Recentemente in OpenDataSicilia è nata l’esigenza di fare il download di una sorgente di dati WMS e produrre in output un file da archiviare. Uno dei modi per farlo è usare la straordinaria GDAL.
Nella fattispecie era utile avere in output un file GeoTIFF che fosse:
- riproiettato secondo uno Spatial Reference System diverso da quello originale;
- ben compresso, al limite anche con un sistema lossy;
- tassellato e piramidato.
Non scriverò (quasi) nulla sui temi di fondo e lascio al lettore eventuali approfondimenti. A fine post ho inserito delle letture consigliate.
GDAL
È la più importante libreria per la gestione di dati cartografici. L’ho usata perché molto semplice, efficiente e perché la conosco un po’. È open source e utilizzabile su qualsiasi sistema operativo: qui i binari e qui il codice sorgente.
Ho usato alcune tra le sue numerose utility a riga di comando.
Ottenere informazioni sul server WMS
A partire dall’URL di unserver WMS è possibile fare una richiesta e avere le informazioni di base sui dati esposti. L’utility di GDAL da usare è gdalinfo con una sintassi che (con la source in questione) sarà:
gdalinfo wms:"http://mapsrv.comune.palermo.it/erdas-iws/ogc/wms/CARTE_STORICHE?"
Nella risposta si ottengono diverse informazioni utili, tra cui quelle sui layer esposti. Uno di questi, quello che verrà usato per questo tutorial, è il subdataset 5.
SUBDATASET_5_NAME=WMS:http://mapsrv.comune.palermo.it/erdas-iws/ogc/wms/CARTE_STORICHE?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=IRTA_1956_2000&SRS=EPSG:3004&BBOX=2368600,4215099.9435559697,2383000,4230700
SUBDATASET_5_DESC=IRTA_1956_2000Ottenere informazioni su un determinato layer
L’utility da utilizzare è sempre gdalinfo, avendo stavolta cura di specificare con il parametro -sdil numero del subdataset:
gdalinfo wms:"http://mapsrv.comune.palermo.it/erdas-iws/ogc/wms/CARTE_STORICHE?" -sd 5Tra le informazioni importanti che verranno restituite quelle sulla proiezione:
PROJCS[“Monte Mario / Italy zone 2”,
GEOGCS[“Monte Mario”,
DATUM[“Monte_Mario”,
SPHEROID[“International 1924”,6378388,297,
AUTHORITY[“EPSG”,”7022"]],
TOWGS84[-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68],
AUTHORITY[“EPSG”,”6265"]],
PRIMEM[“Greenwich”,0,
AUTHORITY[“EPSG”,”8901"]],
UNIT[“degree”,0.0174532925199433,
AUTHORITY[“EPSG”,”9122"]],
AUTHORITY[“EPSG”,”4265"]],
PROJECTION[“Transverse_Mercator”],
PARAMETER[“latitude_of_origin”,0],
PARAMETER[“central_meridian”,15],
PARAMETER[“scale_factor”,0.9996],
PARAMETER[“false_easting”,2520000],
PARAMETER[“false_northing”,0],
UNIT[“metre”,1,
AUTHORITY[“EPSG”,”9001"]],
AXIS[“X”,EAST],
AXIS[“Y”,NORTH],
AUTHORITY[“EPSG”,”3004"]]E anche le coordinate dei vertici del layer, espresse in latitude/longitude e nel sistema di coordinate nativo:
Corner Coordinates:
Upper Left ( 2368600.000, 4230700.000) ( 13d16'14.70"E, 38d12'39.27"N)
Lower Left ( 2368600.000, 4215099.944) ( 13d16'26.63"E, 38d 4'13.36"N)
Upper Right ( 2383000.000, 4230700.000) ( 13d26' 6.58"E, 38d12'47.58"N)
Lower Right ( 2383000.000, 4215099.944) ( 13d26'17.38"E, 38d 4'21.62"N)
Center ( 2375800.000, 4222899.972) ( 13d21'16.32"E, 38d 8'30.56"N)Creazione di un file sorgente XML
Per scaricare questo layer dal server WMS — secondo la procedura descritta in seguito — è necessario creare un file XML, secondo le specifiche ufficiali di GDAL e che deve contenere diverse informazioni sulla sorgente, tra cui:
- l’URL del server;
- lo Spatial Reference System del layer sorgente;
- il nome del layer;
- le coordinate del vertice in alto a sinistra e di quello in basso a destra, secondo lo Spatial Reference System definito prima;
- la dimensione in pixel del file di output, che dipenderà dalla risoluzione desiderata. Se ad esempio volessi una risoluzione di 0.5 metri, per un’area che copre 2000 metri lungo l’asse delle ascisse, la larghezza dell’immagine di output dovrà essere di 4000 pixel.
<GDAL_WMS>
<Service name="WMS">
<Version>1.1.1</Version>
<ServerUrl>http://mapsrv.comune.palermo.it/erdas-iws/ogc/wms/CARTE_STORICHE?</ServerUrl>
<SRS>EPSG:3004</SRS>
<ImageFormat>image/png</ImageFormat>
<Layers>IRTA_1956_5000</Layers>
<Styles></Styles>
</Service>
<DataWindow>
<UpperLeftX>2366800</UpperLeftX>
<UpperLeftY>4231900</UpperLeftY>
<LowerRightX>2382999</LowerRightX>
<LowerRightY>4213900</LowerRightY>
<SizeX>32398</SizeX>
<SizeY>36000</SizeY>
</DataWindow>
<Timeout>3000</Timeout>
<MaxConnections>1</MaxConnections>
<Projection>EPSG:3004</Projection>
</GDAL_WMS>Qui un esempio di file XML.
Download
Costruito il file XML, il download sarà eseguito con gdal_translate:
gdal_translate IRTA_1956_5000.xml IRTA_1956_2000.tifFatevi un caffè, e aspettate che arrivi sino in fondo.
Riproiezione
Era utile riproiettare il file scaricato dal sistema di coordinate originale — EPSG:3004 — al sistema EPSG:32633. È una trasformazione di coordinate che va fatta tenendo conto dei parametri Bursa-Wolf adeguati per quest’area.
L’utility per la riproiezione è gdalwarp:
gdalwarp -s_srs "+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9996 +x_0=2520000 +y_0=0 +ellps=intl +units=m +no_defs +towgs84=-50.2,-50.4,84.8,-0.690,-2.012,0.459,-28.08" -t_srs "+proj=utm +zone=33 +ellps=WGS84 +units=m +datum=WGS84 +no_defs" IRTA_1956_2000.tif IRTA_1956_2000_32633.tifCompressione
L’immagine di output è un file di grandi dimensioni e si può comprimere in varie modalità. Un buon compromesso tra qualità di output e rapporto di compressione è quello descritto in questo post, che si realizza sempre con gdal_translate:
gdal_translate -co "TILED=YES" -co "INTERLEAVE=PIXEL" -co "COMPRESS=JPEG" -co "PHOTOMETRIC=YCBCR" -co "JPEG_QUALITY=80" -a_srs "EPSG:3004" IRTA_1956_2000_32633.tif IRTA_1956_2000_32633_cp.tifIl comando di sopra oltre a comprimere l’immagine, la suddivide in tasselli e ne velocizza (per i client compatibili) la velocità di rendering a monitor.
Piramidazione
Per rendere l’immagine ancora più leggibile in termini di rapidità di operazioni di pan & zoom è consigliato anche piramidarla:
gdaladdo -r cubic --config COMPRESS_OVERVIEW JPEG --config INTERLEAVE_OVERVIEW PIXEL IRTA_1956_2000_32633_cp.tif 2 4 8 16L’output
L’output di questa procedura è finito in questa mappa di Giovan Battista Vitrano, che consente di osservare come è cambiata Palermo dagli anni 30 ad oggi, confrontando cartografie di periodi diversi.
Come faremmo senza GDAL?
È una libreria, uno strumento, di grandissima utilità, la cui conoscenza di base è a mio avviso propedeutica per chiunque voglia “fare cose” bene con i dati spaziali.
Trovate online decine di tutorial, numerosissimi thread in forum e mailing list, e anche dei libri (come questo, che consiglio).
Alcune riferimenti utili (citati anche nel testo):
