Python rozhraní GRASS GIS - pygrass¶
PyGRASS je Python rozhraní open source GIS nástroje GRASS GIS. Kromě PyGRASS ještě využijeme balíček grass.notebooks
, který usnadní práci s Jupyter notebooky v prostředí GRASS GIS.
Pro zpracování dat je nutné nejprve založit nový projekt (v terminologii GRASS "location") a vstupní data do tohoto projektu naimportovat. Po úspěšném importu dat si ukážeme jak, k datovým vrtsvám přistupovat přes rozhraní PyGRASS. Přístup k projektu nastavíme v Jupyter notebooku pomocí gj.init()
.
import grass.jupyter as gj
from grass.pygrass.vector import VectorTopo
gj.init(r"C:Users\landamar\Documents\grassdata\uzpr\PERMANENT")
nuts_path = "NUTS_RG_20M_2021_3035"
nuts = VectorTopo(nuts_path)
nuts.open()
print(nuts.num_primitives())
{'point': 0, 'line': 0, 'boundary': 4669, 'centroid': 1784, 'face': 0, 'kernel': 0, 'area': 0, 'volume': 0}
Atributy a jejich vlastnosti lze zjistit pomocí table.columns
.
Poznámka: porovnejte počet atributů s arcpy.
len(nuts.table.columns)
10
for name, dtype in nuts.table.columns.items():
print(name, dtype)
cat INTEGER NUTS_ID varchar ( 5 ) LEVL_CODE INTEGER CNTR_CODE varchar ( 2 ) NAME_LATN varchar ( 70 ) NUTS_NAME varchar ( 106 ) MOUNT_TYPE INTEGER URBN_TYPE INTEGER COAST_TYPE INTEGER FID varchar ( 5 )
Záznamy v atributové tabulce lze procházet následovně:
for row in nuts.table:
print(row[3])
break
FR
count = {}
for row in nuts.table:
level = row[2]
if level not in count:
count[level] = 0
count[level] += 1
print(count)
{0: 37, 1: 125, 2: 334, 3: 1514}
Získejme přehled všech kódů jednotlivých států.
codes = []
count = 0
for row in nuts.table:
if row[3] not in codes:
codes.append(row[3])
print(sorted(codes))
['AL', 'AT', 'BE', 'BG', 'CH', 'CY', 'CZ', 'DE', 'DK', 'EE', 'EL', 'ES', 'FI', 'FR', 'HR', 'HU', 'IE', 'IS', 'IT', 'LI', 'LT', 'LU', 'LV', 'ME', 'MK', 'MT', 'NL', 'NO', 'PL', 'PT', 'RO', 'RS', 'SE', 'SI', 'SK', 'TR', 'UK']
Přístup ke geometrii prvků umožňuje objekt VectorTopo
.
code = 'CZ'
area = 0
for feat in nuts.viter('areas'):
if feat.attrs is not None and feat.attrs['CNTR_CODE'] == code and feat.attrs['LEVL_CODE'] == 0:
area += feat.area()
print(f"{code} area in km2: {area/1e6:.2f}")
CZ area in km2: 78896.87
PyGRASS umožňuje spoustět jakýkoliv nástroj, který je možný spustit v rámci grafického uživatelského rozhraní GRASS GIS. V příkladu níže v.extract.
Poznámka: Aby PyGRASS fungoval v rámci Jupyter notebooku, tak je nutné nastavit Mapset().current()
.
Poznámka: Úroveň "upovídanosti" nástrojů systému GRASS můžeme ovlivnit pomocí proměnné prostředí GRASS_VERBOSE
.
import os
os.environ["GRASS_OVERWRITE"] = "1"
os.environ["GRASS_VERBOSE"] = "-1"
from grass.pygrass.gis import Mapset
from grass.pygrass.modules import Module
Mapset().current()
for code in codes[:3]:
print("{}...".format(code))
Module("v.extract", input=nuts_path, output=f"nuts_{code}", where=f"CNTR_CODE = '{code}' AND LEVL_CODE = 2")
FR... HR... HU...
Vytvořené datové vrstvy lze zobrazit přímo v Jupyter notebooku pomocí objektu gj.Map
.
V příkladu výše je každý stát uložen v odděleném souboru ve formátu Esri Shapefile. Nyní spojíme NUTS jednotky na úroveň 0.
for code in codes[:3]:
print("{}...".format(code))
Module("v.dissolve", input=f"nuts_{code}", output=f"nuts0_{code}", column="CNTR_CODE")
map = gj.Map()
map.d_vect(map=f"nuts0_{code}")
map.show()
Rastrová data¶
Základní informace o rastrových datech poskytuje objekt RasterRow
.
dem_path = "eu_dem_v11_E40N30"
from grass.pygrass.raster import RasterRow
dem = RasterRow(dem_path)
dem.open()
dem.info.rows, dem.info.cols
(40000, 40000)
Zájmové území vybereme na základě NUTS jednotky a vrstvu digitálního modelu terénu na základě toho ořežeme pomocí nástroje r.mask.
nuts_id = 'CZ042'
aoi_path = "nuts_aoi"
Module("v.extract", input=nuts_path, output=aoi_path, where=f"NUTS_ID = '{nuts_id}'")
dem_aoi_path = "dmt_aoi"
Module("g.region", vector=aoi_path, align=dem_path)
Module("r.mask", vector=aoi_path)
Module("r.mapcalc", expression=f"{dem_aoi_path} = {dem_path}")
dem_aoi = RasterRow(dem_aoi_path)
dem_aoi.open()
dem_aoi.info.rows, dem_aoi.info.cols
(4426, 4619)
Na základě rastrových dat můžeme vytvořit numpy pole a dále pro výpočet používat knihovnu Numpy.
import numpy
array = numpy.array(dem_aoi)
array.shape, array.dtype
((4426, 4619), dtype('float32'))
Získejme minimální a maximální nadmořskou výšku pro dané zájmové území.
numpy.nanmin(array), numpy.nanmax(array)
(49.590366, 1228.7755)
Pomocí numpy můžeme provést např. reklasifikaci hodnot.
array[array<=200] = 0
array[array>200] = 1
numpy.unique(array, return_counts=True) # numpy >= 1.24 equal_nan=True
(array([ 0., 1., nan], dtype=float32), array([ 913759, 7596521, 11933414]))
Výsledek reklasifikace rastrových hodnot zapíšeme na disk.
from grass.pygrass.raster import numpy2raster
numpy2raster(array.astype('int16'), "CELL", "dem200", overwrite=True)
map = gj.Map()
map.d_rast(map="dem200")
map.show()