Rastrová data

Rastrová data mohou být v porovnání s vektorovými daty často řádově objemnější. Tomu je třeba přizpůsobit práci s nimi. Rastrová data jsou většinou uspořádané do matice hodnot v číselné podobě.

Poznámka

Více k rastrové reprezentaci ve školení Úvod do GIS.

Pro práci s rastrovými geodaty se "tradičně" používá knihovna GDAL. Knihovna GDAL je nízkoúrovňová, přistupuje k datům pokud možno efektivním způsobem. Alternativou ke knihovně GDAL je Rasterio, která je nad touto knihovnou postavena. Jedná se o jakousi analogii ke knihovnám OGR a Fiona pro práci s vektorovými daty, viz kapitola Vektorová data.

Poznámka

Téma Rasterio je popsáno ve školení GeoPython pro začátečníky.

V současnosti knihovna GDAL podporuje více než 140 rastrových GIS formátů.

Datový model

Koncept pro rastrová data do značné míry odpovídá přístupu k vektorovým datům:

  • Driver - ovladač pro čtení a zápis dat
  • Datasource - zdroj dat, ze kterého a do kterého se čte a zapisuje
  • RasterBand - rastrový kanál. U něterých zdrojů dat je jenom jedno pásmo, ale může jich mít teoreticky neomezeně (např. u hyperspektrálních dat).
../../_images/gdal-schema.png

Mezi další důležité charakteristiky rastrových dat patří prostorové rozlišení (velikost pixelu v mapových jednotkách) a jeho hraniční souřadnice.

Užitečné odkazy:

Příklad vytvoření nového souboru z matice hodnot

V následující ukázce vytvoříme nový rastrový soubor a vyplníme ho maticí hodnot. Výsledek uložíme do souboru ve formátu GeoTIFF.

Nejprve nastavíme několik výchozích hodnot, velikost matice (mřížky), tj. velikost pixelu, číselná hodnota pro nodata (žádná) data, jméno výsledného souboru a extent (hraniční souřadnice) rastrových dat.

>>> from osgeo import gdal, ogr, osr

>>> # počet pixelů ve směru os x a y, a hodnota pro nodata
>>> pixel_size = 20
>>> NoData_value = -9999

>>> # název výstupního souboru
>>> raster_fn = 'test.tif'

>>> # hraniční souřadnice mřížky
>>> x_min, x_max, y_min, y_max = (0, 100, 0, 100)

V dalším kroku vypočteme prostorové rozlišení, velikost pixelu na základně počtu pixelů ve směru os a rozsahu rastrových dat.

>>> # prostorové rozlišení
>>> x_res = int((x_max - x_min) / pixel_size)
>>> y_res = int((y_max - y_min) / pixel_size)

Nyní můžeme vytvořit datový zdroj pro rastrová data. Nejprve vytvoříme instanci objektu Driver pro požadovaný formát a následně prázdný rastrový soubor. Zde musíme specifikovat

  • jméno výsledného souboru
  • prostorové rozlišení ve směru os x a y
  • počet pásem (kanálů)
  • typ číselné hodnoty

Nakonec nastavíme transformační parametry, které jsou ve stejném formátu v jakém bývají uloženy v tzv. world file souboru:

  • souřadnice levého-horního rohu x
  • rozlišení ve směru osy x
  • naklonění osy x
  • souřadnice levého-horního roku y
  • rozlišení ve směru osy y
  • naklonění osy y
>>> target_driver = gdal.GetDriverByName('GTiff')
>>> target_ds = target_driver.Create(raster_fn, x_res, y_res, 1, gdal.GDT_Byte)
>>> target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))

V dalším kroku zapíšeme data do vybraného pásma (číslování pásem začíná hodnotou 1 a ne více obvyklou 0). Do připraveného rastrového kanálu můžeme nyní zapsat hodnoty jako matici hodnot ve formátu NumPy.

>>> band = target_ds.GetRasterBand(1)
>>> ...
>>> import numpy as np
>>> band.WriteArray(np.array([[0, 0, 0, 0, 0],
...                  [0, 10, 15, 10, 0],
...                  [0, 15, 25, 15, 0],
...                  [0, 10, 15, 10, 0],
...                  [0, 0, 0, 0, 0]]))

Dále definujeme pro data souřadnicový systém. Ten se nastavuje pomocí zápisu ve formátu Well-known text (WKT). Souřadnicový systém definujeme pomocí kódu EPSG a vyexportujeme jako formátu WKT:

>>> outRasterSRS = osr.SpatialReference()
>>> outRasterSRS.ImportFromEPSG(5514)
>>> target_ds.SetProjection(outRasterSRS.ExportToWkt()) # !!! jiné než u vektorových dat

A nakonec uklidíme (pro jistotu) a uzavřeme zápis:

>>> band.FlushCache()

Rasterizace vektorových dat

Další ne zcela obvyklou operací může být převod vektorových dat do rastrové reprezentace. Začátek je stejný jako v předchozím případě:

>>> from osgeo import gdal, ogr, osr
>>> ...
>>> # počet pixelů ve směru os x a y, a hodnota pro nodata
>>> pixel_size = 50
>>> NoData_value = -9999
>>> ...
>>> # název výstupního souboru
>>> raster_fn = 'chko.tif'

Otevřeme vstupní vektorová data:

>>> # název vstupního vektorového souboru
>>> vector_fn = 'chko.shp'
>>> # otevření zdroje dat (DataSource)
>>> source_ds = ogr.Open(vector_fn)
>>> # načtení první vrstvy z datového zdroje
>>> source_layer = source_ds.GetLayer()

A nyní můžeme zjistit potřebné hraniční souřadnice vstupních geodat a vytvořit tak cílový rastrový soubor:

>>> # získat hraniční souřadnice
>>> x_min, x_max, y_min, y_max = source_layer.GetExtent()
>>> ...
>>> # vytvořit data source pro výstupní data
>>> x_res = int((x_max - x_min) / pixel_size)
>>> y_res = int((y_max - y_min) / pixel_size)
>>> tiff_driver = gdal.GetDriverByName('GTiff')
>>> target_ds = tiff_driver.Create(raster_fn, x_res, y_res, 3, gdal.GDT_Byte)
>>> target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))

Zkopírujeme také informaci o souřadnicovém systému (S-JTSK EPSG:5514) ze vstupního datové zdroje na výstup:

>>> outRasterSRS = osr.SpatialReference()
>>> outRasterSRS.ImportFromEPSG(5514)
>>> target_ds.SetProjection(outRasterSRS.ExportToWkt()) # !!! jiné než u vektorů

Zlatým hřebem tohoto příkladu je funkce RasterizeLayer() s následujícími parametry:

  • cílový datový zdroj
  • rastrová pásma (kanály)
  • zdrojová vektorová vrstva
  • hodnoty pro jednotlivá pásma
  • dodatečné parametry
>>> gdal.RasterizeLayer(target_ds,
    [1, 2, 3],
    source_layer,
    burn_values=[255,125,0],
    options=['ALL_TOUCHED=TRUE']) # žádné mezery okolo znaku '='
>>> target_ds.FlushCache()

Tato funkce vektorová data zrasterizuje a zapíše je do výstupního rastrového souboru.

../../_images/chko.png

Obrázek 2: Výsledek rasterizace CHKO.