Virtuaalirasterit

Virtuaalirasterit ovat hyödyllinen GDAL-konsepti suurten rasteriaineistojen hallintaan, kun aineisto on jaettu toisiaan leikkaamattomiin karttalehtiin. Virtuaalirasterit eivät sen sijaan sovellu esimerkiksi aikasarjojen tai päällekkäisten rasterien, kuten kaukokartoitusruutujen, hallintaan.

Teknisesti virtuaalirasteri on vain pieni xml-tiedosto, joka kertoo GDALille, missä varsinaiset datatiedostot sijaitsevat, mutta käyttäjän näkökulmasta virtuaalirastereita voidaan käsitellä pitkälti kuten mitä tahansa muuta rasterimuotoa. Virtuaalirasterit voivat sisältää rasteridataa missä tahansa tiedostomuodossa, jota GDAL tukee. Virtuaalirasterit ovat hyödyllisiä, koska niiden avulla suuria aineistoja voidaan käsitellä ikään kuin ne olisivat yksi tiedosto, jolloin oikeiden tiedostojen etsimistä ei tarvita.

Esimerkiksi Maanmittauslaitoksen 2 m:n ja 10 m:n korkeusmallit ovat saatavilla Puhdissa. Nämä aineistot on jaettu useisiin tif-tiedostoihin (karttalehtiin), ja jos haluaisimme esimerkiksi laskea vyöhyketilastoja alueille, jotka sijaitsevat eri puolilla Suomea, meidän pitäisi jotenkin selvittää, mikä tiedosto kattaa minkäkin alueen, ja laskea tilastot oikeasta tiedostosta. Lisähaasteita syntyisi, jos alue, jolle haluamme laskea tilastoja, sattuu sijaitsemaan kahden tai useamman karttalehden rajalla. Samankaltaisia reunaefekteihin liittyviä ongelmia syntyisi esimerkiksi käytettäessä fokustoimintoja, joissa tarvitaan myös ympäröivien tiedostojen tietoja. Nämä ongelmat voidaan helposti välttää luomalla koko tutkimusalueelle virtuaalirasteri, jolloin GDAL huolehtii automaattisesti edellä mainituista asioista.

Virtuaalirastereita on mahdollista käyttää myös niin, että vain pieni xml-tiedosto tallennetaan paikallisesti ja suuret rasteritiedostot sijaitsevat Altaassa, Amazon S3:ssa, julkisesti palvelimella tai missä tahansa muussa paikassa, jota GDALin virtuaaliohjaimet tukevat. Data siirretään paikallisesti vain siltä alueelta ja siltä zoomaustasolta, jota pyydetään, kun virtuaalirasteri avataan. Parhaiten toimiva tiedostomuoto rasteridatan tallentamiseen etäpalveluun on Cloud optimized GeoTIFF, mutta myös muut muodot ovat mahdollisia.

Virtuaalirasterien käyttö

Virtuaalirastereita tukevat kaikki GDAL-pohjaiset työkalut, mukaan lukien Pythonin ja R:n paikkatietopaketit, ArcGIS, FME, GrassGIS, MapInfo, QGIS ja SagaGIS. Lisäksi ArcGISillä ja MapInfolla on kummallakin myös oma virtuaalirasterimuotonsa, joka muistuttaa GDALin virtuaalirasteria.

Vain osan suuren virtuaalirasterin lukeminen bboxin avulla

Virtuaalirasterien tehokkaan käytön avain on usein tuki sille, että suuresta virtuaalirasterista luetaan vain osa. Tämä onnistuu sekä R:llä että Pythonilla; monissa muissa työkaluissa rajaus täytyy tehdä erillisenä vaiheena ja rajattu data tallentaa tiedostoon.

GDALilla on helppo rajata pieni osa suuresta virtuaalirasterista:

gdal_translate -projwin 614500 6668500 644500 6640500 test.vrt test_clip.tif

R:

Terra:

library(terra)  
vrt <- rast("test.vrt")  
data = crop(vrt , ext(614500, 644500, 6640500, 6668500))

Raster:

library(raster)  
vrt <- raster("test.vrt")  
data = crop(vrt , extent(614500, 644500, 6640500, 6668500))

Python:

