library(sf)
## Linking to GEOS 3.7.1, GDAL 2.4.0, PROJ 5.2.0
library(scales)
options(warn = -1)
# Input data
dat = data.frame(
group1 = c("A", "A", "A", "A", "A", "A", "B", "B", "B", "B", "C", "C", "C", "D", "D", "D", "D"),
group2 = c("a", "b", "b", "c", "c", "c", "c", "d", "e", "f", "f", "g", "h", "h", "i", "a", "a"),
group3 = 1:17,
count = c(2, 2, 1, 1, 1, 1, 2, 1, 2, 3, 2, 3, 3, 2, 4, 1, 1),
stringsAsFactors = FALSE
)
dat
## group1 group2 group3 count
## 1 A a 1 2
## 2 A b 2 2
## 3 A b 3 1
## 4 A c 4 1
## 5 A c 5 1
## 6 A c 6 1
## 7 B c 7 2
## 8 B d 8 1
## 9 B e 9 2
## 10 B f 10 3
## 11 C f 11 2
## 12 C g 12 3
## 13 C h 13 3
## 14 D h 14 2
## 15 D i 15 4
## 16 D a 16 1
## 17 D a 17 1
# Create circle points
ctr = st_point(c(0, 0))
ctr = st_sfc(ctr)
circle = st_buffer(ctr, dist = 5, nQuadSegs = 360 / 4)
circle = st_sfc(circle)
pnt = st_cast(circle, "POINT")
pnt = pnt[1:360]
# Rotate so that zero is on top
rot = function(a) matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)
pnt = (pnt - ctr) * rot(-pi/2)
# Create rays
x = split(pnt, 1:length(pnt))
x = lapply(x, st_geometry)
x = lapply(x, "[[", 1)
x = lapply(x, function(x) c(ctr[[1]], x))
x = lapply(x, st_cast, "LINESTRING")
x = do.call(c, x)
x = st_sfc(x)
x = st_cast(x, "LINESTRING")
# Calculate section polygons
pol = st_buffer(ctr, dist = 3)
dat$start = cumsum(dat$count) - dat$count
dat$end = cumsum(dat$count)
dat$start1 = rescale(dat$start, to = c(1, 361), from = range(c(dat$start, dat$end)))
dat$end1 = rescale(dat$end, to = c(1, 361), from = range(c(dat$start, dat$end)))
dat$start1 = round(dat$start1)
dat$end1 = round(dat$end1)
dat$start = NULL
dat$end = NULL
dat$start1[dat$start1 == 361] = 1
dat$end1[dat$end1 == 361] = 1
result = list()
for(i in 1:nrow(dat)) {
ray1 = x[dat$start1[i]]
ray2 = x[dat$end1[i]]
section = c(ray1, ray2)
section = st_union(section)
section = st_convex_hull(section)
section = st_intersection(pol, section)
result[[i]] = section
}
result = do.call(c, result)
# Combine with group IDs
pol = st_sf(dat, geometry = result)
# Split to groups
group1 = aggregate(pol[, "count"], st_drop_geometry(pol["group1"]), sum)
group2 = aggregate(pol[, "count"], st_drop_geometry(pol["group2"]), sum)
group3 = aggregate(pol[, "count"], st_drop_geometry(pol["group3"]), sum)
# Remove inner circles
group1 =
group1 %>%
st_intersection(st_buffer(ctr, dist = 1))
group2 =
group2 %>%
st_intersection(st_buffer(ctr, dist = 2)) %>%
st_difference(st_buffer(ctr, dist = 1))
group3 =
group3 %>%
st_intersection(st_buffer(ctr, dist = 3)) %>%
st_difference(st_buffer(ctr, dist = 2))
# Combine
group1$level = 1
colnames(group1)[1] = "group"
group2$level = 2
colnames(group2)[1] = "group"
group3$level = 3
colnames(group3)[1] = "group"
dat = rbind(group1, group2, group3)
dat
## Simple feature collection with 30 features and 3 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: -3 ymin: -3 xmax: 3 ymax: 3
## epsg (SRID): NA
## proj4string: NA
## First 10 features:
## group count level geometry
## 1 A 8 1 POLYGON ((1 -6.661338e-16, ...
## 2 B 8 1 POLYGON ((1.160755e-14 -1, ...
## 3 C 8 1 POLYGON ((-1 -2.943235e-14,...
## 4 D 8 1 POLYGON ((6.123234e-17 1, 0...
## 5 a 4 2 POLYGON ((0.7812242 1.84044...
## 6 b 3 2 POLYGON ((1.65757 1.118045,...
## 7 c 5 2 POLYGON ((1.840449 -0.78122...
## 8 d 1 2 POLYGON ((1.65757 -1.118045...
## 9 e 2 2 POLYGON ((1.118045 -1.65757...
## 10 f 5 2 POLYGON ((-0.7812242 -1.840...
#########################################################################
library(ggplot2)
# Set colors
n1 = unique(dat$group[dat$level == 1])
n2 = unique(dat$group[dat$level == 2])
n3 = unique(dat$group[dat$level == 3])
cols = c(
hcl.colors(length(n1), "Pastel 1"),
hcl.colors(length(n2), "Heat 2"),
hcl.colors(length(n3), "Dark 2")
)
names(cols) = c(n1, n2, n3)
# Plot
ggplot() +
geom_sf(data = dat, aes(fill = group)) +
geom_sf_text(data = dat, aes(label = group)) +
scale_fill_manual(values = cols, breaks = names(cols)) +
theme_bw() +
theme(
axis.title = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank(),
panel.border = element_blank()
) +
guides(fill = guide_legend(ncol = 2))
