eks (version 1.1.1)

tidyst_kde: Tidy and geospatial kernel density estimates

Description

Tidy and geospatial versions of kernel density estimates for 1- and 2-dimensional data.

Usage

tidy_kde(data, ...)
st_kde(x, ...)

Value

--For tidy_kde, the output is an object of class tidy_ks, which is a tibble with columns:

x

evaluation points in x-axis (name is taken from 1st input variable in data)

y

evaluation points in y-axis (2-d) (name is taken from 2nd input variable in data)

estimate

kernel estimate values

ks

first row (within each group) contains the untidy kernel estimate from ks::kde

tks

short object class label derived from the ks object class

label

long object class label

group

grouping variable (if grouped input) (name is taken from grouping variable in data).

--For st_kde, the output is an object of class st_ks, which is a list with fields:

tidy_ks

tibble of simplified output (ks, tks, label, group) from tidy_kde

grid

sf object of grid of kernel density estimate values, as polygons, with attributes estimate, group copied from the tidy_ks object

sf

sf object of 5%, 10%, ..., 95% contour regions of kernel density estimate, as multipolygons, with attributes estimate, contlabel, contperc, contregion derived from the contour level, and group copied from the tidy_ks object.

Arguments

data

data frame/tibble of data values

x

sf object with point geometry

...

other parameters in ks::kde function

Details

For tidy_kde, the first columns of the output tibble are copied from aes(x) (1-d) or aes(x,y) (2-d). These columns are the evaluation grid points. The estimate column is the kernel density values at these grid points. The group column is a copy of the grouping variable of the input data. The ks column is a copy of the untidy kernel estimate from ks::kde, since the calculations for the layer functions geom_contour_ks, geom_contour_filled_ks require both the observations data and the kernel estimate as a kde object. For this reason, it is advised to compute a tidy kernel estimate first and then to create a ggplot with this tidy kernel estimate as the default data in the layer.

For st_kde, the output list contains the field tidy_ks which is the output from tidy_ks. The field grid is the kernel estimate values, with rectangular polygons. The field sf is the 5%, 10%, ..., 95% probability contour regions as multipolygons, with the derived attribute contlabel = 100%-cont.

The structure of the tidy_kde output is inherited from the input, i.e. if the input is a data frame/ (grouped) tibble then the output is a data frame/(grouped) tibble. Likewise for the sf object outputs for st_kde.

The default bandwidth matrix is the unconstrained plug-in selector ks::Hpi, which is suitable for a wide range of data sets, since it is not restrained to smoothing along the coordinate axes. This produces a kernel estimate which is more representative of the data than with the default bandwidth in geom_density_2d and geom_density_2d_filled. For further details of the computation of the kernel density estimate and the bandwidth selector procedure, see ?ks::kde.

If breaks is missing, then for geom_contour_ks, geom_filled_contour_ks, the density contours at cont percent are displayed, and contperc is the default aesthetic, which shows legend keys with percentages. For eks \(\geq\) 1.1.0, the aesthetic label_percent(contlabel) is deprecated in favour of contperc. The other derived quantities contregion, estimate computed via stat_filled_contour_ks, stat_contour_ks, are also suitable aesthetics and these display the contour levels and \(\geq\) contour levels in the legend keys respectively. For geom_sf, the default aesthetic is not suitable, and should be chosen from contperc, contregion, estimate, to produce analogous plots.

If breaks is not missing, then the density contours at breaks are displayed, and estimate is the default aesthetic. The other derived quantity contregion is also suitable, whereas contperc is NA so is not suitable.

For facetted plots, it is recommended to set the same breaks for all facets, and to use estimate or contregion as aesthetics. The same density contour percentages cont can lead to different contour levels for different densities in different facets which is not suitable for a facetted plot.

Examples

## tidy density estimates
library(ggplot2)
data(crabs, package="MASS")
## tidy 1-d density estimate per species
crabs1 <- dplyr::select(crabs, FL, sp)
crabs1 <- dplyr::group_by(crabs1, sp)
t1 <- tidy_kde(crabs1)
gt1 <- ggplot(t1, aes(x=FL)) 
gt1 + geom_line(colour=1) + geom_rug_ks(colour=4) + facet_wrap(~sp)

## tidy 2-d density estimate
## suitable smoothing matrix gives bimodal estimate
crabs2 <- dplyr::select(crabs, FL, CW)
t2 <- tidy_kde(crabs2)
gt2 <- ggplot(t2, aes(x=FL, y=CW)) 

## filled contour regions with density percentage labels
gt2 + geom_contour_filled_ks(colour=1) 
## filled contour regions with density level labels
gt2 + geom_contour_filled_ks(colour=1, aes(fill=after_stat(estimate)))
gt2 + geom_contour_filled_ks(colour=1, aes(fill=after_stat(contregion)))

