Virtuaalirasterit

Virtuaalirasterit ovat hyödyllinen GDAL-konsepti suurten rasteriaineistojen hallintaan, kun aineisto on jaettu päällekkäisyyttä sisältämättömiin 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 GDALin tukemassa tiedostomuodossa. Virtuaalirasterit ovat hyödyllisiä, koska niiden avulla suuria aineistoja voidaan käsitellä ikään kuin ne olisivat yksi tiedosto, jolloin oikeiden tiedostojen etsiminen ei ole tarpeen.

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 GDALin virtuaaliohjainten tukemassa paikassa. Data siirretään paikallisesti vain siltä alueelta ja sillä zoomaustasolla, 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 myös omat GDALin virtuaalirasteria muistuttavat virtuaalirasterimuotonsa.

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 tallentaa rajattu data 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 niin, että vain osa 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, ks. 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 esimerkiksi koko Suomen kattavan 2 m:n korkeusmallin analysointi perus-.vrt-tiedostolla on täysin mahdollista Puhdissa, datan tarkastelu esimerkiksi QGISillä ei ole käytännöllistä näin suurelle aineistolle ilman lisäoptimointia. Jos haluat tarkastella suurta virtuaalirasteria sujuvasti, sinun täytyy tehdä muutama asia:

Luo yleiskatselutasot virtuaalirasterillesi gdaladdo-komennolla. Kannattaa huolehtia siitä, etteivät yleiskatselutasot ole niin suuria, että niistä itsestään tulee valtava 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ävistä yleiskatselutasoista tulee todella suuria. 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 otos 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). Tämän avulla voit helposti säätää väriskaalaa näytöllä näkyvän alueen mukaan, jolloin kontrasti pysyy hyvänä 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-analyysitulosten 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 data on jaettu alikansioihin, helpoin vaihtoehto on luoda virtuaalirasteri syötetiedostoluettelon avulla:

gdalbuildvrt -input_file_list file_list.txt virtual_raster.vrt

Katso lisäasetukset gdalbuildvrtin dokumentaatiosta, esimerkiksi no-data-arvon asettamista varten.

Tiedostoluettelon luominen gdalbuildvrtille

Tiedostoluettelossa tulisi mieluiten käyttää täydellisiä polkuja, 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äydelliset 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 kuin suoritat komennon gdalbuildvrt.

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 hakee tiedostoja corine_land_cover_at_geocubes-kokoelmasta määritellyllä aikavälillä ja sijainnilla (bbox) käyttäen Asseteja (=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.