Knihovna Rasterio

Knihovna Rasterio je opět dílo zejména Seana Gilliese (podobně jako Fiona či Shapely), tentokrát v rámci jeho působení ve firmě MapBox. Rasterio je knihovna pro práci s rastrovými geografickými datovými sadami. Na pozadí Rasterio používá spolehlivou knihovnu GDAL.

Rasterio pracuje s objekty knihovny NumPy (podobně jako dříve zmíněná Fiona pracuje s objekty JSON). Autor tvrdí, že Rasterio se vyslovuje [raw-STEER-ee-oh] a měla by práci s rastrovými daty udělat více zábavnou a produktivnější.

V následujícím příkladu otevřeneme rastrový soubor ve formátu GeoTIFF a podíváme se na některá metadata:

import rasterio
src = rasterio.open('data/lsat7_2002_nir.tiff')
print (src.bounds)
BoundingBox(left=596670.0, bottom=185000.0, right=678330.0, top=258500.0)
print (src.crs)
{u'lon_0': -79, u'datum': u'NAD83', u'y_0': 0, u'no_defs': True,
u'proj': u'lcc', u'x_0': 609601.22,
u'units': u'm', u'lat_2': 34.33333333333334, u'lat_1': 36.16666666666666,
u'lat_0': 33.75}
print (src.tags())
{u'AREA_OR_POINT': u'Area'}
print (src.width, src.height)
(1287, 831)
print (src.res)
(10.0, 10.0)
../../_images/rgb.png

Obrázek 1: RGB kompozice.

Načtení barevných kanálů:

data = src.read()
print (len(data))
3

Vidíme, že v rastru jsou obsaženy tři barevné kanály. Vytvoříme nyní nový soubor, obsahující pokus o index NDVI.

Poznámka

Normalizovaný vegetační index je poměr mezi viditelným červeným kanálem a blízkým infračerveným kanálem satelitních dat.

\[NDVI = (NIR - VIS) / (NIR + VIS)\]

Neprve vytvoříme novou matici pro výsledné hodnoty, následně do tohoto pole uložíme výsledek výpočtu pro každý pixel. Pracujeme v prostředí NumPy, které práci s poli významně usnadňuje.

(nir, vis) = (data[0], data[1])
ndvi = (nir - vis) / (nir + vis)
print (ndvi.min(), ndvi.max())
-0.94444442, 0.97435898

Výsledek uložíme do nově vytvořeného souboru. Data budou zkomprimována pomocí LWZ komprese a uložena v číselném formátu float64 (rastrový soubor obsahuje čísla s plovoucí desetinnou čárkou a záporné hodnoty). Výsledný soubor ve formátu GeoTIFF bude mít pouze jeden kanál.

kwargs = src.meta
kwargs.update(dtype=rasterio.float64, count=1, compress='lzw')
with rasterio.open('data/ndvi.tif', 'w', **kwargs) as dst:
    dst.write_band(1, ndvi.astype(rasterio.float64))
../../_images/ndvi.png

Obrázek 2: Výsledný soubor s NDVI indexem.