Nedávno jsem řešil docela zajímavý problém, který se týkal jádrového vyhlazování (KDE) na prostorových datech. Ačkoliv R obsahuje docela širokou škálu balíků, které KDE umí (já je hledal takhle), žádný z nich neřeší všechna nezbytná nastavení - zejména volbu jádrové funkce, dosahu a rastru na němž se KDE počítá. Pro prostorové aplikace jsou bohužel všechna zmíněná nastavení poměrně důležitá (viz třeba zde).

Existuje poměrně slušná implementace KDE, která bere ohledy na specifika jádrového vyhlazování na prostorových datech, a to v QGISu. Nástroj má poměrně slušnou nápovědu. Nabízela by se varianta použití balíku RQGIS3, který umí z R přímo používat nástroje QGISu, ale vzhledem k tomu, že mojí snahou bylo dát dohromady řešení kompletně v R, které by šlo použít i interaktivně (skrze Shiny), tak tahle varianta byla pryč ze stolu poměrně dost rychle.

Naštěstí v případě QGISu není problém nakouknout do zdrojáků konkrétních nástrojů, zde konkrétně nástroj Heatmap. Kód je celkem jednoduchý a hlavně skvěle čitelný, takže nebyl problém poměrně rychle vyrobit jeho verzi v R. Ta ale bohužel byla pro potřeby dané analýzy neskutečně pomalá. Problémem bylo velké množství bodů pro výpočet KDE a i hodně podrobný rastr, což bylo dáno zadaním úlohy. Aby vzniklo něco použitelného, bylo potřeba kód rapidně zrychlit, což se dá pouze přepsáním do nějakého “rychlejšího” jazyka. V případě R je to C++, v němž je R napsáno a k němuž mí díky balíku Rcpp skvělé rozhraní. Vyzkoušet Rcpp jsem se stejně chystal už nějakou dobu tohle byla ideální motivace. Výsledkem jednoodpoledního programování je balík SpatialKDE, který poskytuje implementaci KDE shodnou s QGIS přímo v R.

Balík obsahuje pouze 3 funkce (viz web). Přičemž jedna počítá KDE a dvě pomáhají s tvorbou podkladového rastru pro výpočet.

Ukázka

Malá ukázka, co balík umí a jak funguje. Instalace je možná buď přímo z R, pokud jste na UNIXovém systému, nebo máte nainstalované RTools.

remotes::install_github("JanCaha/SpatialKDE")

Na Windows lze použít předchystané binárky z GitHubu, které lze nainstalovat příkazem:

install.packages('SpatialKDE_x.x.x.zip', repos = NULL, type = "win.binary")

Zde je potřeba jen nahradit x.x.x aktuální verzí balíku.

Krátká ukázka, kde vyrobíme KDE pro vodní plochy v okresech Brno-město a Brno-venkov. Použijeme tu balík CzechData, který je též nutné instalovat z GitHubu.

library(tidyverse)
library(sf)
library(raster)
library(CzechData)
library(leaflet)
library(SpatialKDE)

Nejdříve získáme vrstvu vodních ploch.

water_points <- load_Data50(layer = "VodniPlocha") %>% 
  st_centroid()

A posléze dva zájmové okresy.

okresy_Brno <- load_RUIAN_state(layer = "OKRESY_P") %>% 
  filter(nazev %in% c("Brno-město", "Brno-venkov"))

Obě vrstvy přeložíme přes sebe a získáme jejich průnik, abychom získali pouze vodní plochy v zájmových okresech.

water_points_Brno <- water_points %>% 
  st_intersection(okresy_Brno)

Nastavíme rozlišení podkladového rastru a dosahu jádrové vyhlazování. V tomto případě ani jedna z hodnot nedává nějak extra smysl, jsou zvoleny poměrně náhodně, pouze s cílem vyvořit použitelnou vizualizaci.

band_width = 7500

cell_size = 5000

Vytvoříme podkladovou vrstvu nad níž budeme KDE počítat a následně vizualizovat. V tomto případě to budou hexagony a nad touto pravidelnou sítí hexagonů spočteme jádrové vyhlazení.

kde_grid <- create_raster_hexagonal(okresy_Brno, cell_size)

kde_hexagony <- kde(water_points_Brno, band_width, cell_size, 
                    kernel = "quartic", grid = kde_grid)  

## Using centroids instead of provided geometry to calculate KDE.

Vyše uvedená varovná hláška nám napovídá, že KDE se počítají pro body nikoliv plochy.

Jak výslednou vrstu KDE tak i samotné body převedeme do WGS84, pro vizualizaci skrze Leaflet.

water_points_Brno <- st_transform(water_points_Brno, 4326)
kde_hexagony <- st_transform(kde_hexagony, 4326)

Složíme výslednou vizualizaci v podobě interaktivní mapy.

col_pal <- colorNumeric("viridis", kde_hexagony$kde_value, na.color = "transparent")

leaflet(data = water_points_Brno, width = "100%") %>%
  addTiles() %>%
  addPolygons(data = kde_hexagony,
              color = "#444444", weight = 0.15, smoothFactor = 0.5,
              opacity = 1.0, fillOpacity = 0.75,
              fillColor = col_pal(kde_hexagony$kde_value)) %>% 
  addCircleMarkers(radius = 1,
                   color = "black") %>% 
  addLegend(pal = col_pal, 
            values = kde_hexagony$kde_value,
            title = "KDE vodní plochy")

Výsledná vrstva se dá samozřejmě i rasterizovat, například pomocí výborného balíku fasterize.