Using {sf} operations

Because {densityarea} can return {sf} polygons, this can allow you to use its functionality (see this cheat sheet for some examples).

Density polygons as simple features

As a first example, we’ll estimate how much different data clusters overlap in the s01 dataset.

data(s01)

s01 |> 
  mutate(lF1 = -log(F1),
         lF2 = -log(F2)) ->
  s01

We’ll focus on the vowels iy, ey, o and oh which correspond to the following lexical classes:

vowel label	lexical class
`iy`	Fleece
`ey`	Face
`o`	Lot
`oh`	Thought

s01 |> 
  filter(
    plt_vclass %in% c("iy", 
                      "ey", 
                      "o", 
                      "oh")
  ) ->
  vowel_subset

Within this subset of vowel categories, we’ll get the 80% probability density estimate as sf::st_polygon()s.

vowel_subset |> 
  group_by(plt_vclass) |> 
  reframe(
    density_polygons(lF2, 
                     lF1, 
                     probs = 0.8,
                     as_sf = TRUE)
  ) |> 
  st_sf()->
  vowel_polygons

vowel_polygons
#> Simple feature collection with 4 features and 3 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -8.050893 ymin: -7.037138 xmax: -6.996794 ymax: -5.820745
#> CRS:           NA
#> # A tibble: 4 × 4
#>   plt_vclass level_id  prob                                             geometry
#>   <chr>         <int> <dbl>                                            <POLYGON>
#> 1 ey                1   0.8 ((-7.896047 -6.222708, -7.897749 -6.224405, -7.9032…
#> 2 iy                1   0.8 ((-7.833113 -5.836786, -7.843041 -5.828972, -7.8529…
#> 3 o                 1   0.8 ((-7.351372 -6.23354, -7.365699 -6.233222, -7.38002…
#> 4 oh                1   0.8 ((-7.240256 -6.419138, -7.249486 -6.417758, -7.2587…

We can plot these directly by using the sf::geom_sf() geom for ggplot2.

ggplot(vowel_polygons) +
  geom_sf(
    aes(fill = plt_vclass),
    alpha = 0.6
  ) +
    scale_fill_brewer(palette = "Dark2")

Initial `{sf}` operations

All of the {sf} operations for geometries are available to use on vowel_polygons. For example, we can get the area of each polygon, with sf::st_area() and use it in plotting.

vowel_polygons |> 
  mutate(
    area = st_area(geometry)
  ) ->
  vowel_polygons

ggplot(vowel_polygons) +
  geom_sf(
    aes(fill = area),
    alpha = 0.6
  )+
  scale_fill_viridis_c()

Or, we can get the polygon centroids and plot them.

vowel_polygons |> 
  st_centroid() |> 
  ggplot()+
    geom_sf_label(
      aes(label = plt_vclass,
          color = plt_vclass,
          size = area)
    )+
    scale_color_brewer(palette = "Dark2")+
    coord_fixed()
#> Warning: st_centroid assumes attributes are constant over geometries

Getting overlaps

To use the density polygons like “cookie cutters” on each other, we need to use st_intersections().

vowel_polygons |> 
  st_intersection() -> 
  vowel_intersections

vowel_intersections
#> Simple feature collection with 6 features and 6 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -8.050893 ymin: -7.037138 xmax: -6.996794 ymax: -5.820745
#> CRS:           NA
#> # A tibble: 6 × 7
#>   plt_vclass level_id  prob   area n.overlaps origins                   geometry
#>   <chr>         <int> <dbl>  <dbl>      <int> <list>                   <POLYGON>
#> 1 ey                1   0.8 0.0865          1 <int>   ((-7.63623 -6.326645, -7.…
#> 2 ey                1   0.8 0.0865          2 <int>   ((-7.896883 -6.458929, -7…
#> 3 iy                1   0.8 0.202           1 <int>   ((-7.902604 -6.461049, -7…
#> 4 o                 1   0.8 0.285           1 <int>   ((-7.236759 -6.940033, -7…
#> 5 o                 1   0.8 0.285           2 <int>   ((-7.093493 -6.493589, -7…
#> 6 oh                1   0.8 0.224           1 <int>   ((-7.092569 -6.49308, -7.…

