GRASS GIS / Propojení s R / Interpolace povrchu

Z FreeGIS portál
Přejít na: navigace, hledání

Tento článek vychází z kurzu geostat-course.org, viz prezentace, data a lokace.

Příprava

Nastavíme výpočetní region:

g.region -p rast=watc
g.region -p align=ahn5m

Dále vyplníme díry v DMT:

r.to.vect -z input=ahn5m output=ahn5m type=point
v.surf.bspline input=ahn5m raster=ahn5m_interp sie=15 sin=15 method=cubic lambda_i=0.05
Porovnání DMT před a po zaplnění děr

Oba DMT spojíme, interpolované hodnoty se použijí pouze v místech děr:

r.patch input=ahn5m,ahn5m_interp output=ahn5m_patched

Extrahování říční sítě

Extrahuje z dat říční síť:

r.mapcalc "river = if(watc == 5, 1, 0)"

Říční síť v rastrové podobě vyhladíme:

r.neighbors in=river out=river_filt method=mode size=15

Hodnoty '0' převedeme na 'no-data'.

r.null map=river_filt setnull=0

Nastavíme si pomocí modulu r.colors či interaktivně vhodnou tabulku barev Raster → Manage color rules interactively.

Nastavení tabulky barev pro vrtsvu 'river'

Dále vytvoříme rastr vzdálenosti k řece:

r.grow.distance in=river_filt dist=river_dist
Rastr vzdálenosti od řeky (v metrech)

Aktualizace bodových měření

Do atributové tabulky měření přidáme dva nové sloupečky kam uložíme hodnotu výšky a vzdálenosti od řeky.

g.copy vect=meuse@PERMANENT,meuse
v.db.addcolumn map=meuse columns="elev2 double precision,dist2 double precision"

Hodnoty vypočítáme pomocí modulu v.what.rast.

v.what.rast map=meuse raster=ahn5m_patched column=elev2
v.what.rast map=meuse raster=river_dist column=dist2

Výpočet koncentrace zinku v půdě

Spustíme y příkazové konzole GRASS GIS prostředí R a naimportujeme balíček spgrass6:

$ R
> library(spgrass6)

Načteme vektorová data:

> meuse <- readVECT6("meuse")
> meuse$zinc_log <- log1p(meuse$zinc) 
> lm.zinc <- lm(zinc_log~elev2+dist2, meuse)
> summary(lm.zinc)
> step.zinc <- step(lm.zinc)
> summary(step.zinc)
> meuse$zinc_res <- residuals(lm.zinc)
> writeVECT6(meuse, "meuse_res")

Interpolace zbytkových chyb

r.mapcalc "zinc_mask = if(isnull(dist), null(), 1)"
v.surf.rst input=meuse_res elev=meuse_res_rst maskmap=zinc_mask zcolumn=zinc_res segmax=150 npmin=600

Zpětná interpolace do skutečných hodnot

r.mapcalc "meuse_zinc = exp(b0+b1* ahn5m_patched +b2* river_dist + meuse_res_rst) - 1"

Porovnání:

v.surf.rst input=meuse_res elev=meuse_zinc_standard maskmap=zinc_mask zcolumn=zinc segmax=150 npmin=600h