r R pca绘图ggplot2

Posted 2021-05-24

# plot multi-bands features
library(tidyr)
library(ggplot2)
library(ggalt)
source('R/multiplot.r')
csvfile = paste('datasets/', gsub(".R",".csv", basename(rstudioapi::getActiveDocumentContext()$path)), sep='')
print(csvfile)
df = read.csv(csvfile)

df$AI = (df$A8 - df$A4) / (df$A8 + df$A4)
df$BI = (df$B8 - df$B4) / (df$B8 + df$B4)

pca_df = prcomp(df[3:12])
df_pc = data.frame(pca_df$x, class=df$class)
df_pc_crop = df_pc[df_pc$class=='Croplands',]
df_pc_nc = df_pc[df_pc$class=='Non-Croplands',]

p2 = ggplot(df_pc, aes(PC1, PC2, col=class)) + 
  geom_point(aes(shape=class), size=2) +   # draw points
  labs(subtitle="(b)Principal components PC1 and PC2 as X and Y axis") + 
  coord_cartesian(xlim = 1.2 * c(min(df_pc$PC1), max(df_pc$PC1)),
                  ylim = 1.2 * c(min(df_pc$PC2), max(df_pc$PC2))) +   # change axis limits
  theme_bw() +
  theme(legend.position="top") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
  # geom_encircle(data = df_pc_crop, aes(x=PC1, y=PC2)) +   # draw circles
  # geom_encircle(data = df_pc_nc, aes(x=PC1, y=PC2))

# print(p2)
xlabels = strsplit('B2.B3.B4.B8.NDVI', '[.]')
df = df %>% gather(bands,values,-id,-class)
df$class=factor(df$class)
p1 = ggplot(df, aes(x = bands, y = values, color=class)) +
  geom_boxplot(outlier.shape = NA)  + 
  labs(subtitle="(a) Knowledge Feature") + 
  ylim(0,.7) +
  scale_x_discrete(labels=xlabels) +
  annotation_custom(
    grob = textGrob(label = 'season 1', hjust = 0, gp = gpar(cex = 1.5)),
    ymin = 0,      # Vertical position of the textGrob
    ymax = 0,
    xmin = 1,         # Note: The grobs are positioned outside the plot area
    xmax = 3) +
  annotation_custom(
    grob = textGrob(label = 'season 2', hjust = 0, gp = gpar(cex = 1.5)),
    ymin = 0,      # Vertical position of the textGrob
    ymax = 0,
    xmin = 1,         # Note: The grobs are positioned outside the plot area
    xmax = 12) +
  theme_bw() +
  theme(legend.position="top") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
  # geom_point(position = position_jitter(width = 0.05))

multiplot(p1, p2, cols=2)