植物与生境关联的检验: Torus translation简介

2019-05-01 修订

1 前言

由于进化以及适应性的差异，很多物种对生境有一定的偏好。在亚热带森林生态系统中，不同生境中可能生长着不同的物种，如山脊上常有杜鹃花属(Rhododendron)、南烛属(Lyonia)植物生长； 而阴湿的沟谷中，则常见杨梅叶蚊母树(Distylium myricoides)、山茶属(Camellia)植物、茜草科植物等；在经常受干扰的地段，常会有马尾松(Pinus massoniana)、枫香(Liquidambar formosana)、木荷(Schima superba)等较为喜光的树种。

物种对生境的偏好，需要生境关联检验进行分析。假设某样地的生境分为： A.沟谷，B.山坡，C. 山脊三种类型，要检验一种目标植物是否倾向于在C山脊上生长，却不能直接用卡方检验 (Chi Square Test)。这是因为卡方检验的前提是物种的个体之间是完全独立的，换句话说说，就是个体在群落中是完全随机分布的，个体与个体之间没有什么联系。但这个假设是没有办法满足的：由于传播限制等因素的影响， 不同个体在空间上往往表现出一定的聚集性。例如，幼苗和小树往往只能在母树附近生长。如果此时仍然按照卡方检验去判断物种是否对某种生境具有偏好， 就容易高估物种对生境的关联。从这个角度上看，过去大部分探讨物种和生境关联的研究，都需要用新方法检验。

Harms, K. E等人在2001年的一篇论文里提出了Torus translation的方法，以检验物种和生境的关联，本文就要介绍这个检验。虽然该方法已经提出十几年了，但在大样地物种分布格局的研究中仍然十分重要。

Torus translation其实是一种随机化方法，但是是有条件的随机化。Torus是拓扑学专业词汇，中文为“环面”，torus translation是指在随机化过程中，所研究的样方位置，好像在环面上滚动一样，保存样方之间的空间关系。

Torus translation大体过程如下：

假设现在有一块320mX500m的森林样地，如果以20m边长的样方为研究尺度，那么就有400个边长为20m的样方。每一个20m样方的生境已经按照地形划分为A、B、C、D四种类型。现在要研究某个种是否与生境存在关联，也就是说，这个种的个体数是否在某种生境中出现的更多，而在另外一种生境中出现地更少。

首先，计算目标种在每种生境中的个体数，同时计算样地中所有种在每种生境中的个体数，两者相除，得到目标种在每种生境中的实际相对多度（或者翻译为相对密度）。

现在，让每个样方中目标种的个体数固定，而生境类型按照行或者列，每次增加20m。如果生境的样方超过了样地的坐标范围，则放在最前面一排，或者最前面一列，以此类推。在这个过程中，任何一个样方在任意变换中，除了一条或者两条边以外，和其他周围样方的相对关系都不改变，从而最大限度上保留了该物种所有个体在空间上的聚集程度。 经过这样的位置变换之后，重新计算目标种在每种生境中的个体数，就可以计算出在该次torus translation格局下，目标种在每种生境中的相对多度。

Torus translation可以进行多少次？ 举例来说，320mX500m的样方，以20m为尺度，可以进行399次torus translation。不过真实的生境分布，还可以衍生出以下三种情况：横轴相反的镜像 （X mirror image）， 纵轴相反的镜像 （180° rotation of mirror image），横纵轴同时反向（相当于旋转180度， 180° rotation）。每种情况都可以再进行一套完整的torus translation。 这样一来， 对于320mX500m的样地，就有(400*4 - 1), 即1599种torus translation。每一次torus translation， 针对每一种生境，都可以比较物种实际相对多度相比基于torus translation转换以后获得的1599个相对多度，再基于95%置信区间做出进一步判断。

在某一种生境中，可以把这1599个基于torus translation的计算出来的相对多度按照从高到低排序。如果实际观测值小于torus translation 2.5%分位数，也就是排名在40位之前，那么我们就在alpha = 0.05的显著性水平断定目标物种不喜欢这种生境， 或者说排斥这种生境。如果实际观测值大于torus translation 97.5%分位数，也就是排名在1559位之后，那么我们就能在alpha = 0.05的显著性水平断定目标物种偏好这种生境。如果该物种的相对多度位于torus translation 这1599个值的95%置信区间中，就可以认为物种对这种生境没有什么特别排斥或者偏好，换言之，是随机的。

当然，如果物种的个体数比较少，经过torus translation，或者有时候某种生境的面积很少，就会出现在某种生境中，能够比较的相对多度数量很少的情况，此时permutation test的检验能力可能会不足。因此，一般来说，torus translation test都需要物种的个体数比较多。如Kyle Harms等人的研究中， 就只选用了多过50个体的种来检验。

torus translation 检验的本质， 是基于条件控制下的随机排列来计算显著性水平的方法。 正是对样方之间相对空间关系的严格控制，才保证零假设不是完全随机。如果让生境在样方中完全随机排列，所得结果将与卡方检验相差无几。

torus translation的算法虽然比较简单，但是直接编写函数，对于R初学者仍然比较困难。虽然生态学家Kyle E. Harms在他的个人网站也公布了相应的函数和相关数据 （https://sites01.lsu.edu/faculty/kharms/wp-content/uploads/sites/23/2014/05/R_TorusTranslationTest_HarmsK_011911_For_txt_Trees.doc ），但是从数据格式上看，这个函数需要一定的操作技巧。同时，由于函数一直没有整合到程序包中，R初学者在使用时也有一定困难。

前不久，CTFS下辖的ForestGEO森林监测网络重新整理了这些函数，将样地数据分析相关的重要函数整合到一系列程序包中，其中fgeo.habitat程序包（ https://github.com/forestgeo) ）就可以进行torus translation检验，使用起来也较为方便。

下面的例子是fgeo.habitat程序包tt_test函数的示例代码。笔者修改了该示例代码，增加了绘图部分，希望结果torus translation能更直观， 以方便理解。示例代码之后，是tt_test和torusonesp.all函数的详细注释。为了格式统一，笔者对源代码进行了一些调整，如参数中的逻辑变量，F，T统一改为 FALSE，TRUE； if条件判断中， x，y条件的顺序调整为一致；if函数后仅有一行命令的，也增加了花括号等。

由于本人水平有限，错误之处在所难免，敬请读者批评指正。

2 R代码

2.1 导入程序包

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
library(fgeo.habitat)library(fgeo.data)
## 
## Attaching package: 'fgeo.data'
## The following objects are masked from 'package:fgeo.habitat':
## 
##     luquillo_habitat, luquillo_stem6_random, luquillo_tree6_random
library(fgeo.base)
## 
## Attaching package: 'fgeo.base'
## The following objects are masked from 'package:fgeo.data':
## 
##     luquillo_stem_random_tiny, luquillo_vft_4quad
library(fgeo.tool)
## 
## Attaching package: 'fgeo.tool'
## The following objects are masked from 'package:fgeo.base':
## 
##     %||%, check_crucial_names, detect_insensitive, drop_if_na,
##     drop_status, extract_insensitive, flag_if, guess_plotdim,
##     is_duplicated, is_multiple, pick_dbh_max, pick_dbh_min,
##     pick_dbh_over, pick_dbh_under, pick_status, rename_matches
## The following objects are masked from 'package:fgeo.habitat':
## 
##     extract_gridsize, extract_plotdim, guess_plotdim
## The following object is masked from 'package:stats':
## 
##     filter
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
### luquillo 样地位于波多黎各，这里用样地三种植物的分布数据，检验生境关联
habitat <- luquillo_habitat
census  <- luquillo_top3_sp
``

## 2.2 绘制生境分布图

```R
m <- ggplot()+ 
  geom_raster(data = habitat, aes(x = gx + 10, y = gy + 10, fill = as.factor(habitats))) + 
  geom_hline(yintercept = sort(unique(habitat$gy)), color = "grey") + 
  geom_vline(xintercept = sort(unique(habitat$gx)), color = "grey") + 
  scale_y_discrete(expand = c(0,0)) + 
  scale_x_discrete(expand = c(0,0)) + 
  xlab("Easting (m)") + 
  ylab("Northing (m)") + 
  scale_fill_discrete(name = "Habitat", labels = c("A", "B", "C", "D")) 

m  + 
  ggtitle(paste("Actual Habitat Distribution of Luquillo Plot"))
img

2.3 物种分布图

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
# 筛选出胸径10以上的树
pick <- filter(census, status == "A", dbh >= 10)
# 必须要达到50个个体以上再进行检验
pick2 <- add_count(pick, sp)
pick3 <- filter(pick2, n > 50)

species <- unique(pick3$sp)

sp1 <- subset(pick3, sp == species[1] )
sp2 <- subset(pick3, sp == species[2] )
sp3 <- subset(pick3, sp == species[3] )

m + geom_point(data = sp1, mapping = aes(x = gx, y = gy, size = dbh ), shape = 1) + 
scale_shape(solid = FALSE) + ggtitle(paste("Distribution of ", species[1]))
img

m + geom_point(data = sp2, mapping = aes(x = gx, y = gy, size = dbh ), shape = 1) + 
scale_shape(solid = FALSE) + 
ggtitle(paste("Distribution of ",species[2]))
img

m + geom_point(data = sp3, mapping = aes(x = gx, y = gy, size = dbh ), shape = 1) + 
scale_shape(solid = FALSE) + 
ggtitle(paste("Distribution of ",species[3]))
img

2.4 三种变换

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
有三种情况, 是没有在循环中考虑的
1. X Mirror 左右颠倒
2. flip Y 上下颠倒
3. 180° rotate 180°旋转，左右和上下同时颠倒

2.4.1 X-mirror
```R
    newx <- -(habitat$gx - max(habitat$gx))
    newy <-   habitat$gy
    mirror_pattern <- cbind(habitat, newx = newx, newy = newy)
ggplot()+  geom_raster(data = mirror_pattern, aes(x = newx + 10, y = newy + 10, 
fill = as.factor(habitats))) + 
  geom_hline(yintercept = sort(unique(mirror_pattern$newy)), color = "grey") + 
  geom_vline(xintercept = sort(unique(mirror_pattern$newx)), color = "grey") + 
  scale_y_discrete(expand = c(0,0)) + 
  scale_x_discrete(expand = c(0,0)) + 
  xlab("Easting (m)") + 
  ylab("Northing (m)") + 
  scale_fill_discrete(name = "Habitat", labels = c("A", "B", "C", "D")) + 
  ggtitle(paste("Torus translation", "X-Mirror")) + 
  scale_colour_hue(c=45, l=80)
img

2.4.2 Y-flip

1
2
3
4
5
6
7
8
9
10
11
12
13
14
    newx <- habitat$gx
    newy <- -(habitat$gy - max(habitat$gy))
    mirror_pattern <- cbind(habitat, newx = newx, newy = newy)
ggplot()+  geom_raster(data = mirror_pattern, aes(x = newx + 10, y = newy + 10, fill = as.factor(habitats))) + 
  geom_hline(yintercept = sort(unique(mirror_pattern$gy)), color = "grey") + 
  geom_vline(xintercept = sort(unique(mirror_pattern$gx)), color = "grey") + 
  scale_y_discrete(expand = c(0,0)) + 
  scale_x_discrete(expand = c(0,0)) + 
  xlab("Easting (m)") + 
  ylab("Northing (m)") + 
  scale_fill_discrete(name = "Habitat", labels = c("A", "B", "C", "D")) + 
  ggtitle(paste("Torus translation", "Y-Flip")) + 
  scale_colour_hue(c=45, l=80)
img

2.4.3 180度旋转

1
2
3
4
5
6
7
8
9
10
11
12
13
14
    newx <- -(habitat$gx - max(habitat$gx))
    newy <- -(habitat$gy - max(habitat$gy))
    mirror_pattern <- cbind(habitat, newx = newx, newy = newy)