This data frame contains a polygon for each unique intersection of the input polygons, with a new n.overlaps column.

ggplot(vowel_intersections) +
  geom_sf(
    aes(fill = n.overlaps)
  )+
  scale_fill_viridis_c()

The labels of the new overlapping regions aren’t very informative, but we can create some new labels by using the indices in the origins column.

new_label <- function(indices, labels){
  str_c(labels[indices],
        collapse = "~")
}

vowel_intersections |> 
  mutate(
    groups = map_chr(
      origins,
      .f = new_label,
      labels = vowel_polygons$plt_vclass
    ) 
  ) |> 
  relocate(groups, .after = plt_vclass)->
  vowel_intersections

vowel_intersections
#> Simple feature collection with 6 features and 7 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -8.050893 ymin: -7.037138 xmax: -6.996794 ymax: -5.820745
#> CRS:           NA
#> # A tibble: 6 × 8
#>   plt_vclass groups level_id  prob   area n.overlaps origins  
#>   <chr>      <chr>     <int> <dbl>  <dbl>      <int> <list>   
#> 1 ey         ey            1   0.8 0.0865          1 <int [1]>
#> 2 ey         ey~iy         1   0.8 0.0865          2 <int [2]>
#> 3 iy         iy            1   0.8 0.202           1 <int [1]>
#> 4 o          o             1   0.8 0.285           1 <int [1]>
#> 5 o          o~oh          1   0.8 0.285           2 <int [2]>
#> 6 oh         oh            1   0.8 0.224           1 <int [1]>
#> # ℹ 1 more variable: geometry <POLYGON>

ggplot(vowel_intersections)+
  geom_sf(
    aes(fill = groups)
  )+
  scale_fill_brewer(palette = "Dark2")