import rasterio
from rasterio.windows import from_bounds
vrt_path = "test.vrt"
with rasterio.open(vrt_path) as src:
    rst = src.read(window=from_bounds(614500, 6640500, 644500, 6668500, src.transform))

Analyysin suorittaminen siten, että vain joitakin osia virtuaalirasterista luetaan

On mahdollista työskennellä erittäin suurten virtuaalirasterien kanssa silloin, kun analyysin ei tarvitse tuottaa samankokoista rasteritulosta. Hyvä esimerkki tästä on vyöhyketilastojen laskeminen polygonialueille, jotka sijaitsevat laajalla alueella, katso esimerkiksi CSC training github. Tässä tapauksessa rasteridataa luetaan vain niiltä alueilta, jotka menevät päällekkäin polygonien kanssa.

Suurten virtuaalirasterien visuaalinen käyttö

On hyvä huomata, että vaikka jonkin analyysin suorittaminen koko Suomen kattavalla 2 m:n korkeusmallilla on täysin mahdollista Puhdissa perus-.vrt-tiedoston avulla, datan tarkastelu esimerkiksi QGISillä ei ole käytännöllistä näin suurelle aineistolle ilman lisäoptimointia. Jos haluat tarkastella suurta virtuaalirasteria helposti, sinun täytyy tehdä muutama asia:

Luo yleiskatselutasot virtuaalirasterillesi gdaladdo-komennolla. Kannattaa huolehtia siitä, ettet luo niin suuria yleiskatselutasoja, että niistä itsestään tulee valtavan suuri tiedosto.
Jos virtuaalirasterisi on todella suuri, on järkevää luoda hierarkkinen virtuaalirasterirakenne, jossa ylimmän tason virtuaalirasteri viittaa pienempiin virtuaalirastereihin, jotka puolestaan viittaavat vielä pienempiin virtuaalirastereihin ja niin edelleen, kunnes viimeinen virtuaalirasteri viittaa varsinaisiin tiedostoihin. Tämän lähestymistavan syy on se, että muuten myös käytettävät yleiskatselutasot kasvavat todella suuriksi. Huomaa, että tällaisen hierarkkisen rakenteen käyttö voi aiheuttaa joitakin artefakteja analyysiä ajettaessa, joten sitä kannattaa käyttää vain visualisointitarkoituksiin.
Laske tilastot valmiiksi virtuaalirastereillesi ja lähdetiedostoillesi. Tämä nopeuttaa tiedostojen avaamista esimerkiksi QGISissä. QGISin täytyy ottaa näytteitä datan minimi- ja maksimiarvoista, jotta se voi asettaa väriskaalan oikein, ja tämä vie aikaa suurilla virtuaalirastereilla. Tämän välttämiseksi voit laskea tilastot etukäteen erilliseen XML-tiedostoon komennolla gdalinfo --stats.
Hyvä niksi QGISissä suurten rasterien kanssa työskenneltäessä on ottaa rasterityökalupalkki käyttöön (View->Toolbars->Raster Toolbar). Sen avulla voit helposti säätää väriskaalan näytöllä näkyvän alueen mukaan, jolloin saat hyvän kontrastin zoomaustasosta riippumatta.
QGIS näyttää käsittelevän suuria aineistoja varsin hyvin, kun edellä mainitut vaiheet on tehty. Jopa koko Suomen 2 m:n korkeusmallilla zoomaus ja kartan siirtäminen ovat melko sujuvia.

Virtuaalirasterien luominen

Seuraavat työkalut tukevat virtuaalirasterien luomista:

GDAL gdalbuildvrt komentorivityökalu.
GDAL gdal_translate STAC-haulla
Pythonissa ja R:ssä on kääreet GDALin gdalbuildvrtille, ks. pidempi esimerkki R:lle StackOverflow-vastauksessa.
QGIS, GrassGIS ja SagaGIS tarjoavat graafisen käyttöliittymän gdalbuildvrtille
lidR tukee lidar-data-analyysin tulosten kirjoittamista suoraan virtuaalirasteriksi
vrt_creator.py Puhdissa mukautetuille alueille 2 m:n tai 10 m:n korkeusmallista