ggplot()+  geom_raster(data = mirror_pattern, aes(x = newx + 10, y = newy + 10, fill = as.factor(habitats))) + 
  geom_hline(yintercept = sort(unique(mirror_pattern$gy)), color = "grey") + 
  geom_vline(xintercept = sort(unique(mirror_pattern$gx)), color = "grey") + 
  scale_y_discrete(expand = c(0,0)) + 
  scale_x_discrete(expand = c(0,0)) + 
  xlab("Easting (m)") + 
  ylab("Northing (m)") + 
  scale_fill_discrete(name = "Habitat", labels = c("A", "B", "C", "D")) + 
  ggtitle(paste("Torus translation", "Rotate 180")) + 
  scale_colour_hue(c=45, l=80)
img

2.5 Torus translation 核心算法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
### 变换样方原点的坐标， 如超出样地横纵坐标x或者y的最大值， 就减样地样地x或者y的最大值。
increment <- unique(habitat[2:nrow(habitat), c("gx", "gy")])
nrow(increment)
## [1] 399
res <- list()for (i in 1:nrow(increment)){
    newx <- habitat$gx + as.numeric(increment[i,1]) 
    newy <- habitat$gy + as.numeric(increment[i,2])    
    for (j in 1:length(newx)){        
    if(newx[j] > max(habitat$gx)){
            newx[j] <- newx[j] - max(habitat$gx)
        } 
    }    
    for(k in 1:length(newy)){        
    if(newy[k] > max(habitat$gy)){
            newy[k] <- newy[k] - max(habitat$gy)
        }
    }
    res[[i]] <- cbind(habitat, newx = newx, newy = newy)
}

2.6 绘制经过Torus Translation随机化之后的生境分布图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
### ### 这里从400种Torus Translation随机化中，随机抽出10个查看
for (i in sort(sample(1:length(res), 10))){
temp <- res[[i]]
p <- ggplot()+  geom_raster(data = temp, aes(x = newx + 10, y = newy + 10, 
fill = as.factor(habitats))) + 
  geom_hline(yintercept = sort(unique(temp$gy)), color = "grey") + 
  geom_vline(xintercept = sort(unique(temp$gx)), color = "grey") + 
  scale_y_discrete(expand = c(0,0)) + 
  scale_x_discrete(expand = c(0,0)) + 
  xlab("Easting (m)") + 
  ylab("Northing (m)") + 
  scale_fill_discrete(name = "Habitat", labels = c("A", "B", "C", "D")) + 
  ggtitle(paste("Torus translation", i)) + 
  scale_colour_hue(c=45, l=80)
print(p)
}

2.7 tt_test函数实例

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
head(census)
## # A tibble: 6 x 19
##   treeID stemID tag   StemTag sp    quadrat    gx    gy MeasureID CensusID
##    <int>  <int> <chr> <chr>   <chr> <chr>   <dbl> <dbl>     <int>    <int>
## 1    180    225 1001… 100174  CASA… 921     165.   410.    617049        6
## 2    631    775 10069 10069   PREM… 213      38.3  245.    598429        6
## 3   1380   1702 1015… 101560  CASA… 820     142.   386.    614023        6
## 4   1840   2240 10208 10208   PREM… 613     116.   245.    607825        6
## 5   2849   3421 1031… 103156  CASA… 420      79.2  389.    603814        6
## 6   3354   4054 1037… 103756  CASA… 220      32.0  385.    599273        6
## # … with 9 more variables: dbh <dbl>, pom <chr>, hom <dbl>,
## #   ExactDate <dbl>, DFstatus <chr>, codes <chr>, nostems <dbl>,
## #   status <chr>, date <dbl>
head(species)
## [1] "CASARB" "PREMON" "SLOBER"
head(habitat)
## # A tibble: 6 x 3
##      gx    gy habitats
##   <dbl> <dbl>    <int>
## 1     0     0        1
## 2     0    20        1
## 3     0    40        1
## 4     0    60        1
## 5     0    80        1
## 6     0   100        1
habitat$x <- habitat$gx
habitat$y <- habitat$gy# Test any number of species (output a list of matrices)
tt_lst <- tt_test(census, sp = unique(census$sp), habitat)

str(tt_lst, give.attr = FALSE)
## List of 3
##  $ : num [1, 1:24] 29 1242 356 2 0 ...
##  $ : num [1, 1:24] 91 1093 504 3 0 ...
##  $ : num [1, 1:24] 18 273 1324 3 0 ...
tt_lst[[1]]
##        N.Hab.1 Gr.Hab.1 Ls.Hab.1 Eq.Hab.1 Rep.Agg.Neut.1 Obs.Quantile.1
## CASARB      29     1242      356        2              0        0.77625
##        N.Hab.2 Gr.Hab.2 Ls.Hab.2 Eq.Hab.2 Rep.Agg.Neut.2 Obs.Quantile.2
## CASARB      20      390     1206        4              0        0.24375
##        N.Hab.3 Gr.Hab.3 Ls.Hab.3 Eq.Hab.3 Rep.Agg.Neut.3 Obs.Quantile.3
## CASARB      12      778      817        5              0        0.48625
##        N.Hab.4 Gr.Hab.4 Ls.Hab.4 Eq.Hab.4 Rep.Agg.Neut.4 Obs.Quantile.4
## CASARB       5      932      658       10              0         0.5825
# Try also: View((to_df(tt_lst)))
head(to_df(tt_lst))
##           metric     sp      value
## 1        N.Hab.1 CASARB   29.00000
## 2       Gr.Hab.1 CASARB 1242.00000
## 3       Ls.Hab.1 CASARB  356.00000
## 4       Eq.Hab.1 CASARB    2.00000
## 5 Rep.Agg.Neut.1 CASARB    0.00000
## 6 Obs.Quantile.1 CASARB    0.77625

3 函数注释

3.1 tt_test函数注解

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
## check_tt 函数用于检查
## census数据集内，所需信息是否齐全
## 提供的参数必须包括census, sp, habitat, plotdim, gridsize四个参数
## 而torusonesp.all 函数所需要的abundance矩阵，tt_test()会根据census, plotdim 和 gridsize自动生成
tt_test <- function(census,
                    sp,
                    habitat,
                    plotdim = extract_plotdim(habitat),
                    gridsize = extract_gridsize(habitat)) {


  check_tt(
    census = census, sp = sp, habitat = habitat, plotdim = plotdim,
    gridsize = gridsize
  )

  abundance <- abund_index(census, plotdim, gridsize)
  out <- lapply(
    X = sp,
    FUN = torusonesp.all,
    allabund20 = abundance,
    hab.index20 = habitat,
    plotdim = plotdim,
    gridsize = gridsize
  )
  new_tt_lst(out)
}

3.2 torusonesp.all函数注解

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
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182

