用 ggplot2 畫 boxplot

這個簡單的例子是用 ggplot2 畫 boxplot

setwd('~/Documents/PhD/thesis/figs/draw_env_boxplot/')
library(RPostgreSQL)
library(ggplot2)

建立 database connection,我想把 nvdimp_uplots 這個資料表中的 plotid, twsvf, soil, mtrl 這四個欄位的值取出,並且關連至我論文分類出來的分類群資料表(syntaxa_groups),並根據 plotid 的點位,去抽取 dtm, slope 和 aspect 這三個 raster 的數值,並輸出為 hc_env

conn = dbConnect(PostgreSQL(), user= "username", dbname="database")
db_query = "
select 
    p.plotid, 
    h.group_no,
    st_value(s.rast, h.geom_twd97) AS elevation,
    st_value(s2.rast, h.geom_twd97) AS slope,
    st_value(s3.rast, h.geom_twd97) AS aspect,
    twsvf openness,
    soil soil_rockiness,
    mtrl rock_cover
from 
    public.nvdimp_uplots p, 
    hc.syntaxa_groups h,
    env.tw_dtm_asterv2 s,
    env.tw_slope_asterv2 s2,
    env.tw_aspect_asterv2 s3
where 
    p.plotid=h.plotid AND
    st_intersects(s.rast, h.geom_twd97) AND
    st_intersects(s2.rast, h.geom_twd97) AND
    st_intersects(s3.rast, h.geom_twd97)
order by group_no;"

qres = dbSendQuery(conn, db_query)
hc_env = fetch(qres)

上面的 SQL 查詢出來的表大概長得像下面(資料有 mask 過)

plotid group_no elevation slope aspect openness soil_rockiness rock_cover
01-0000 1 3358 2.5 150 998 50 5
02-0000 1 3537 9.6 70 964 60 10
03-0000 1 3326 45.6 340 975 100 15
04-0000 1 3625 11.9 345 995 30 20
05-0000 1 3554 17.3 346 1000 40 25
06-0000 1 3863 5 106 989 50 30
07-0000 1 3412 15.5 293 999 60 35
08-0000 1 3611 16.6 10 996 70 40

因為查詢出的表中 group_no 欄位類型是類別資料(factor),我想把這欄的類別資料轉成文字,所以我用 factor() 來做

# factorize the association group_no and rename to abbreviation
hc_env$group_no =factor(hc_env$group_no, levels=c(1,2,3,4,5,6,7,8,9),
       labels=c("Ger-Jun", "Aco-Jun", "Jun-Abi", "Yus-Abi",
                "Tsu-Abi", "Yus-Tsu", "Pin-Tsu", "Rho-Tsu", "Ell-Pic"))

接下來,就是用 ggplot 畫圖,先把顏色和欄位名稱都設定好

# draw boxplot
col_palette = terrain.colors(9)
env_names = c("elevation", "slope", "aspect", "openness",
              "terrain", "soil_rockiness", "rock_cover")
env_labels = c("Elevation (m)", "Slope (degree)", "Aspect (degree)",
              "Openness", "Terrain", "Soil Rockiness (%)", "Rock cover (%)")

因為我每個 boxplot y 軸尺度都有點不太一樣,因此我寫了個簡單的 function 去呼叫 ggplot,一般來說可以用 qplot() 來直接畫圖,但是我希望有更多客製化的選項,所以用疊圖層的方式來繪製

draw_boxplot = function(x, y, y_lower,y_upper) {
    x$y = x[,y]
    # 個別呼叫 ggplot 的元件
    p = ggplot(x, aes(x=group_no, y = y, fill=group_no))
    # 希望 x 軸的標籤轉 90 度,顏色設定為黑色,並且不要畫出圖例
    p + theme(axis.text.x = element_text(angle=90, colour="black"),
          axis.text.y = element_text(colour="black"),
          legend.position="none") +
    \# 設定 y 的數值上下界
    ylim(y_lower,y_upper)+
    xlab('') +
    # 標上 y 軸說明
    ylab(env_labels[y-2]) +
    # 繪製 boxplot,先統計出最大最小值、中位數等資料,再呼叫 geom_boxplot() 畫圖
    stat_boxplot(geom ='errorbar')  +
    geom_boxplot() + 
    # 手動填上指定的顏色
    scale_fill_manual(values=col_palette)
}

因為 ggplot 似乎無法用 par(mfrow()) 去設定,所以要另外寫 function 來繪製,請參考 Cookbook for R: multiple graphs on one page (Reference 2)

# 輸出成 pdf
pdf(file='env_factors.pdf', width=8, height=8)
# multiplot() 這個 function 
multiplot(draw_boxplot(hc_env, 3, 2500, 4000),
          draw_boxplot(hc_env, 4, 0, 60),
          draw_boxplot(hc_env, 5, 0, 360),
          draw_boxplot(hc_env, 6, 700, 1000),
          draw_boxplot(hc_env, 8, 0, 100),
          draw_boxplot(hc_env, 9, 0, 100),
          cols=3)
dev.off()

匯出圖的成果就像這樣:

References:

  1. Use ggplot2 within another function in R
  2. R cookbook multiple graphs on one page
  3. R cookbook - ggplot2 Colors
  4. put whisker on ggplot2 boxplot

Comments

comments powered by Disqus