Puhdissa glalbuildvrt sisältyy kaikkiin GDALin sisältäviin moduuleihin, Pythonin BuildVRT geocondaan, QGIS QGISiin sekä R:n terra ja lidR r-env -moduuliin.

Virtuaalirasterin luominen GDALin gdalbuildvrtillä

Jos datasi on jaettu alikansioihin, helpoin vaihtoehto on luoda virtuaalirasteri syötetiedostoluettelon avulla:

gdalbuildvrt -input_file_list file_list.txt virtual_raster.vrt

Katso gdalbuildvrtin dokumentaatiosta lisävalintoja, esimerkiksi no-data-arvon asettamiseen.

Tiedostoluettelon luominen gdalbuildvrtille

Tiedostoluettelon tulisi sisältää mieluiten täydet polut, mutta paikallisille tiedostoille voidaan käyttää myös suhteellisia polkuja.

Linux

find /appl/data/geo/mml/dem2m/ -name "*.**tif**" > file_list.txt

Windows

dir \data\dem2m\*.**tif** /S /B > file_list.txt

Rasteritiedostot Altaassa / jossakin muussa S3-palvelussa

Jos teet tämän Puhdissa, lataa allas-moduuli.

Listaa tiedostonimet siinä muodossa kuin ne ovat ämpärissä käyttäen rclonea tai jotain muuta työkalua:

rclone lsf --include '*.**tif**' allas:<your_bucket_name > file_list.txt

Lisää seuraavaksi tiedostoluetteloon täydet polut siinä muodossa kuin GDAL niitä vaatii käyttäen vsicurl-, vsis3- tai vsiswift-ajureita. Katso pidemmät selitykset GDAL-ajureista ja Altaasta Puhdin GDAL-sivulta.

sed -i -e 's-^-/**vsicurl**/https://a3s.fi/<your_bucket_name>/-' file_list.txt

sed -i -e 's-^-/**vsis3**/<your_bucket_name>/-' file_list.txt

sed -i -e 's-^-/**vsiswift**/<your_bucket_name>/-' file_list.txt

Määritä tunnistetietosi GDALille ennen gdalbuildvrt-komennon suorittamista.

Virtuaalirasterin luominen STAC-haun tuloksista

STAC (Spatio-Temporal Asset Catalog) on tapa kuvata (rasteri)aineistoja siten, että dataa voidaan hakea ajan ja sijainnin perusteella. Esimerkiksi Paituli STAC sisältää useita suomalaisia aineistoja ja selittää STAC-käsitteitä tarkemmin.

GDAL tukee STAC-hakua ja virtuaalirasterien luomista hakutulosten perusteella.

Esimerkiksi seuraava komento luo virtuaalirasterin Paituli STACin hakutulosten perusteella. Kysely pyytää tiedostoja corine_land_cover_at_geocubes-kokoelmasta määritellyllä aikavälillä ja sijainnilla (bbox) käyttäen Assetteja (=tiedostoja), joiden nimi on COG.

Linuxissa ja Macissa:

gdal_translate "STACIT:\"https://paituli.csc.fi/geoserver/ogc/stac/v1/search?collections=corine_land_cover_at_geocubes&datetime=2017-01-05/2019-02-14&bbox=19.5,61.5,28.7,63.0\":asset=COG" -of VRT corine.vrt

Windowsissa (lisää ^ ennen &-merkkejä):

gdal_translate "STACIT:\"https://paituli.csc.fi/geoserver/ogc/stac/v1/search?collections=corine_land_cover_at_geocubes^&datetime=2017-01-05/2019-02-14^&bbox=19.5,61.5,28.7,63.0\":asset=COG" -of VRT corine.vrt

Huomaa, että GDAL poistaa virtuaalirasterista Assetit, jotka uudemmat Assetit peittävät kokonaan. Esimerkiksi Paituli STACissa CORINE-dataa on saatavilla vuosilta 2012 ja 2018. Jos yllä olevaa hakua muutetaan kattamaan myös vuosi 2012, STAC-haku löytäisi myös CORINE 2012 -Itemit, mutta niitä ei silti olisi luodussa VRT:ssä. Jos vuoden 2012 data tarvitaan virtuaalirasteriin, aikaväliä täytyy säätää niin, ettei vuosi 2018 sisälly siihen.