### species     要进行tt检验的是哪一个种
### hab.index20 20m方格生境组成类型
### allabund20  物种多度矩阵，通常应该由abund_index 函数生成
### plotdim     样地的大小，长和宽
### gridsize    样地方格的大小，默认为20
mtorusonesp.all <- function(species, hab.index20, allabund20, plotdim, gridsize) {  
# 计算x轴方向上，20m样方的数量
  plotdimqx <- plotdim[1] / gridsize  # 计算y轴方向上，20m样方的数量
  plotdimqy <- plotdim[2] / gridsize  # 生境类型数
  num.habs <- length(unique(hab.index20$habitats))  
  ### GrLsEq: 是Greater, Less, Equal的缩写。用来保存最终的结果
  ### 每种生境类型，都进一步分为六种情况，放在6个列中，分别为：
  ### N.Hab  在真实生境和模拟生境中都没有出现， 出现的次数都是0， 不能比较
  ### Gr.Hab 在真实生境中个体数高于模拟生境中个体数的次数
  ### Ls.Hab 在真实生境中个体数低于模拟生境中个体数的次数
  ### Eq.Hab 在真实生境中个体数等于模拟生境中个体数的次数
  ### Rep.Agg.Neut 排斥-0， 偏好1， 无显著关系：0
  ### Obs.Quantile 该种在该生境下的相对多度
  ###
  GrLsEq <- matrix(0, 1, num.habs * 6) # 空白矩阵

  ### 只有一行，用来保存目标种的结果
  rownames(GrLsEq) <- species 

  ### 依次为每一列创建列名， 生境按照出现的次序，用i表示。所以
  for (i in 1:num.habs) 
  {    if (i == 1) { ## 利用生境1初始化各列
      cols <- c(paste("N.Hab.", i, sep = ""), 
      paste("Gr.Hab.", i, sep = ""), 
      paste("Ls.Hab.", i, sep = ""), 
      paste("Eq.Hab.", i, sep = ""), 
      paste("Rep.Agg.Neut.", i, sep = ""), 
      paste("Obs.Quantile.", i, sep = ""))
    }    if (i > 1) { 
      ## 用 test <- c(test, add) 的方法，逐步增加向量的长度，
      ## 一直到所有的生境类型都加进去了为止。
      cols <- c(cols, paste("N.Hab.", i, sep = ""), 
                paste("Gr.Hab.", i, sep = ""), 
                paste("Ls.Hab.", i, sep = ""), 
                paste("Eq.Hab.", i, sep = ""), 
                paste("Rep.Agg.Neut.", i, sep = ""), 
                paste("Obs.Quantile.", i, sep = ""))
    }
  }
  colnames(GrLsEq) <- cols 
  # 要输出的结果更改列名，这样计算结果就可以逐步写入该matrix了。

  # 基于真实生境分布计算实际相对多度

  # 从所有种的多度数据中提取目标种在每个20m样方中的个体数。
  allabund20.sp <- allabund20[which(rownames(allabund20) == species), ]  
  # 将目标物种的多度数据转换为矩阵， 
  # 矩阵的行数为y轴上20m的样方数， 
  # 列为x轴方向上的20m样方数。
  # 值为每个20m样方中，该种的个体数
  spmat <- matrix(as.numeric(allabund20.sp), 
                  nrow = plotdimqy, ncol = plotdimqx, byrow = FALSE)  
                  # 将群落中所有物种的多度数据转换为矩阵，
  #行列设置与spmat相同。值为每个20m样方中，所有种的个体数
  totmat <- matrix(apply(allabund20, MARGIN = 2, FUN = "sum"), 
                   nrow = plotdimqy, ncol = plotdimqx, byrow = FALSE)  
                   # 每个种的生境类型转换为矩阵。行列设置与spmat相同
  habmat <- matrix(hab.index20$habitats, nrow = plotdimqy, 
                   ncol = plotdimqx, byrow = FALSE)  
                   # 空向量，用于保存每种生境中，每个种的全部个体数
  spstcnthab  <- numeric()  
  # 空向量，用于保存每种生境中，全部种的个体数
  totstcnthab <- numeric()  
  for (i in 1:num.habs)
  {    # 各生境中，全部种的个体数
    totstcnthab[i] <- sum(totmat[habmat == i])    
    # 各生境中，目标种的个体数
    spstcnthab[i] <- sum(spmat[habmat == i])
  }  # 在真实生境中，目标种的相对多度
  spprophab <- spstcnthab / totstcnthab  
  # 计算Torus translation转换以后的相对多度
  habmat.template <- habmat  for (j in 1:4) {    
  # 当j == 1时， habmat <- habmat, 所以不用写这一句， 
    # 直接进入torus translation阶段即可。

    ### x轴方向反转
    if (j == 2) {
      habmat <- apply(habmat.template, 2, rev)
    }    ### y轴方向反转
    if (j == 3) {
      habmat <- t(apply(habmat.template, 1, rev))
    }    ### x和y轴同时反转
    if (j == 4) {
      habmat <- t(apply(apply(habmat.template, 2, rev), 1, rev))
    }    #### 针对以上三种情形，分别进行torus translation随机化
    # 从x轴开始translation
    for (x in 0:(plotdimqx - 1)){      
    # 从y轴开始translation
      for (y in 0:(plotdimqy - 1)){        
      # 创建一个全部为0的矩阵， 用于保存生境随机化后的生境
        newhab <- matrix(0, plotdimqy, plotdimqx)        
        # 作者是用下标提取矩阵的一部分，再通过合并下标，组成新的矩阵。
        if (x == 0 & y == 0) {
          newhab <- habmat
        }        if (x > 0  & y == 0) {
          newhab <- habmat[c(1:plotdimqy),
              c((plotdimqx - x + 1):plotdimqx, 1:(plotdimqx - x))]
        }        if (x == 0 & y > 0) {
          newhab <- habmat[c((plotdimqy - y + 1):plotdimqy, 1:(plotdimqy - y)),
              c(1:plotdimqx)]
        }        if (x > 0 & y > 0) {
          newhab <- habmat[c((plotdimqy - y + 1):plotdimqy, 1:(plotdimqy - y)),
              c((plotdimqx - x + 1):plotdimqx, 1:(plotdimqx - x))]
        }        
        # 空向量， 用于保存各生境在torus translation后， 20m样方中每个种的相对多度.
        Torspstcnthab  <- numeric() 
        
        # 空向量， 用于保存各生境在torus translation后， 20m样方中所有种总的相对多度.
        Tortotstcnthab <- numeric() 
        
        for (i in 1:num.habs){          
        # 各生境中，所有种的个体数
          Tortotstcnthab[i] <- sum(totmat[newhab == i])          
          # 各生境中，目标种的个体数
          Torspstcnthab[i] <- sum(spmat[newhab == i])
        }        
        # 经过torus translation转换之后，每种生境中，目标种的相对多度
        Torspprophab <- Torspstcnthab / Tortotstcnthab        
        ## 这部分为每种生境中物种真实个体数与随机化之后的
        ## 个体数比较的结果分别计数
        ## 情况如下
        for (i in 1:num.habs)
        {          
        ### 每个种在某一种生境中的相对多度，
          ### 如果不能和所有种的相对多度比较，则这种比较无意义
          if (is.na(spprophab[i] > Torspprophab[i])) {
            warn_invalid_comparison(spprophab[i], Torspprophab[i])
          }          
          ### 如果目标种的在某生境中的相对多度高于
          ### torus translation转换的结果， 则为greater than这一列加1
          if (spprophab[i] > Torspprophab[i]) {
            GrLsEq[1, (6 * i) - 4] <- GrLsEq[1, (6 * i) - 4] + 1
          }          ### 如果目标种的在某生境中的相对多度低
          ### 于torus translation转换的结果， 则为less than这一列加1
          if (spprophab[i] < Torspprophab[i]) {
            GrLsEq[1, (6 * i) - 3] <- GrLsEq[1, (6 * i) - 3] + 1
          }          ### 如果目标种的在某生境中的相对多度等于
          ### torus translation转换的结果， 则为equal to这一列加1
          if (spprophab[i] == Torspprophab[i]) {
            GrLsEq[1, (6 * i) - 2] <- GrLsEq[1, (6 * i) - 2] + 1
          }
        }
      } #  x方向的20m translation结束
    } # y方向的20m translation结束} # x镜像， y镜像，180度旋转


  for (i in 1:num.habs)
  {   ## 每一种生境的第二列
    GrLsEq[1, (6 * i) - 5] <- spstcnthab[i] 
    
    ### 这里的结果是输出到Rep.Agg.Neut的

    ### 若实际生境中，物种个体密度小于torus translation
    ### 转换后随机个体密度的2.5%
    ### 就表示这个物种在统计上显著排斥这种生境，结果用-1表示
    if (GrLsEq[1, (6 * i) - 4] / (4 * (plotdimqx * plotdimqy)) <= 0.025){
      GrLsEq[1, (6 * i) - 1] <- -1 #
    }    ### 若实际生境中，物种个体密度大于torus translation 
    ### 转换后随机个体密度整体分布的的97.5%
    ### 表示这个物种在统计上显著偏好这种生境，结果用1表示

    if (GrLsEq[1, (6 * i) - 4] / (4 * (plotdimqx * plotdimqy)) >= 0.975){
      GrLsEq[1, (6 * i) - 1] <- 1
    }    ### 如果实际生境中， 物种个体数介于torus translation 
    ### 转换后95%的置信区间之中,
    ### 则在统计上没有显著差异，此时就看不出排斥或者偏好这种生境，用0表示
    # otherwise it's neutral (not different from random dist)
    if (GrLsEq[1, (6 * i) - 4] / (4 * (plotdimqx * plotdimqy)) < 0.975 & GrLsEq[1, (6 * i) - 4] / (plotdimqx * plotdimqy) > 0.025) {
        GrLsEq[1, (6 * i) - 1] <- 0
    }   # 下面是最后一列
   # 物种的实际相对多度
    GrLsEq[1, (6 * i)] <- GrLsEq[1, (6 * i) - 4] / (4 * (plotdimqx * plotdimqy))
  }  # 返回值
  return(GrLsEq)
}

4 参考文献：

Harms, K. E., Condit, R., Hubbell, S. P., & Foster, R. B. (2001). Habitat associations of trees and shrubs in a 50-ha Neotropical forest plot. Journal of Ecology, 89(6), 947–959.

Mauro Lepore (NA). fgeo.habitat: Analize Soils and Tree-Habitat Data. R package version 0.0.0.9006. https://github.com/forestgeo/fgeo.habitat

天文年历中几种常见天象图的绘制

中科院紫金山天文台和北京天文馆《天文爱好者》杂志社曾连续多年编纂《天文普及年历》， 其中就有不少行星周年视运动、日月和大行星的升没时刻和中天时刻图、木星伽利略卫星动态等。类似的图形，在《香港天文年历》中也有出现， 如 http://www.weather.gov.hk/gts/astron2018/2018transit.pdf ， http://www.weather.gov.hk/gts/astron2018/2018rise-set.pdf

这里提供一些R代码， 用于绘制类似图形。这些图形有些已经用于近几年的《爱牧夫天文台历》等 （参见 http://bbs.imufu.cn/thread-689955-1-1.html ）。由于学识有限，错误和疏漏在所难免， 恳请广大读者指正。

用到的R程序包及其安装
天象图全部是用R语言程序包 skycalc, moonsun和ephem， 假设已安装了devtools程序包以及Rtools(针对Windows用户)， 可以用以下命令安装上述三个程序包，需要按照如下顺序安装：

先安装Rcpp

install.packages(“Rcpp”)

skycalck安装：

devtools::install_github(“helixcn/skycalc”)

ephem包安装：

devtools::install_github(“helixcn/ephem”)

moonsun程序包: install.packages(“moonsun”)

R脚本
library(ephem)

Loading required package: moonsun

Loading required package: skycalc

Loading required package: Rcpp

library(skycalc)
options(longitude = 120) ## 假设经度为东经120度
金星和水星在日出日落时的高度和方位
venus_sunset(CAADate_DateToJD(2018, 1, 1, TRUE),
latitude = 40, days = 365, label_interval = 10 )
2018年金星在日落时的高度和方位

2018年金星在日落时的高度和方位

venus_sunrise(CAADate_DateToJD(2018, 1, 1, TRUE),
latitude = 40, days = 365, label_interval = 10)
2018年金星在日出时的高度和方位

2018年金星在日出时的高度和方位

mercury_sunrise( CAADate_DateToJD(2018, 1, 1, TRUE) ,
latitude = 40, days = 365, label_interval = 10)
2018年水星在日出时的高度和方位

2018年水星在日出时的高度和方位

mercury_sunset(CAADate_DateToJD(2018, 1, 1, TRUE),
latitude = 40, days = 365, label_interval = 10 )
2018年水星在日落时的高度和方位

2018年水星在日落时的高度和方位

2018年日月和大行星每月中旬在星空背景中的位置
plot_location_atlas(CAADate_DateToJD(Year = 2018,
Month =1,
Day =15,
TRUE))

Warning in log(constellation$MAG): NaNs produced

Warning in log(constellation$MAG): NaNs produced

Warning in log(constellation$MAG): NaNs produced

2018年1月15日月和大行星在星空背景中的位置

2018年1月15日月和大行星在星空背景中的位置

ra d phase angle dist size mag

2018-01-15-Sun 19h 46m 53s -22* 49’ 43’’ NA NA 0.98 0.5 NA

2018-01-15-Moon 18h 9m 11s -24* 38’ 15’’ 3.9 89.2 1.05 29.6 NA

2018-01-15-Mercury 18h 23m 22s -24* 33’ 59’’ 85.4 87.4 1.25 5.4 -1.5

2018-01-15-Venus 19h 51m 40s -22* 4’ 13’’ 100.0 -55.9 1.71 9.9 -3.9

2018-01-15-Mars 15h 18m 55s -18* 32’ 58’’ 92.1 105.3 1.82 5.1 0.9

2018-01-15-Jupiter 15h 5m 39s -17* 41’ 9’’ 99.3 106.4 5.76 34.2 -1.9

2018-01-15-Saturn 18h 11m 18s -23* 28’ 40’’ 100.0 91.1 10.99 15.1 1.3

2018-01-15-Uranus 1h 33m 44s 9* 11’ 8’’ 99.9 -111.4 19.90 3.3 5.8

2018-01-15-Neptune 22h 52m 20s -9* 50’ 40’’ 100.0 -111.5 30.66 2.0 7.9

2018-01-15-Pluto 18h 45m 15s -21* 55’ 52’’ 100.0 97.1 33.40 0.2 14.2

2018-01-15-Sun Sgr

2018-01-15-Moon Sgr

2018-01-15-Mercury Sgr

2018-01-15-Venus Sgr

2018-01-15-Mars Lib

2018-01-15-Jupiter Lib

2018-01-15-Saturn Sgr

2018-01-15-Uranus Psc

2018-01-15-Neptune Aqr

2018-01-15-Pluto Sgr

2018年日月和大行星中天时刻
plot_transit_table(transit_table(
CAADate_DateToJD(Year = 2018,
Month =1,
Day =1,
TRUE),
days =365),
legend = TRUE)
2018年日月和大行星的中天时刻

2018年日月和大行星的中天时刻

2018年行星的周年视运动
library(moonsun)
data(starcat)
ephem.mercury <- mercury(jd(year = 2018,
month = 1,
day = 1,
length = 365))
track(ephem.mercury,
mag = 4,
interval.lab = 30,
bright.lab = TRUE,
starcat = starcat,
bright = bright)
2018年水星的周年视运动

2018年水星的周年视运动

ephem.venus <- venus(jd(year = 2018,
month = 1,
day = 1,
length = 365))
track(ephem.venus, mag = 4,
interval.lab = 30,
bright.lab = TRUE,
starcat = starcat,
bright = bright)
2018年金星的周年视运动

2018年金星的周年视运动

ephem.mars <- mars(jd(year = 2018,
month = 1,
day = 1,
length = 365))
track(ephem.mars,
mag = 4,
interval.lab = 30,
bright.lab = TRUE,
starcat = starcat,
bright = bright)
2018年火星的周年视运动

2018年火星的周年视运动

ephem.jupiter <- jupiter(jd(year = 2018,
month = 1,
day = 1,
length = 365))
track(ephem.jupiter,
mag = 6,
interval.lab = 30,
bright.lab = TRUE,
starcat = starcat,
bright = bright)
2018年木星的周年视运动

2018年木星的周年视运动

ephem.saturn <- saturn(jd(year = 2018,
month = 1,
day = 1,
length = 365))
track(ephem.saturn,
mag = 8,
interval.lab = 30,
bright.lab = TRUE,
starcat = starcat,
bright = bright)
2018年土星的周年视运动

2018年土星的周年视运动

ephem.uranus <- uranus(jd(year = 2018,
month = 1,
day = 1,
length = 365))
track(ephem.uranus,
mag = 8,
interval.lab = 30,
bright.lab = TRUE,
starcat = starcat,
bright = bright)
2018年天王星的周年视运动

2018年天王星的周年视运动

ephem.neptune <- neptune(jd(year = 2018,
month = 1,
day = 1,
length = 365))
track(ephem.neptune,
mag = 8,
interval.lab = 30,
bright.lab = TRUE,
starcat = starcat,
bright = bright)
2018年海王星的周年视运动

2018年海王星的周年视运动

ephem.pluto <- pluto(jd(year = 2018,
month = 1,
day = 1,
length = 365))
track(ephem.pluto,
mag = 8,
interval.lab = 30,
bright.lab = TRUE,
starcat = starcat,
bright = bright)
2018年冥王星的周年视运动

2018年冥王星的周年视运动

木星伽利略卫星的位置变化
plot_juptermoons(CAADate_DateToJD(Year = 2018,
Month = 1,
Day =1,
TRUE),
days = 31, increment = 1/24)
2018年1月木星伽利略卫星的位置变化

2018年1月木星伽利略卫星的位置变化

参考文献

Komsta, Lukasz. 2013. Moonsun: Basic Astronomical Calculations with R. https://CRAN.R-project.org/package=moonsun.

Zhang, Jinlong. 2015. Skycalc: Computational Astronomy in R Based on Pj Naughter’s Aa+ Library. https://github.com/helixcn/skycalc/.

Zhang, Jinlong.. 2016. Ephem: Creating Figures for Astronomical Almanac.https://github.com/helixcn/.

利用R提取海拔、温度和降水等信息

全球范围内的温度降水等信息， 常用于物种潜在分布区预测研究中。与此同时，在生态学研究中，常需要某些地点的温度降水等信息。 较早公开的较高精度的气候数据，包括Worldclim（http://www.worldclim.org/）。 2017年，另一套高精度数据CHELSA （Climatologies at high resolution for the earth’s land surface areas， http://chelsa-climate.org/） 公开发布。 CHELSA数据整合了局域尺度的降雨等信息， 部分气象因子比Worldclim更准确。 同一年，Worldclimate也发布了2.0版本。

已知某些地点的经纬度，要提取这些地点的温度降水等信息，用DIVA-GIS(http://www.diva-gis.org/ ) 和Worldclim数据很容易做到， 具体方法参见 http://blog.sciencenet.cn/blog-255662-679159.html

不过由于DIVA-GIS是依靠鼠标操作，很难批量化处理。 本文给出提取气候数据的R代码。相关操作主要基于Robert Hijmans的raster程序包， 更多详细信息请参考 （http://rspatial.org/index.html）

海拔数据下载
http://edcintl.cr.usgs.gov/downloads/sciweb1/shared/topo/downloads/GMTED/Grid_ZipFiles/mn30_grd.zip

CHELSA提供的19个Bioclimate 数据下载
https://www.wsl.ch/lud/chelsa/data/bioclim/integer/

各环境变量
Bio1 = Annual Mean Temperature 年均温

Bio2 = Mean Diurnal Range，即 Mean of monthly (max temp - min temp) 平均月温差

Bio3 = Isothermality 等温性

Bio4 = Temperature Seasonality 温度的季节性

Bio5 = Max Temperature of Warmest Month 最暖月最高温

Bio6 = Min Temperature of Coldest Month 最冷月最低温

Bio7 = Temperature Annual Range 气温年较差

Bio8 = Mean Temperature of Wettest Quarter 最湿季均温

Bio9 = Mean Temperature of Driest Quarter 最干季均温

Bio10 = Mean Temperature of Warmest Quarter 最暖月均温

Bio11 = Mean Temperature of Coldest Quarter 最冷月均温

Bio12 = Annual Precipitation 年降水

Bio13 = Precipitation of Wettest Month 最湿润月降水

Bio14 = Precipitation of Driest Month 最干旱月降水

Bio15 = Precipitation Seasonality 气温的季节性

Bio16 = Precipitation of Wettest Quarter 最湿润季降水

Bio17 = Precipitation of Driest Quarter 最干旱季降水

Bio18 = Precipitation of Warmest Quarter 最暖季降水

Bio19 = Precipitation of Coldest Quarter 最冷季降水

CHELSA气候图层， 选择整数还是浮点数图层？
可选择整数图层 integer ，这是因为整数图层占用空间较小， 同时也有相对较高的精度， 如温度可精确到0.1度。 注意，这里的温度是*10的值， 也就说，某地年均温为15.2度， 则提取出的数值为152。

30s数据能达到什么精度？
edge.length <- 63782pi1000/(36060*60)*30
edge.length
表明最大栅格边长约为928米。

计算对角线长度
edge.length *1.414

[1] 1311.686

所以如果两个地点相距不足1312m， 那么两个点有可能落在同一个栅格内， 那么提取的气候要素也就没有任何区别。

图层用的是什么投影？
所有图层都使用WGS1984投影， 用于提取数据的坐标， 也必须是WGS1984坐标

解压缩
下载的7z文件用 7-Zip（https://www.7-zip.org/ ）

R代码

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
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
setwd("C:/Users/jlzhang/Desktop/CHELSA")
library(raster)
## Loading required package: sp
library(openxlsx)
#### locations <- read.csv("https://simplemaps.com/static/data/world-cities/basic/simplemaps-worldcities-basic.csv", header = TRUE)
#### write.csv(locations, "locations.csv")
locations <- read.csv("locations.csv")
locations30 <- locations[sample(1:nrow(locations), 30),]
colnames(locations30)
##  [1] "X"          "city"       "city_ascii" "lat"        "lng"       
##  [6] "pop"        "country"    "iso2"       "iso3"       "province"
locations30
##         X            city      city_ascii        lat           lng
## 3081 3081       Bhagalpur       Bhagalpur  25.229996   86.98000321
## 6838 6838           Homer           Homer  59.642934 -151.54827970
## 7257 7257           Lahij           Lahij  13.058181   44.88376135
## 1214 1214            Dori            Dori  14.033997   -0.02799852
## 6716 6716         Douglas         Douglas  31.507778  -82.85068994
## 4313 4313      Bluefields      Bluefields  11.999977  -83.76494938
## 1743 1743          Fuyang          Fuyang  32.900407  115.82000000
## 6831 6831          Ambler          Ambler  67.086485 -157.85140910
## 6249 6249     Southampton     Southampton  50.900031   -1.39997685
## 4920 4920     Vladikavkaz     Vladikavkaz  43.050381   44.66997595
## 5423 5423     Al Mubarraz     Al Mubarraz  25.429054   49.56590450
## 5784 5784    Schaffhausen    Schaffhausen  47.706003    8.63299852
## 1304 1304         Creston         Creston  49.099960 -116.51669700
## 7110 7110        Guliston        Guliston  40.495731   68.79072587
## 5309 5309     Pervouralsk     Pervouralsk  56.910026   59.95503780
## 6658 6658        Woodward        Woodward  36.433421  -99.39769027
## 1405 1405        Cochrane        Cochrane  49.067017  -81.01656417
## 7264 7264          Sayhut          Sayhut  15.210504   51.24544023
## 187   187   Troll Station   Troll Station -72.016290    2.53332312
## 7034 7034      Cincinnati      Cincinnati  39.161885  -84.45692265
## 4395 4395          Sokoto          Sokoto  13.060015    5.24003129
## 3078 3078         Lucknow         Lucknow  26.855039   80.91499874
## 3877 3877       Tidjikdja       Tidjikdja  18.550016  -11.41660059
## 4130 4130           Rabat           Rabat  34.025299   -6.83613082
## 7000 7000      Montgomery      Montgomery  32.361602  -86.27918868
## 1990 1990            Bose            Bose  23.899716  106.61332680
## 6893 6893       Kalispell       Kalispell  48.197767 -114.31597860
## 3587 3587           Rudny           Rudny  52.952697   63.13003780
## 1030 1030        Mineiros        Mineiros -17.569536  -52.55998071
## 3114 3114 Padangsidempuan Padangsidempuan   1.388738   99.27336137
##            pop                  country iso2 iso3
## 3081  361548.0                    India   IN  IND
## 6838    5021.5 United States of America   US  USA
## 7257   44831.5                    Yemen   YE  YEM
## 1214   37806.0             Burkina Faso   BF  BFA
## 6716   12159.0 United States of America   US  USA
## 4313   40033.0                Nicaragua   NI  NIC
## 1743  170023.0                    China   CN  CHN
## 6831     258.0 United States of America   US  USA
## 6249  384417.0           United Kingdom   GB  GBR
## 4920  341000.0                   Russia   RU  RUS
## 5423  294682.0             Saudi Arabia   SA  SAU
## 5784   33863.0              Switzerland   CH  CHE
## 1304    4816.0                   Canada   CA  CAN
## 7110   74446.5               Uzbekistan   UZ  UZB
## 5309  127236.0                   Russia   RU  RUS
## 6658   12339.5 United States of America   US  USA
## 1405    4441.0                   Canada   CA  CAN
## 7264     189.0                    Yemen   YE  YEM
## 187       25.5               Antarctica   AQ  ATA
## 7034  971191.0 United States of America   US  USA
## 4395  648019.5                  Nigeria   NG  NGA
## 3078 2583505.5                    India   IN  IND
## 3877   19981.0               Mauritania   MR  MRT
## 4130 1680376.5                  Morocco   MA  MAR
## 7000  194491.5 United States of America   US  USA
## 1990  132942.5                    China   CN  CHN
## 6893   25040.0 United States of America   US  USA
## 3587  104235.5               Kazakhstan   KZ  KAZ
## 1030   28247.0                   Brazil   BR  BRA
## 3114  183721.5                Indonesia   ID  IDN
##                            province
## 3081                          Bihar
## 6838                         Alaska
## 7257                          Lahij
## 1214                          Séno
## 6716                        Georgia
## 4313                 Atlántico Sur
## 1743                          Anhui
## 6831                         Alaska
## 6249                    Southampton
## 4920                  North Ossetia
## 5423                  Ash Sharqiyah
## 5784                   Schaffhausen
## 1304               British Columbia
## 7110                       Sirdaryo
## 5309                     Sverdlovsk
## 6658                       Oklahoma
## 1405                        Ontario
## 7264                      Al Mahrah
## 187                                
## 7034                           Ohio
## 4395                         Sokoto
## 3078                  Uttar Pradesh
## 3877                         Tagant
## 4130 Rabat - Salé - Zemmour - Zaer
## 7000                        Alabama
## 1990                        Guangxi
## 6893                        Montana
## 3587                       Qostanay
## 1030                         Goiás
## 3114                 Sumatera Utara
coordinates(locations30) = c("lng","lat") # specify column names
### 年均温
bio1 <- raster("CHELSA_bio10_1.tif")
plot(bio1, main="Annual Mean Temperature", 
xlab="Longitude", ylab="Latitude", col=rev(heat.colors(100)))
img

res_bio1 <- extract(x = bio1, y = locations30)


### 海拔
elevation <- raster("mn30_grd")
plot(elevation, main="Altitude", 
xlab="Longitude", ylab="Latitude",col= terrain.colors(100))
img

res_elevation <- extract(x = elevation, y = locations30)

res <- data.frame(locations30, res_bio1, res_elevation)
res
##         X            city      city_ascii        lat           lng
## 3081 3081       Bhagalpur       Bhagalpur  25.229996   86.98000321
## 6838 6838           Homer           Homer  59.642934 -151.54827970
## 7257 7257           Lahij           Lahij  13.058181   44.88376135
## 1214 1214            Dori            Dori  14.033997   -0.02799852
## 6716 6716         Douglas         Douglas  31.507778  -82.85068994
## 4313 4313      Bluefields      Bluefields  11.999977  -83.76494938
## 1743 1743          Fuyang          Fuyang  32.900407  115.82000000
## 6831 6831          Ambler          Ambler  67.086485 -157.85140910
## 6249 6249     Southampton     Southampton  50.900031   -1.39997685
## 4920 4920     Vladikavkaz     Vladikavkaz  43.050381   44.66997595
## 5423 5423     Al Mubarraz     Al Mubarraz  25.429054   49.56590450
## 5784 5784    Schaffhausen    Schaffhausen  47.706003    8.63299852
## 1304 1304         Creston         Creston  49.099960 -116.51669700
## 7110 7110        Guliston        Guliston  40.495731   68.79072587
## 5309 5309     Pervouralsk     Pervouralsk  56.910026   59.95503780
## 6658 6658        Woodward        Woodward  36.433421  -99.39769027
## 1405 1405        Cochrane        Cochrane  49.067017  -81.01656417
## 7264 7264          Sayhut          Sayhut  15.210504   51.24544023
## 187   187   Troll Station   Troll Station -72.016290    2.53332312
## 7034 7034      Cincinnati      Cincinnati  39.161885  -84.45692265
## 4395 4395          Sokoto          Sokoto  13.060015    5.24003129
## 3078 3078         Lucknow         Lucknow  26.855039   80.91499874
## 3877 3877       Tidjikdja       Tidjikdja  18.550016  -11.41660059
## 4130 4130           Rabat           Rabat  34.025299   -6.83613082
## 7000 7000      Montgomery      Montgomery  32.361602  -86.27918868
## 1990 1990            Bose            Bose  23.899716  106.61332680
## 6893 6893       Kalispell       Kalispell  48.197767 -114.31597860
## 3587 3587           Rudny           Rudny  52.952697   63.13003780
## 1030 1030        Mineiros        Mineiros -17.569536  -52.55998071
## 3114 3114 Padangsidempuan Padangsidempuan   1.388738   99.27336137
##            pop                  country iso2 iso3
## 3081  361548.0                    India   IN  IND
## 6838    5021.5 United States of America   US  USA
## 7257   44831.5                    Yemen   YE  YEM
## 1214   37806.0             Burkina Faso   BF  BFA
## 6716   12159.0 United States of America   US  USA
## 4313   40033.0                Nicaragua   NI  NIC
## 1743  170023.0                    China   CN  CHN
## 6831     258.0 United States of America   US  USA
## 6249  384417.0           United Kingdom   GB  GBR
## 4920  341000.0                   Russia   RU  RUS
## 5423  294682.0             Saudi Arabia   SA  SAU
## 5784   33863.0              Switzerland   CH  CHE
## 1304    4816.0                   Canada   CA  CAN
## 7110   74446.5               Uzbekistan   UZ  UZB
## 5309  127236.0                   Russia   RU  RUS
## 6658   12339.5 United States of America   US  USA
## 1405    4441.0                   Canada   CA  CAN
## 7264     189.0                    Yemen   YE  YEM
## 187       25.5               Antarctica   AQ  ATA
## 7034  971191.0 United States of America   US  USA
## 4395  648019.5                  Nigeria   NG  NGA
## 3078 2583505.5                    India   IN  IND
## 3877   19981.0               Mauritania   MR  MRT
## 4130 1680376.5                  Morocco   MA  MAR
## 7000  194491.5 United States of America   US  USA
## 1990  132942.5                    China   CN  CHN
## 6893   25040.0 United States of America   US  USA
## 3587  104235.5               Kazakhstan   KZ  KAZ
## 1030   28247.0                   Brazil   BR  BRA
## 3114  183721.5                Indonesia   ID  IDN
##                            province optional res_bio1 res_elevation
## 3081                          Bihar     TRUE      252            44
## 6838                         Alaska     TRUE       46            56
## 7257                          Lahij     TRUE      294           128
## 1214                          Séno     TRUE      290           281
## 6716                        Georgia     TRUE      185            81
## 4313                 Atlántico Sur     TRUE      266             5
## 1743                          Anhui     TRUE      158            34
## 6831                         Alaska     TRUE      -46            34
## 6249                    Southampton     TRUE      109             6
## 4920                  North Ossetia     TRUE      103           663
## 5423                  Ash Sharqiyah     TRUE      273           155
## 5784                   Schaffhausen     TRUE      100           446
## 1304               British Columbia     TRUE       88           622
## 7110                       Sirdaryo     TRUE      163           273
## 5309                     Sverdlovsk     TRUE       20           329
## 6658                       Oklahoma     TRUE      146           592
## 1405                        Ontario     TRUE       16           276
## 7264                      Al Mahrah     TRUE      283            10
## 187                                     TRUE     -212          1242
## 7034                           Ohio     TRUE      126           193
## 4395                         Sokoto     TRUE      285           291
## 3078                  Uttar Pradesh     TRUE      251           128
## 3877                         Tagant     TRUE      283           409
## 4130 Rabat - Salé - Zemmour - Zaer     TRUE      193            22
## 7000                        Alabama     TRUE      183            88
## 1990                        Guangxi     TRUE      223           137
## 6893                        Montana     TRUE       78           903
## 3587                       Qostanay     TRUE       39           145
## 1030                         Goiás     TRUE      229           769
## 3114                 Sumatera Utara     TRUE      250           325
write.xlsx(res, "extracted.xlsx")

在R中分析物种和系统发育多样性

背景

本课程中， 我们将学习如何分析阿尔伯塔草原的一组植物群落数据。数据中包括阿尔伯塔几个地点草原植物盖度、功能性状和系统发育关系等，

数据详细介绍参见以下论文: S.W. Kembel and J.F. Cahill, Jr. 2011. Independent evolution of leaf and root traits within and among temperate grassland plant communities. PLoS ONE 6(6): e19992. (doi:10.1371/journal.pone.0019992).

怎样使用课程材料

本课程覆盖了导入数据，在R中分析生物多样性的全部内容。如果要学习全部课程， 可以从R入门部分开始。如果已经有R的基础， 则可以直接学习本讲义。本课程需要用到picante程序包，以及数据映像（即.Rdata），后者包括练习用的全部数据，有了这些数据，运行后面的命令就可直接得到结果。

生物多样性数据导入R

首先， 查看一下我们需要的程序包是否能加载。需要用的程序包有ape, picante, vegan. 由于picante是依赖另外两个程序包的，所以加载picante的时候，另外两个程序包也会自动加载。

设定工作路径，请注意Windows和MacOS以及Linux操作系统上，路径斜杠的方向

1
2
3
#setwd("/Users/jinlong/Dropbox/04\ to\ review/kembel/")
setwd("C:\\Users\\jlzhang\\Dropbox\\04 to review\\kembel")
library(picante)

群落数据

群落数据包括每个地点，或者每块样地物种的个体数，或者相对多度等。在本课程的数据中，多度数据是每个不同生境20m*20m样方中植物的相对盖度。

群落数据的格式为data.frame, 行表示样地，列表示物种。 练习数据已经包括comm数据了， 所以直接用就行了。 如果数据是保存在csv文件中， 则需要用read.csv函数读取，格式如下：

1
2
3
4
# 每一行代表一块样地，因此行名就是样地编号，
# 第一行作为各列名称，所以这里设定header = TRUE。

comm <- read.csv("grassland.community.csv", header = TRUE, row.names = 1)

用以上方法将数据读取到R中，注意行名和列名并不是数据本身，而是各行列的名称。行列有了名称， 在后续就更容易操作和展示。群落数据的各行列有了名称，才方便和其他类型的数据关联。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# 先查看comm的数据类型
class(comm)

# 查看该data.frame是由多少行列组成的。(rows x columns)
dim(comm)

# 查看行名
rownames(comm)

# 查看前六列的名
head(colnames(comm))

# 提取数据的子集，1到5行，1到5列数据
comm[1:5, 1:5]

本数据中， 多度其实是样地中某种植物的盖度。多元统计中的很多方法对于多度很敏感，因此可以将绝对多度转换为相对多度。用vegan中的函数即可完成。 

1
2
3
4
5
6
7
8
9
10
11
12
# 查看每个样方中盖度之合， 如果数据为个体数， 
# 则表示每个样方中所有种的总个体数
apply(comm, 1, sum)

# 除以每个样方的盖度和， 将盖度转换为相对盖度
comm <- decostand(comm, method = "total")

# 每个样方的总盖度
apply(comm, 1, sum)

# 查看转换后的数据
comm[1:5, 1:5]

性状数据

trait包括每个种叶和根的数据。读取到R中的方法与群落数据一样，但是在性状数据中，每个种是一行，每个性状是一列。

1
2
3
4
5
6
7
8
9
10
11
12
13
traits <- read.csv("species.traits.csv", header = TRUE, row.names = 1)

# 查看性状数据的前6列
head(traits)

# 绘制两两相关图
pairs(traits)

# 由于一些变量显著偏向一侧，因此进行对数转换
traits <- log10(traits)

# 查看转换后的数据
pairs(traits)

样地的总体情况Metadata

主要包括样地的生境，采集地点，环境数据，如坡度、干湿度等

1
2
3
4
metadata <- read.csv("plot.metadata.csv", header = TRUE, row.names = 1)

# 查看前6行
head(metadata)

进化树（或称为系统树）

newick格式的进化树，用read.tree读取. Nexus格式的进化树，用read.nexus读取。

1
2
3
phy <- read.tree("grassland.phylogeny.newick")
class(phy)
phy

读取到R中的进化树为phylo为格式。 phylo格式的详细说明参见ape的主页(http://ape.mpl.ird.fr/)。 phylo对象的本质是list，包括进化树末端分类单元名称（tip label），枝长（edge length）等等， ape内置的相关函数可以针对phylo这种数据格式打印， 汇总， 绘图等。

1
2
3
4
5
6
7
8
9
10
11
# 列出phy对象的所有名称， 这些名称是list内部对象的名
names(phy)

# 查看末端分类单元名称， 取前5个
phy$tip.label[1:5]

# 进化树有多少末端分类单元？
Ntip(phy)

# 绘制进化树，用cex参数调整物种名字体大小
plot(phy, cex = 0.5)

数据整理（cleaning）及匹配

本分析用到的数据包括群落、性状、进化树以及metadata。

1
ls()

示例中的数据，因为已经经过了认真整理，样地和样品中的物种名能够完全对应。但是很多情况下，物种名并不能完全匹配。例如， 样地数据可能只有一部分种出现在进化树中，有些种可能有性状数据，但未出现在进化树中。一些分析中， R假设物种在群落数据和进化树中出现的顺序是一致的。但是有可能出现输入错误等。 以上问题都要找出来。picante程序包中有几个函数来检测不同数据中物种名是否相同。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# 检查物种是否匹配
combined <- match.phylo.comm(phy, comm)
# 返回结果是一个包含 phy和com的list，
# 未在两个数据集中都出现的物种都已经去除，
# 并且顺序都排好了。
# 这里用combined中的数据，替换原始数据。

phy <- combined$phy
comm <- combined$comm
# 性状数据也需要进行同样的处理
combined <- match.phylo.data(phy, traits)

# 用结果中检查过的数据，替换掉原来的数据。
# 由于之前的返回结果是list， 这里用$调取phy和data
phy <- combined$phy
traits <- combined$data

#再检查一下metadata与群落数据出现的顺序是否一致？
all.equal(rownames(comm), rownames(metadata))

# 若顺序不同， metadata数据中样方出现的顺序就应该按照群落顺序排序
metadata <- metadata[rownames(comm), ]

# 这样，数据预处理就完成了。 后面即可开始各种生物多样性分析。

生物多样性数据的可视化及汇总

物种丰富度和多样性

每种生境上的物种数是否相同？

1
2
3
4
5
6
7
8
# 比较 fescue 和 mixedgrass 生境中物种数的差异
boxplot(specnumber(comm) ~ metadata$habitat, ylab = "# of species")

# 用t检验比较不同生境中物种数是否存在差异
t.test(specnumber(comm) ~ metadata$habitat)

# 物种累计曲线，用于查看样方总面积是否足够大
plot(specaccum(comm), xlab = "# of samples", ylab = "# of species")

群落多元分析

不同样地植物群落物种组成是怎样变化的？生境类型和环境变量与群落植物组成有什么关系？

要研究这些问题， 需要用到多元排序方法。 相关函数在vegan程序包中，相关的函数的帮助文件和使用指南也非常详尽。 也可以参考 Borcard等人编写的 Numerical Ecology in R， 这本书在2018年出的第二版。

等级聚类

通过聚类，物种组成相近的群落能聚合到一起。 先计算样方之间的Bray-Curtis距离，Bray-Curtis不但考虑了物种组成的差异， 同时考虑了每个种的多度。之后，用聚合等级聚类（agglomerative hierachical clustering algorithm）算法即可实现聚类。 

1
2
3
4
5
6
7
8
# 计算样地之间的 Bray-Curtis距离 
comm.bc.dist <- vegdist(comm, method = "bray")

# 用UPGMA方法聚类
comm.bc.clust <- hclust(comm.bc.dist, method = "average")

# 绘制聚类图
plot(comm.bc.clust, ylab = "Bray-Curtis dissimilarity")

两种类型的草地， 由于有不同的物种组成，分别聚成两类

排序

R中的排序方法比较多，这里我们采用NMDS （non-metric multidimensional scaling）方法。 NMDS的方法是基于距离矩阵的。

1
2
3
4
5
# metaMDS函数会自动进行数据转换，并检查结果是否已经收敛。
comm.bc.mds <- metaMDS(comm, dist = "bray")

# 绘制stressplot，检查排序结果
stressplot(comm.bc.mds)

排序结果可通过多种方式展示

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
# 在Biplot中显示样地
ordiplot(comm.bc.mds, display = "sites", type = "text")

# 添加物种名
ordipointlabel(comm.bc.mds)

# ordiplot可以只打开绘图设备和画布，不绘制任何内容
mds.fig <- ordiplot(comm.bc.mds, type = "none")

# 只添加样地，颜色为样地类型，pch这里设定的是点的形状。
# 试试 plot(1:30, pch = 1:30)
points(mds.fig, "sites", pch = 19, col = "green", 
       select = metadata$habitat == "Fescue")
points(mds.fig, "sites", pch = 19, col = "blue", 
       select = metadata$habitat == "Mixedgrass")

# 添加生境置信椭圆
ordiellipse(comm.bc.mds, metadata$habitat, 
            conf = 0.95, label = TRUE)

# 将聚类结果在NMDS排序图中显示
ordicluster(comm.bc.mds, comm.bc.clust, col = "gray")

也可以用ordisurf函数显示物种多度

1
2
3
4
# 显示球葵属Sphaeralcea 的多度. cex参数可以调整点的大小
ordisurf(comm.bc.mds, comm[, "Sphaeralcea_coccinea"], 
         bubble = TRUE, main = "Sphaeralcea coccinea abundance",
         cex = 3)

添加环境信息和性状数据到排序中

环境变量与排序轴之间是怎样的关系？这里使用envfit函数拟合

1
2
ordiplot(comm.bc.mds)
plot(envfit(comm.bc.mds, metadata[, 3:6]))

用vegan中的cca或者rda等排序方法也是可以的。cca或者rda直接考虑了物种组成和多度， 并不用先计算样方之间的距离。

性状与进化树

系统发育信号

最近，系统发育保守性受到很多人关注。群落系统发育结构分析就假设物种的性状是系统发育保守的。

系统发育信号是度量不同物种的性状在进化树上相似程度的指数。Blomberg K统计量就是比较性状的观察值与布朗运动模型预测值的统计量 （Blomberg et al. 2003）。K 值接近1表明性状进化过程接近布朗运动， 表明有一定程度的系统发育信号或者呈现一定的保守性。K接近于0表明性状进化倾向于随机，K>1 表明性状保守。

K的显著性检验可以用以下方法： 将进化树上的物种名随机打乱，计算每次打乱物种名后，性状的系统发育独立差（Phylogenetic independent contrast）的方差, 在多次打乱物种名生成零分布后，真实进化树的系统发育独立差方差与之进行比较。相关检验在Kcalc，phylosignal和multiPhylosignal 函数中可实现。

下面用练习数据计算系统发育信号

1
2
3
4
5
6
# 可以用Kcalc依次计算traits数据框的所有列
apply(traits, 2, Kcalc, phy)

# multiPhylosignal函数可以检验多列的P值，
# 但该函数需要的进化树为严格的二分叉树
multiPhylosignal(traits, multi2di(phy))

结果中包括K和PIC.variance.P, 后者表示性状的非随机程度。大部分变量都完全随机表现出更强的系统发育信号。

性状进化的可视化

1
2
3
4
plot(phy, direction = "up", show.tip.label = FALSE, 
     show.node.label = TRUE, cex = 0.7)
tiplabels(pch = 19, col = "black", 
          cex = 3 * (traits[, "LeafArea"]/max(traits[,"LeafArea"])))

性状相关性的系统发育分析

若功能性状表现出较强的系统发育信号，就违背了数据完全独立的假设。此时可以在gls中设定corBrownian 关联矩阵，考虑物种之间的系统发育关系。（注：在分析过程中， 可以使用caper的pgls回归。）

其中gls是一般最小二乘(Generalised least squares)法的简称，方法类似ANOVA 或线性模型。无序类别变量以及连续变量都可以通过这种方式检验。

下面检验检测一下， 考虑物种之间系统发育关系时， specific root length (SRL) 和 root tissue density 之间的关系。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# 普通的gls模型， 不考虑系统发育关系
root.gls <- gls(RootTissueDens ~ SRL, data = traits)
anova(root.gls)

# 考虑系统发育关系
root.pgls <- gls(RootTissueDens ~ SRL, 
                 correlation = corBrownian(value = 1, phy),
                 data = traits)
anova(root.pgls)
# 绘图
plot(RootTissueDens ~ SRL, data = traits, 
     xlab = "SRL (specific root length)",
     ylab = "Root tissue density")
abline(coef(root.gls), lwd = 2, col = "black")
abline(coef(root.pgls), lwd = 2, col = "red")
legend("bottomleft", 
       legend = c("GLS fit", "Phylogenetic GLS fit"), 
       lwd = 2, col = c("black", "red"))

当不考虑系统发育关系时，SRL和root tissue density的关系很弱。考虑系统发育信号后，两者的相关性就明显了。

系统发育与性状多样性

系统发育多样性

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# 计算Faith's PD
comm.pd <- pd(comm, phy)
head(comm.pd)

# 不同生境的Faith's PD
boxplot(comm.pd$PD ~ metadata$habitat, 
        xlab = "Habitat", ylab = "Faith's PD")

# 不同生境PD是否相同
t.test(comm.pd$PD ~ metadata$habitat)

# 查看PD和物种丰富度的关系
plot(comm.pd$PD ~ comm.pd$SR, 
     xlab = "Species richness", ylab = "Faith's PD")

picante可以用来计算
(MPD), (MNTD), (SES_{MPD}) and (SES_{MNTD})
具体定义， 请参见Phylocom的说明书。

群落系统发育最重要指数之一就是标准效应值。公式如下
(SES_{metric} = \frac{ Metric_{observed} - mean(Metric_{null}) }{sd(Metric_{null})})

一般包括：
(SES_{MPD}) 等于 -NRI， (SES_{MNTD}) 等于 -NTI

picante中提供了不同的零模型，用来进行随机化，以计算SES。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# 系统发育距离矩阵
phy.dist <- cophenetic(phy)

# ses.mpd
comm.sesmpd <- ses.mpd(comm, phy.dist, null.model = "richness", 
                       abundance.weighted = FALSE, 
                       runs = 999)
head(comm.sesmpd)

# ses.mntd
comm.sesmntd <- ses.mntd(comm, phy.dist, null.model = "richness",
                         abundance.weighted = FALSE,
                         runs = 999)
head(comm.sesmntd)

输出结果
ntaxa - 物种丰富度
mpd.obs - 群落中mpd观察值
mpd.rand.mean - 零群落mpd平均值
mpd.rand.sd - 零模型mpd的标准差
mpd.obs.rank - 观察值在零模型中的顺序rank（秩）
mpd.obs.z - mpd标准效应指数 (等于 -NRI)
mpd.obs.p - P-value (quantile) of observed mpd vs. null communities (= mpd.obs.rank / runs + 1)
runs - 随机化次数

SES.mpd>0，且mpd.obs.p>0.95时，表示群落系统发育均匀。 SES.mpd<0，且mpd.obs.p<0.05时，表明与零模型相比系统发育聚集。

一般认为， MPD对对系统发育整体的聚集性或均匀性较为敏感，MNTD对靠近进化树末端的均匀性和聚集性更为敏感。

1
2
3
4
5
6
7
8
9
10
11
# 不同生境的ses.mpd
plot(comm.sesmpd$mpd.obs.z ~ metadata$habitat, 
     xlab = "Habitat", ylab = "SES(MPD)")
abline(h = 0, col = "gray")
t.test(comm.sesmpd$mpd.obs.z ~ metadata$habitat)

# 不同生境的ses.mntd
plot(comm.sesmntd$mntd.obs.z ~ metadata$habitat, 
     xlab = "Habitat", ylab = "SES(MNTD)")
abline(h = 0, col = "gray")
t.test(comm.sesmntd$mntd.obs.z ~ metadata$habitat)
1
2
3
4
5
6
# 某fescue群落中的物种
plot(phy, show.tip.label = FALSE, 
     main = "Fescue community fes-K-11")
tiplabels(tip = which(phy$tip.label %in% 
                        colnames(comm)[comm["fes-K-11", ] > 0]), 
          pch = 19)

‘mix-H-23’ Mixedgrass 群落中包含了一系列系统发育聚集的种。

1
2
3
4
5
6
# 某mixedgrass群落中出现的个体
plot(phy, show.tip.label = FALSE, 
     main = "Fescue community mix-H-23")
tiplabels(tip = which(phy$tip.label %in% 
                        colnames(comm)[comm["mix-H-23", ] >
                                         0]), pch = 19)

性状多样性

性状的分析方法与进化树分析的方法类似，一般都是先进行标准化，再计算欧式距离，再计算功能性状相似性，计算MPD或者MNTD

1
2
3
4
5
6
7
# 性状先用scale标准化， 再用dist计算欧式距离计算
trait.dist <- as.matrix(dist(scale(traits), method = "euclidean"))
comm.sesmpd.traits <- ses.mpd(comm, trait.dist, null.model = "richness", 
                              abundance.weighted = FALSE, runs = 999)
plot(comm.sesmpd.traits$mpd.obs.z ~ metadata$habitat, 
     xlab = "Habitat", ylab = "Trait SES(MPD)")
abline(h = 0, col = "gray")

vegan中的treedive函数与picante中的pd类似。

系统发育beta多样性

unifrac函数和phylosor函数，都是计算样方之间的Faith’s PD。 comdist和comdistnt的计算类似于MPD和MNTD，但是是针对样方之间的距离。计算过程中可以考虑多度，由于之前我们计算的Bray-Curtis距离已经考虑了多度。为了方便比较，后续的计算也考虑多度， 即abundance.weighted = TRUE

1
2
3
4
5
6
7
8
9
10
11
12
13
# 系统发育beta多样性
comm.mntd.dist <- comdistnt(comm, phy.dist, 
                            abundance.weighted = TRUE)
# 功能性状beta多样性
comm.mntd.traits.dist <- comdistnt(comm, trait.dist, 
                                   abundance.weighted = TRUE)

# 计算两个距离矩阵的相关性，用mantel检验
# 此处计算Bray-Curtis距离和mntd距离的相关性
mantel(comm.bc.dist, comm.mntd.dist)

# 此处计算Bray-Curtis距离和功能性状距离的相关性
mantel(comm.bc.dist, comm.mntd.traits.dist)

系统发育和功能性状排序

之前，我们用Bray-Curtis进行了NDMS排序，实际上，系统发育MNTD以及功能性状MNTD都可以用类似的方法排序， 以展示样地之间的关系

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
# 由于只有样地之间的距离，这里用monoMDS. 
# 具体请认真阅读 monoMDS的说明书。

# 系统发育距离
comm.mntd.mds <- monoMDS(comm.mntd.dist)
mds.fig <- ordiplot(comm.mntd.mds, type = "none")
points(mds.fig, "sites", pch = 19, col = "green", 
       select = metadata$habitat == "Fescue")
points(mds.fig, "sites", pch = 19, col = "blue", 
       select = metadata$habitat == "Mixedgrass")
ordiellipse(comm.mntd.mds, metadata$habitat, 
            conf = 0.95, label = TRUE)

# 功能性状
comm.mntd.traits.mds <- monoMDS(comm.mntd.traits.dist)
mds.fig <- ordiplot(comm.mntd.traits.mds, type = "none")
points(mds.fig, "sites", pch = 19, col = "green", 
       select = metadata$habitat == "Fescue")
points(mds.fig, "sites", pch = 19, col = "blue", 
       select = metadata$habitat == "Mixedgrass")
ordiellipse(comm.mntd.traits.mds, metadata$habitat, 
            conf = 0.95, label = TRUE)

无论是对物种beta多样性， 系统发育beta多样性还是功能性状beta多样性进行排序方式，两种生境类型的差别都比较大。

##解释beta多样性

为了解释哪些因子对样方之间的系统发育关系有影响， 可以用permutational MANOVA （adonis）计算，R代码如下。输出结果类似ANOVA方差分析。

以下分别检验生境对于Bray-Curtis距离、系统发育距离、功能性状距离的解释能力。

1
2
3
4
5
6
7
8
# Bray-Curtis距离
adonis(comm.bc.dist ~ habitat, data = metadata)

# 系统发育多样性距离
adonis(comm.mntd.dist ~ habitat, data = metadata)

# 功能性状距离
adonis(comm.mntd.traits.dist ~ habitat, data = metadata)

清晨

冬季的清晨六点， 华北平原上的小村还沉浸在黑暗里， 东边的天空一丝鱼肚白也没有。

从家里出发到等车的地方，只需要五分钟，走在路上， 却觉得如此漫长。周围漆黑一片，夜空仿佛被黑幕笼罩着，什么都看不清，星星也没有。一切沉浸在冬天那安静的雾气里，只能偶尔听到公鸡报晓。早起人们的走路声打破了拂晓的沉寂。

等车的人裹着厚厚的棉衣，穿着粗笨的棉鞋，头上或是包着围巾或是戴着帽子。零下七八度，已经足够严寒。 不得不把双手插进口袋或袖口，站在穿村而过的乡间公路边， 吸进带着白霜的寒冷空气，身体也不知不觉瑟瑟发抖。 呼出的“哈气”只有偶尔在汽车经过时才能看到，“哈气”在呼出的瞬间向上飘散，消失在眉宇间，湿润了睫毛和头发，很快又结成霜。

这就是二十年前一天清晨，等车时的一切。

车来了。

一辆白色的破旧中巴车缓缓驶来。狭窄的乡间公路似乎刚刚容得下一辆中巴。公路并不平， 车跟着公路的坑坑洼洼一直在喘。 车开过来了，车灯照得人睁不开眼。

上车后拣一个靠窗的位子坐下。好不容易等到了车， 车里却汽油味弥漫。 还好我不晕车，不然短短的旅途又不知道长了多少倍。窗玻璃上结满一层白霜，阻挡了视线，白霜用窗帘都擦不掉，如果真要看，只能用手的温度化开一小片，向往望望。 窗外的世界还沉浸在黑暗里， 偶尔有一两盏灯光， 不看也罢。

车窗玻璃由深蓝逐渐变成了浅蓝。 最后外面天光大亮，窗玻璃后面仿佛也在演着电影。天亮了，窗玻璃上的霜花也越来越清楚。汽车前方的玻璃上，远远感觉只是薄薄一层雾气， 靠近看却是霜。车窗的四角霜多些，多到看不到车外的景物。 距离自己近的玻璃上霜总是厚些，霜花又叫霜雪（雪读轻声）。霜雪里面有各种想要的图案：有山，有河，有农田，有树，有云彩，有大雁，有房子，有背着书包上学的小学生，有汽车，有大桥……有一切你能想到的东西——结霜的玻璃，总是让人忍不住冥想。

天边红起来了。 这就是冬天的日出，华北平原上一个普通得不能再普通的日出。那时候，一切好像都还没有准备好：炊烟，薄雾都还没有，小河也在冰封，树叶也都落了，冬天的大地透着一望无际的荒芜，小河沟里的芦苇早已枯黄，仿佛告诉你，这一切真的还没准备好迎接新的一天。然而太阳还是升起来了，如天边的一盏红灯，发出温暖而柔和的光。 太阳总是在你不经意的时候出来，就在最红最红的那片云彩下面。

太阳一出来，一切就有了温暖。

太阳出来，大地上绵延千里的白霜仿佛一瞬间就散了。大地暖了，风也暖了。公路两旁的树木，把影子投在公路上，中巴车的影子就在树木的影子上缓缓驶过。 玻璃上的霜花，不知不觉间就已经消融： 那一切，或者平凡，或者普通，或者壮丽辉煌的事物，还有一个人全部的想象，都在那一刻戛然而止。窗外的路和景透进来， 一切又变得俗不可耐。

新的一天了。

车窗外，有人在上学路上，有人骑车去赶集，有人在吃早点……一切的一切，就这样，在这个冬日的清晨，以这种方式开始。 从那时候开始，每个人的生活都以同样的方式开始了几千次。

一切有情皆过往。那一个寒冬离开了，那一个清晨的雾气散去了，那一块窗玻璃上的冰雪消融了。

那一天的风，吹到哪里了？那一天的目的地， 有没有到达？ 那一天里的梦，有没有醒来？

2018年5月25日凌晨

香港大埔

用R绘制采样点地图——带指北针、比例尺、图例

用R绘制采样点地图——带指北针、比例尺、图例

已有 4286 次阅读 2018-4-28 01:03 |系统分类:科研笔记 推荐到群组

本例子中的代码已经不推荐使用，tmap版本的地图示例参见 http://blog.sciencenet.cn/blog-255662-1144743.html

如果所示： 用R绘制的采样点地图——带指北针、比例尺、图例等

更多请参考 https://news-at.zhihu.com/story/9275383

源代码和数据下载 用R绘制采样点地图.zip

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
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
setwd("/Users/jinlong/world map")

library(maps)     ## 提供绘制世界地图的函数

library(mapdata)  ## 提供世界地图

library(openxlsx) ## 读取xlsx文件

library(ggplot2)  ## 绘图核心函数

library(ggsn)     ## 提供指北针

library(legendMap) ## 提供比例尺

\### devtools::install_github("3wen/legendMap")

library(maptools) ## 提供readShapePoly函数， 读取中国多边形数据



\### Download https://simplemaps.com/static/data/world-cities/basic/simplemaps-worldcities-basic.csv



china0   <-readShapePoly("bou1_4p") ## 读取中国地图

china_df <-fortify(china0) ## 将list转换为ggplot可以使用的dataframe



\##读取数据https://simplemaps.com/static/data/world-cities/basic/simplemaps-worldcities-basic.csv

cities <-read.csv("simplemaps-worldcities-basic.csv")

cities <-cities[sample(1:nrow(cities),size =200), ] ## 元数据有7000条， 这里随机选出100条



pdf(file ="sampling sites location.pdf",width =11.5,height =5.8)

\### tiff(filename = "sampling sites location.tiff",

\###      width = 7200, height = 3600, res = 600,

\###      compression = "lzw")  ## 如果输出tiff文件，600dpi，宽度7200像素，高度3600像素

rrr <-ggplot()+

borders(colour='darkgrey')+                 ##添加世界地图

geom_polygon(data =china_df,

 aes(x =long, y =lat, group =group), colour ="darkgrey",

 fill ="white")+                          ##添加中国地图

geom_point(data =cities,

 aes(x=lng,y=lat,colour=pop),size =2)+ ##点的颜色代表人口

 theme_bw()+                               ##简约风格主题

 ggtitle("Population of World Cities")+    ##标题

 scale_bar(lon =100,lat =-60,

​           distance_lon =1000,distance_lat =200,

​           distance_legend =500,dist_unit ="km",

​           arrow_length =0,arrow_distance =0,

​           arrow_north_size =0)+    ##比例尺

 xlab("Longitude (Degree)")+         ##横坐标

 ylab("Latitude (Degree)")+          ##纵坐标

 scale_colour_gradientn(colours=c("green","yellow","orange","red"))


north2(rrr,.82,.88)                   ## 数字是指定指北针的位置

dev.off()

大样地海拔数据的克里格插值与地形图绘制

最近收到几位同学求助， 咨询如何通过全站仪测得的海拔数据绘制大样地的地形图以及计算地形因子，本文作简要说明。本文中，“样方”是大样地中边长为20m的样方。前面的博文 （http://blog.sciencenet.cn/blog-255662-1081363.html ）简要介绍了大样地中样方的命名方法以及物种调查数据的整理过程，可供参考。

大样地研究除了收集物种分布数据，还需要获得每个样方对应的环境因子才能进行进一步分析。环境因子又包括地形变量和土壤因子变量，这里只说地形数据。地形数据一般包括某尺度下（如5m, 10m, 20m, 50m, 100m等）每个样方中心点的海拔、坡度、坡向以及凹凸度等。其中坡度、坡向和凹凸度又都是通过样方中心点海拔矩阵计算得出的。

每个样方中心点的海拔有两种获取方式：（1） 用全站仪直接测量，但是这种情况极少见；（2）对全站仪测得的海拔数据进行克里格插值之后，再在海拔模型上从中心点坐标提取。

为什么不全都用全站仪在每个20m样方的中心点直接测量呢？ 这是由于全站仪在打点时一般都是测量各样方顶点的海拔。有些样方的顶点正好位于深沟内或巨石上，难以获得，此时只能在附近的点测量。因此大多数情况下，全站仪所测得的数据并总不出现在样方顶点。所以第二种处理方式就很常见了。

有些研究中，样方的海拔没有经过克里格插值之后再提取中心点，而是直接用四个顶点海拔的平均值，这也是一种近似的处理方法。仍然是上述原因，相应条件不容易满足。

下面以西双版纳补蚌大样地的海拔数据为例， 介绍海拔数据读取，克里格插值，并提取20m中心点的海拔。获得样方中心点海拔之后，即可计算出坡度、坡向和凹凸度等（http://blog.sciencenet.cn/blog-255662-690600.html ）

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
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
setwd("C:/Users/jlzhang/Desktop/bubeng plot")

library(gstat)

library(sp)

library(raster)

library(openxlsx)

library(lattice)

 

#### 边长为20m米的点阵

bubeng.grid <- expand.grid(x = seq(0, 400, 20), y = seq(0, 500, 20))

coordinates(bubeng.grid) <- ~ x + y

gridded(bubeng.grid) <- TRUE

 

## 数据来源 https://doi.org/10.1371/journal.pone.0108450.s008

bubeng.topo20160318 <- read.xlsx("bubeng_topo.xlsx")

 

## 指定横纵坐标

coordinates(bubeng.topo20160318) = ~ x + y

 

## 变异函数

bubeng_variogram <- variogram(alt ~ 1, data = bubeng.topo20160318)

plot(bubeng_variogram)

 

### 拟合高斯模型， 初始值选择参见 （How Kriging works—Help）

bubeng_variogram.fit = fit.variogram(bubeng_variogram, model = vgm(psill = 10, "Gau", range = 200, nugget = .1))

plot(bubeng_variogram, bubeng_variogram.fit)

 

### 进行克里格插值

bubeng.krig <- krige(alt ~ 1, bubeng.topo20160318, bubeng.grid, model = bubeng_variogram.fit)

 

### 绘图

levelplot(bubeng.krig$var1.pred ~ bubeng.krig$x+bubeng.krig$y, cuts = 50,

          main="Ordinary kriging predictions",contour = TRUE,

          labels = TRUE, col.regions = terrain.colors(255),  pretty = TRUE)

 

## 三维图

wireframe(bubeng.krig$var1.pred ~ bubeng.krig$x+bubeng.krig$y,

          main="Ordinary kriging predictions",

          labels = TRUE, col.regions = terrain.colors(255),  pretty = TRUE, drape = TRUE)

 

## 转换为raster

dat <- data.frame(bubeng.krig$var1.pred, bubeng.krig$x, bubeng.krig$y)

 

### Convert the object x to a raster object.

rdat <- raster(bubeng.krig)

 

### plot the raster

plot(rdat, col = terrain.colors(255))

contour(rdat, add = TRUE, nlevel = 15)

 

### 获得每个20m样方的中心点

x10 <- seq(10, 400-10, by = 20)

y10 <- seq(10, 500-10, by = 20)

grids10 <- expand.grid(x10, y10)

plot(Var2 ~ Var1, data = grids10)

coordinates(grids10) <- c(1, 2)

 

###从raster图层中提取海拔

alt_center <- extract(x = rdat, y = grids10)

 

### 20m样方的中心点海拔数据

res_alt_center <- data.frame(grids10, alt_center)

参考文献

How Kriging works—Help | ArcGIS for Desktop. (2018). Available at: http://desktop.arcgis.com/en/arcmap/10.3/tools/3d-analyst-toolbox/how-kriging-works.htm. Last accessed 24 April 2018.

Hu, Y.-H., Kitching, R.L., Lan, G.-Y., Zhang, J.-L., Sha, L.-Q. & Cao, M. (2014). Size-Class Effect Contributes to Tree Species Assembly through Influencing Dispersal in Tropical Forests. PLoS ONE, 9, e108450.

网络时光与网络时光机

如果你学过代码管理， 就一定知道Subversions， 知道Tortoise SVN， 知道git。 这些软件能记录你对代码的修改， 每做出一定的修改后， 通常是修改了一些Bug或者增加了一些新的功能， 就保存一次修改意见， 这样就增加了一个版本。 好比攀岩， 每向上一步， 就加上一个锚定点， 确保不会从山崖上滑落。如果进行了错误的更改， 那么最好的办法就是回到之前的版本， SVN以及git都有相应的机制帮你回到之前某一个版本。所以SVN或者git称得上是软件源代码的时光机， 不但保证你稳步向前， 还能带你回到过去。

那么已经关闭的网站， 能否回到过去呢？ 答案可能是否定的， 如果网站没有进行很好地备份，数据损坏后， 再回去， 几乎是不可能的。不过有一种方法， 可以让你查看过去某一个时间点的一部分网页， 就是网络时光机 https://web.archive.org/。通过搜索网址， 就可以查看该网址指向的网站， 在某年某月某日的网页。

网络时光机并不是真正网站回到了过去， 而是在某一个时间点， 通过下载软件， 下载了某一个指定网站的完整网页， 备份存档。 这样在将来就可以访问，再回到这个界面，仿佛回到了过去一样。

1998年， 我第一次接触电脑， 那时候PC的操作系统还是DOS。有一次上计算机课， 启动电脑时我不小心选择了Windows95，老师看到后赶紧过来吼到： “谁让你开Windows了？ Windows界面化操作， 那么复杂， 属于“超纲”的内容，不能打开”。 所谓“超纲”， 就是超出教学大纲的范围， 一般是比较难的教学内容， 学生可能很难理解。 当时我们学习DOS和QuickBasic， 对世界上其他地方发生的事情一无所知。

其实早在1986年， 美国就建成了连接主要大学的NSFNET网络， 网络通讯协议在更早的七十年代就已经确立， 也就是说，互联网在那时候就已经诞生了。1998年， 互联网已经在蓬勃发展。重要的机构陆续建立了自己的网站。

上大学后当然也少不了学习计算机， 这就要提到学校的机房。要进入机房， 先要买一副鞋套，要保证机房里面一尘不染。 每人带着一张1.44M的软盘，机箱上有一个软驱， 重要的资料，一般是计算机作业， 或者要打印的文档， 就都保存在这张软盘内，以便在机房的电脑上修改。

不过机房里的电脑实在是不争气。因为是投射显示器，眼睛看一会儿就回觉得发酸流泪， Windows98系统也经常无缘无故的死机， 再加上机房用电负荷增加，经常跳闸断电更为修改电脑作业增加了无数烦恼： 要么文件损坏， 要么几个小时的打字内容因为断电没能保存。 宏病毒是Office 软件的座上客， Word、PowerPoint和Excel可能在毫无征兆下中毒， 再打开后内容被清空， 也无法恢复。由于计算机配置太低，要读取PDF，打开Adobe Reader 都要几分钟的时间。2005年前后 ，1.44寸软盘终于逐渐被U盘取代， 谁那时候有一个128M的U盘， 也能足够风光一阵了。

QQ号成了非常紧俏的商品， 帮忙申请QQ号是讨女孩子欢心的重要法宝。 坐在网吧里，跟相邻的电脑网聊居然会成为一种享受。 当然， 更多人是在打游戏， 暗黑、魔兽、CS等等， 不过由于我不打游戏， 所以也说不上来很多。网吧里， 除了打游戏的，网络对骂的，看电影的，呼呼大睡的，二手烟更是充斥了网吧的每个房间， 电脑桌甚至键盘上，都不难找到被烟灰烫伤的地方。 有些大学生甚至连续几个月不上课， 也不回宿舍， 饿了叫外卖，完全在网吧里过日子。 我的一位室友就因为抑郁， 在网吧打游戏， 连续五个月没有回宿舍。 把他找回来的时候， 他脸色惨白， 人几乎瘦成了骷髅， 我都担心他能不能站起来， 担心一旦站起来，风会不会把他吹倒。虽然沉迷在网吧的学生有不少， 学校门口的“超越”和“飞跃”两个网吧还是二十四小时开放，后来据说加强监管，不过因为我去的很少，就不清楚情况了。

那时候我特别迷恋一个植物学网站， 网站叫做 “原本山川，极命草木”， 这是中科院植物所的植物爱好者开设的网站，因为网址为 http://www.emay.com.cn, 大家都称它为 “义妹” 论坛。也就是在那个论坛， 我认识了山东大学的刘冰、北京大学的刘夙、北京林业大学的林秦文、植物所高贤明等几位老师，那时候都是网友。 我的注册网名是 mrstarlitsky，头像是一个星系。

管理员“静生”在论坛里非常活跃， 帖子数长期居第一位， 发言也毫无忌惮， 我有时候怀疑他帖子的真实性。静生对在同一个单位工作的同事多有指责， 虽然有些时候不无道理， 但是现在想起来有些内容是不便在网络发表的。

emay论坛里面最吸引人的， 就是有人贴上天南海北的植物求助鉴定， 然后谁能提供相应的线索， 也许就能帮上大忙， 也有人分享自己考察的照片， 展示当地特别的植物。我在大一以后痴迷于认植物，常常抱着书本不放， 抄写了《东北植物检索表》的全部拉丁名和中文名， 抄写了《中国种子植物科属检索表》的各属， 并且努力背诵各科属之间的区别。图书馆的植物图鉴和地方植物志，被我借了个遍。 emay论坛也是我一展身手的好地方，除了北方的植物， emay上南方的照片很多， 这让我感觉南方的植物更丰富， 所以一直非常想到南方看看。也许这就是为什么我现在会在南方工作。

emay网站好像有一种魔力， 让人不断分享， 也不断收获。 后来的普兰塔论坛也有这种魔力。 由于种种因素， 这两个论坛都没有能够坚持下来， 让人觉得很可惜。这两个网站是承载了自己成长过程中的重要回忆。

日月如梭， 过去的时光是回不去了。 要了解网站在当时的样子，就在web.archive.org看看吧。

修改综述遇到的一些问题

修改论文、综述或者申请书，真的很累， 甚至感觉比自己写还要累。总结起来， 大概有以下方面的问题：

1. 缺乏中心思想，看不出作者的观点。 看了半天，不知道作者想说什么。

2. 缺乏清晰的结构。结构不清晰，想要看懂不容易。

3. 行文冗长，絮絮叨叨。唉，作者高兴就好……

4. 内容空泛。XXX将有力促进生物多样性保护，改善环境，保护绿水青山。是的，我也同意，但是是怎么保护的？ 说得太大， 我有点儿不信呢。

5. 引入不必要的概念。本来能用一句话说清楚的事， 非要弄出一个新概念或新的专业词汇。抽象出新词汇或者引用过时的概念都不太好。如，XXX将促进“生态平衡”……问题是现在谁还用生态平衡？

6. 观点前后矛盾。刚刚说的XXX很重要， 下一段马上说， XXX不重要了。真是180°大转弯，变化太快，小心刹不住车啊。

7. 逻辑不清。本来大部分文字读起来就挺枯燥的， 你再来个不按套路出牌，后果可想而知。 如果按照，（1）是什么问题，（2）问题的后果，（3）问题的成因，（4）怎么解决问题的顺序写，读者就容易理解。 没错， 思维是发散的， 思想是立体的， 但文字是线性的。既然如此， 行文时最好能把要说的事按照一定顺序串起来， 而且必须以一种容易理解的顺序串起来。如果一个人先说（1）问题该怎么解决，再说（2）问题的后果，（3）什么问题， 最后再（4）解释问题的成因， 读者心里恐怕要骂人了。

8. 缺乏论据。论述的要求就是有理有据，还是摆点儿事实比较好，因为这些是你的论据啊 。虽然论点可能是你最早意识到的，但是还是要收集证据去支持论点， 做到有理有据， 否则又流于空泛了。

9. 论据选择不当。举例：要说明物种分布区预测的算法存在缺陷， 但后面引用的是关于物种分布数据不足的论文， 然后一直在强调数据是怎样不足……这，我还真不知道说点儿什么好。

10.用词不准。专业人士，用词还是准确一点儿好。比如phylogenetic diversity, 就翻译为“系统发育多样性”吧， 至于“谱系多样性”，额，麻烦再改改吧，倒不是说完全不行， 而是“谱系”已经另有所指，不要再和“系统发育多样性”掺和了。这让我想起几种植物：苦苣菜、苦荬菜、苣荬菜、苦菜、小苦荬……命可真够苦的。做研究本来就不容易，不要再用中文增加困难了……

11.语法错误。 英文没有学好的童鞋们， 请继续努力哦

12.英式中文和中式英文。“DNA条形码图书馆的组建与物种丰富的自然群落的关系特别密切，物种准确鉴定将促进详细的生态学司法研究。”……好了， 你爱说什么说什么， 这中文我看起来不太爽……别拦我，我下次直接百度翻译。在这里，我要友善提示一下，要“小心落水”哦：Take care to fall into water….

13.图表质量太差。重要的事情说三遍： 三线表， 三线表， 三线表；分辨率， 分辨率， 分辨率……

14.引用文献不足或缺少重要细节。做进化树，还是要介绍用什么软件建的树吧；同理， 物种分布区预测， 也要说清楚是用什么软件预测的。All the statistical analysis were carried out in R， 就这么一句话就概括了数据分析部分，真的很厉害……但是，你说的我做不出来， 我面壁思过去了

用RSS追踪研究植物生态学研究进展简明教程

  1. 下载附件 plant_ecology rss.zip 解压缩，有两个文件 RSSOwl.2.2.1.Setup.exe 和 plant_ecology.opml 文件

  2. 双击安装RSSOwl.2.2.1.Setup.exe开始安装RSS浏览器， 选默认就好， 安装过程没什么好说的。因为RSSOwl是用Java写的， 因此如果没有安装Java Runtime的电脑， 需要先安装Java SE （http://www.oracle.com/technetwork/java/javase/downloads/index.html）。

  3. 批量导入RSS数据源

点击 File> Import > 在Choose Source 界面选择 Import Feeds from a file or website > Browse , 选择 RSS for plant ecology 20180403.opml, 点击open > next > Finish

如果之前已经导入过opml文件， 则不能重复导入。

  1. 单独导入RSS数据源， 以Journal of Ecology为例

首先要找到该杂志的网页， 找到RSS图标， 点击图标， 可以见到RSS网址，或者RSS的内容， 这个网址称为Feed。按住Contrl + C 复制该网址。

点击File下方的 + ，在新弹出的New Feed 窗口选择 Create a new feed to read news,复制的网址应该已经自动填写好了，点击Finish即可。

plant_ecology rss.zip

RSS追踪科研进展.pdf