## contour lines only with density percentage labels
gt2 + geom_contour_ks(aes(colour=after_stat(contperc)))
## contour lines only with density level labels
gt2 + geom_contour_ks(aes(colour=after_stat(estimate)))
gt2 + geom_contour_ks(aes(colour=after_stat(contregion)))

## default smoothing in geom_density_2d_filled() gives unimodal estimate
gt3 <- ggplot(crabs2, aes(x=FL, y=CW)) 
gt3 + geom_density_2d_filled(bins=4, colour=1) +
    colorspace::scale_fill_discrete_sequential(palette="Heat") +
    guides(fill=guide_legend(title="Density", reverse=TRUE))

## facet wrapped geom_sf plot with fixed contour levels for all facets
crabs3 <- dplyr::select(crabs, FL, CW, sex)
t3 <- tidy_kde(dplyr::group_by(crabs3, sex))
tb <- contour_breaks(t3, group=FALSE)
gt4 <- ggplot(t3, aes(x=FL, y=CW)) + facet_wrap(~sex)
## filled contours
gt4 + geom_contour_filled_ks(colour=1, breaks=tb, aes(fill=after_stat(estimate))) 
gt4 + geom_contour_filled_ks(colour=1, breaks=tb, aes(fill=after_stat(contregion))) 
## contour lines 
gt4 + geom_contour_ks(breaks=tb, aes(colour=after_stat(estimate)))
gt4 + geom_contour_ks(breaks=tb, aes(colour=after_stat(contregion))) 

## geospatial density estimate
data(wa)
data(grevilleasf)
hakeoides <- dplyr::filter(grevilleasf, species=="hakeoides")
hakeoides_coord <- data.frame(sf::st_coordinates(hakeoides))
s1 <- st_kde(hakeoides)

## base R plot
## filled contour regions with density percentage labels
xlim <- c(1.2e5, 1.1e6); ylim <- c(6.1e6, 7.2e6)
plot(wa, xlim=xlim, ylim=ylim)
plot(s1, add=TRUE)

## geom_sf plot
## suitable smoothing matrix gives optimally smoothed contours
gs1 <- ggplot(s1) + geom_sf(data=wa, fill=NA) + ggthemes::theme_map() 
glim <- coord_sf(xlim=xlim, ylim=ylim)
## filled contour regions with density percentage labels
gs1 + geom_sf(data=st_get_contour(s1), aes(fill=contperc)) + glim 
## filled contour regions with density level labels
gs1 + geom_sf(data=st_get_contour(s1), aes(fill=estimate)) + glim
gs1 + geom_sf(data=st_get_contour(s1), aes(fill=contregion)) + glim

## default smoothing in geom_density_2d_filled() is oversmoothed
ggplot(hakeoides_coord) + geom_sf(data=wa, fill=NA) + 
    ggthemes::theme_map() + 
    colorspace::scale_fill_discrete_sequential(palette="Heat") + 
    geom_density_2d_filled(aes(x=X, y=Y), bins=4, colour=1) +
    guides(fill=guide_legend(title="Density", reverse=TRUE)) + glim
    
## facet wrapped geom_sf plot with fixed contour levels for all facets
grevillea_gr <- dplyr::filter(grevilleasf, species=="hakeoides" |
    species=="paradoxa")
grevillea_gr <- dplyr::mutate(grevillea_gr, species=factor(species))  
grevillea_gr <- dplyr::group_by(grevillea_gr, species)
s2 <- st_kde(grevillea_gr)
sb <- contour_breaks(s2, group=TRUE)

gs3 <- ggplot(s2) + geom_sf(data=wa, fill=NA) + ggthemes::theme_map() +
    facet_wrap(~species)
## filled contour regions with density percentage labels w/o breaks
gs3 + geom_sf(data=st_get_contour(s2), aes(fill=contperc)) +
    glim 
## filled contour regions with density level labels with breaks
gs3 + geom_sf(data=st_get_contour(s2, breaks=sb), aes(fill=estimate)) +
    glim
gs3 + geom_sf(data=st_get_contour(s2, breaks=sb), aes(fill=contregion)) +
    glim 

if (FALSE) ## export as geopackage for external GIS software
sf::write_sf(wa, dsn="grevillea.gpkg", layer="wa")
sf::write_sf(hakeoides, dsn="grevillea.gpkg", layer="hakeoides")
sf::write_sf(gs1_cont, dsn="grevillea.gpkg", layer="hakeoides_cont")
sf::write_sf(s1$grid, dsn="grevillea.gpkg", layer="hakeoides_grid")