Rdev/statistics/cluster_analysis_spat.R at master · zoometh/Rdev · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
library(sf) # pour le spatial
library(ggplot2) # pour les graphiques
library(ggrepel) # pour écarter les labels
library(dplyr) # méthode

# data: n-columns
depots <- read.table('https://raw.github.com/zoometh/Rdev/master/data/data.csv',
                     header = T,
                     sep = ";",
                     row.names = 1)
symbols.col <- c("shape", "color", "x", "y")
depots.symbols <- depots[ , (names(depots) %in% symbols.col)]
depots.ca <- depots[ , !(names(depots) %in% symbols.col)]
# créer le dendrogramme
dend1 <- depots.ca %>%
  dist %>% hclust(method = "complete")
# coupe en 3 groupes
groups <- as.data.frame(cutree(dend1, k = 3))
colnames(groups)[1] <- "group"
depots.group <- merge(depots, groups, by = "row.names")
# download FRA.zip (le shapefile)
temp <- tempfile()
download.file("https://raw.github.com/zoometh/Rdev/master/data/FRA.zip",
              destfile = temp, quiet = TRUE)
td <- tempdir()
lfiles <- unzip(temp, exdir = td) # tous les fichiers (.shp, .dbf, etc.)
FRA <- st_read(dsn = td, layer = "FRA_adm0")
# recoupe FRA sur l'emprise des depots +/- buffer
buff <-  0.5
xmin <- min(depots.group$x) - buff
xmax <- max(depots.group$x) + buff
ymin <- min(depots.group$y) - buff
ymax <- max(depots.group$y) + buff
m <- rbind(c(xmin,ymin), c(xmax,ymin), c(xmax,ymax), c(xmin,ymax), c(xmin,ymin))
roi <- st_polygon(list(m))
roi <- st_sfc(roi)
st_crs(roi) <- "+init=epsg:4326"
FRA.roi <- st_intersection(FRA, roi)
# titre
tit <- paste("CAH sur", nrow(depots.ca), "individus et", ncol(depots.ca), "variables")
# graphique
gcah.sp <- ggplot() +
  ggtitle(tit) +
  geom_sf(data = FRA.roi, fill = 'gray90') +
  geom_point(data = depots.group, aes(x = x, y = y, color = as.factor(group))) +
  geom_text_repel(data = depots.group,
                  aes(x = x, y = y, label = Row.names, color = as.factor(group)),
                  cex = 2,
                  segment.size = 0.1,
                  segment.alpha = 0.5) +
  theme(plot.title = element_text(size = 8, face = "bold")) +
  theme_bw() +
  theme(legend.position="bottom")
gcah.sp

# sauver
png("out/cah_depots_spat.png", width = 12, height = 10, units = "cm", res = 300)
gcah.sp
dev.off()
shell.exec(paste0(getwd(), "/out/cah_depots_spat.png"))