Použití oken pro čtení a zápis

Rastrová data lze číst a zapisovat i pomocí tzv. „oken“. Tento postup se hodí v případě, že pracujete s objemnými rastrovými daty, které se do RAM vašeho počítače nevejdou. Pomocí okna lze načíst část matice rastrových dat, tu zpracovat a zapsat do výstupního souboru.

Okna (rasterio.Window) jsou pravidelné matice vstupního rastru. Lze je také popsat jako 2 páry souřadnic pixelů:

((první_řádek, poslední_řádek), (první_sloupec, poslední_sloupec))

nebo:

Window(první_sloupec, první_řádek, šířka, výška)

Příklad načtení dat:

import rasterio
from rasterio.windows import Window

with rasterio.open('data/B04-2018-05-06.tiff') as red:
    w = red.read(1, window=Window(0, 0, 512, 256))

print(w.shape)

Bloky

Rastrové soubory ukládají většinou data po blocích. Pokud chcete využívat okna pro čtení, je efektivní, aby velikost oken odpovídala velikosti bloků. Následujícím způsobem můžeme zpracovat první kanál našeho rastrového souboru po blocích:

with rasterio.open('data/B04-2018-05-06.tiff') as red:
    for ji, window in red.block_windows(1):
        r = red.read(1, window=window)
        print(r.shape)

To jakým způsobem jsou bloky pro rastrový soubor definovány můžete zjistit z příkazové řádky nástrojem gdalinfo (součást knihovny GDAL).

gdalinfo data/B04-2018-05-06.tiff

Driver: GTiff/GeoTIFF
Files: B04-2018-05-06.tiff
    B04-2018-05-06.tiff.aux.xml
Size is 3117, 1716

...
Band 1 Block=3117x8 Type=UInt16, ColorInterp=Gray

Vidíme, že náš rastr používá bloky o velikosti 3117×8 pixel - tedy 8 řádků. Někdy může být efektivnější celý rastr převzorkovat, než se pustíte do jeho blok-po-bloku processingu. To můžeme udělat např. nástrojem knihovny GDAL gdalwarp (viz školení Úvod Open Source GIS).

gdalwarp -r mode -co TILED=YES -co BLOCKXSIZE=256 -co BLOCKYSIZE=256 data/B04-2018-05-06.tiff outputs/B04-2018-05-06-256block.tiff
...

gdalinfo outputs/B04-2018-05-06-256block.tiff
...
Band 1 Block=256x256 Type=UInt16, ColorInterp=Gray

A otevřít v Pythonu:

with rasterio.open('outputs/B04-2018-05-06-256block.tiff') as red:

    for ji, window in red.block_windows(1):
        r = red.read(1, window=window)
        print(r.shape)

NDVI pomocí bloků

Následující příklad prezentuje opět výpočet NDVI, tentokrát ale bude postupovat po blocích (čtení i zápis). Tento přístup bude fungovat rychleji při větším objemu dat.

import rasterio
from rasterio.windows import Window

with rasterio.open('data/B04-2018-05-06.tiff') as red:
    with rasterio.open('data/B08-2018-05-06.tiff') as nir:

        step = 256
        kwargs = red.meta
        kwargs.update(dtype=rasterio.float32, count=1, compress='lzw')

        with rasterio.open('outputs/ndvi_w.tif', 'w', **kwargs) as dst:
            slices = [(col_start, row_start, step, step) \
                        for col_start in list(range(0, red.width, 256)) \
                        for row_start in list(range(0, red.height, 256))
            ]

            # nepoužijeme block windows, protože data používají 1x1287 velké bloky
            # for ji, window in src.block_windows(1):
            for slc in slices:
                win = Window(*slc)

                nir_data = nir.read(1, window=win).astype(float)
                vis_data = red.read(1, window=win).astype(float)

                ndvi = (nir_data - vis_data) / (nir_data + vis_data)

                write_win = Window(slc[0], slc[1], ndvi.shape[1], ndvi.shape[0])

                dst.write_band(1, ndvi.astype(rasterio.float32), window=write_win)