We can also calculate the areas of these new polygons, and compare them to the original areas (which have been preserved in areas.

vowel_intersections |> 
  mutate(
    group_area = st_area(geometry),
    overlapped_proportion = 1-(group_area/area)
  ) |> 
  filter(n.overlaps == 1) |> 
  ggplot(
    aes(plt_vclass, overlapped_proportion)
  )+
    geom_col()+
    ylim(0,1)

Spatial filters

There are also a number of spatial filters and merges that can be used interestingly if the original data points are also converted to sf objects.

library(sfheaders)

library(forcats)

s01 |> 
  sfheaders::sf_point(
    x = "lF2",
    y = "lF1",
    keep = TRUE
  ) ->
  s01_sf

s01_sf
#> Simple feature collection with 4245 features and 10 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -8.095446 ymin: -7.335764 xmax: -6.541463 ymax: -5.439817
#> CRS:           NA
#> First 10 features:
#>    name age sex     word vowel plt_vclass ipa_vclass    F1     F2  dur
#> 1   s01   y   f     OKAY    EY        eyF        ejF 763.5 2088.1 0.20
#> 2   s01   y   f       UM    AH         uh          ʌ 699.5 1881.3 0.19
#> 3   s01   y   f      I'M    AY         ay         aj 888.8 1934.1 0.07
#> 4   s01   y   f    LIVED    IH          i          ɪ 555.5 1530.5 0.05
#> 5   s01   y   f       IN    IH          i          ɪ 612.2 2323.4 0.06
#> 6   s01   y   f COLUMBUS    AH          @          ə 612.4 1903.7 0.07
#> 7   s01   y   f       MY    AY         ay         aj 578.4 1959.3 0.09
#> 8   s01   y   f   ENTIRE    IH          i          ɪ 529.9 2332.1 0.08
#> 9   s01   y   f   ENTIRE    ER        *hr         ə˞ 538.4 1682.8 0.18
#> 10  s01   y   f     LIFE    AY        ay0        aj0 744.6 1702.1 0.15
#>                       geometry
#> 1   POINT (-7.64401 -6.637913)
#> 2  POINT (-7.539718 -6.550366)
#> 3  POINT (-7.567397 -6.789872)
#> 4   POINT (-7.33335 -6.319869)
#> 5  POINT (-7.750787 -6.417059)
#> 6  POINT (-7.551555 -6.417386)
#> 7  POINT (-7.580343 -6.360266)
#> 8  POINT (-7.754524 -6.272688)
#> 9  POINT (-7.428214 -6.288602)
#> 10 POINT (-7.439618 -6.612847)

s01_sf |> 
  ggplot()+
    geom_sf()

Next, we can get the density polygon for a single vowel,.

s01 |> 
  filter(plt_vclass == "iy") |> 
  reframe(
    density_polygons(lF2, lF1, probs = 0.8, as_sf =T)
  ) |> 
  st_sf()->
  iy_sf

Spatial filter

Let’s get all of the points in s01_sf that are “covered by” the iy_sf polygon.

s01_sf |> 
  st_filter(
    iy_sf,
    .predicate = st_covered_by
  )->
  covered_by_iy

covered_by_iy |> 
  mutate(plt_vclass = plt_vclass |> 
           fct_infreq() |> 
           fct_lump_n(5)) |> 
  ggplot()+
    geom_sf(data = iy_sf)+
    geom_sf(aes(color = plt_vclass))+
    scale_color_brewer(palette = "Dark2")

Obviously, the 80% probability density area for iy is not a homogeneous region.

covered_by_iy |> 
  mutate(plt_vclass = plt_vclass |> 
           fct_infreq() |> 
           fct_lump_n(5)) |> 
  count(plt_vclass) |> 
  ggplot(aes(plt_vclass, n))+
    geom_col(
      aes(fill = plt_vclass)
    )+
    scale_fill_brewer(palette = "Dark2")

Spatial join

Let’s see which vowel category a random vowel token is close to.

set.seed(100)
s01_sf |> 
  slice_sample(n = 1)->
  rand_vowel

rand_vowel
#> Simple feature collection with 1 feature and 10 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -7.440205 ymin: -6.508322 xmax: -7.440205 ymax: -6.508322
#> CRS:           NA
#>   name age sex word vowel plt_vclass ipa_vclass    F1     F2  dur
#> 1  s01   y   f LIKE    AY        ay0        aj0 670.7 1703.1 0.09
#>                      geometry
#> 1 POINT (-7.440205 -6.508322)

Then, we’ll get density polygons at a few different probability points for all vowels.

s01 |>
  group_by(plt_vclass) |>
  reframe(
    density_polygons(
      lF2,
      lF1,
      probs = ppoints(5),
      range_mult = 0.5,
      as_sf = T
    )
  ) |> 
  st_sf() ->
  vowel_probs

There are 5 probability level polygons for each vowel category in vowel_probs. We join the random vowel’s data onto this set of polygons with st_join().

vowel_probs |> 
  st_join(
    rand_vowel,
    .predicate = st_covers
  ) |> 
  drop_na()->
  vowel_within

Now, for each vowel category in this new data frame, let’s get the smallest probability polygon (i.e. where the random point is closest to the center probability mass).

vowel_within |> 
  group_by(plt_vclass.x) |> 
  filter(prob == min(prob)) |> 
  ungroup() |> 
  mutate(plt_vclass = fct_reorder(plt_vclass.x, prob)) ->
  vowel_min_prob

vowel_min_prob |> 
  ggplot()+
    geom_sf(aes(fill = prob)) +
    geom_sf(data = rand_vowel |> mutate(plt_vclass = NULL),
            color = "red") +
    facet_wrap(~plt_vclass)