通过重命名的照片整理植物名录

在生物多样性调查中， 数码照片是快速记录物种信息的手段之一。 调查人员如果对一个地点的植物已经足够熟悉， 无需将所有植物踏查结果写到记录本上， 而是对见到的每种植物拍照， 可回到办公室再进行鉴定和整理。 这样既提高了野外工作的效率， 也提高了鉴定的准确性。 植物照片， 作为的一广义的标本，若有疑问， 可以供专家进一步研究和分析， 并且， 实在有必要， 可以回到野外， 采集标本。

拍摄的照片， 一般是按照地点， 放到文件夹中。按照相机的品牌和型号不同， 照片的名称以IMG_，DSC或DSCN开头，后面是照片的顺序号， 除专业摄影师， 方便后期的处理， 保存为别的格式之外， 一般情况下， 照片均保存为.JPG格式。.JPG的扩展名一般不要修改， 否则在Windows中打开文件会遇到困难。

对照片的重命名， 一般用标准中文名称（中国植物志的正式中文名）， 这样便于管理。依据照片中的物种， 为植物重命名， 一般加在.JPG之前， 而前面的IMG以及顺序号均保留。 一些浏览照片的程序包，如ACDSee， XnView等， 提供重命名的快捷方式， 以及批量重命名的方法；天南星等网友准备的植物名称词库， 也可以快速准确的输入中文植物名，因此， 如果对照片中的植物熟悉的话， 对照片重命名一般可以很快完成。

有了重命名的照片， 可以用以下方法，快速生成植物名录。这里仍然用本人比较熟悉的R软件：

（注意： 本功能已经置于 R 程序包 plantlist中, http://blog.sciencenet.cn/home.php?mod=space&uid=255662&do=blog&id=813796 ， 用户只需给出路径， 即可获得结果， 该结果同时包括APGIII系统的分科）

以下代码仅供参考， 请使用最新版（Ver 0.0.4）的 tpllist
定义一个R函数：

1
2
3
4
5
6
7
8
9
10
11
12
cnplist <- function(path){
   setwd(path)
   temp.photo.list <- list.files()
   photo.list <- temp.photo.list[grepl(".JPG", temp.photo.list)]
   list.res <- list()
   sp.list.cn <- unique(gsub(".JPG","",gsub("[0-9]","",gsub("DSC", "", gsub(pattern = "DSCN", "", gsub(pattern = "IMG_", replacement = "", x = photo.list))))))
   sp.list.cn <- sp.list.cn[!sp.list.cn == ""]
   sp.list.cn.data.frame <- as.data.frame(sp.list.cn)
   res <- merge(x = sp.list.cn.data.frame,  y = spdat, by.x = "sp.list.cn", by.y = "NAME_CN", all.x = TRUE )
   write.csv(res, paste("checklist.csv", sep = ""))
   return(res)
}

读取中文名学名对照表（参见本文附件 plantnames.csv ）：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
spdat <- read.csv("plantnames.csv", header = TRUE)

#### 函数的使用方法举例
#### 列出梧桐寨瀑布文件夹下， 所有物种的名录
cnplist("C:\Plant photos database\2014\20140120梧桐寨瀑布")


#### 列出2014目录下，每个文件夹中包括的物种名：

dirs <- "C:\Plant photos database\2014"
aaa <- list.dirs(dirs)
for(i in 1:length(aaa)){
   cnplist(aaa[i])
}

2014年1月版纳园 群落系统发育数据分析研讨班幻灯片

Community Phylogenetics Workshop: Data analysis using Phylomatic, Phylocom and R

第一天 从DNA序列到进化树 xtbg day 1.pdf

第二天 群落数据处理和进化树的整合分析 xtbg day 2.pdf

本套幻灯片，是2014年1月， 在西双版纳热带植物园的讲座时准备的。主要内容如下：

探讨群落物种组成的成因， 一直是群落生态学的核心问题之一。 近些年来， 随着将进化信息和群落数据整合，

人们得以从一些新的角度探讨群落物种组成的成因。 但在群落系统发育的数据分析中， 从DNA序列如何建立进化树？ 群落数据应如何处理？ 在数据分析的过程应该注意哪些问题？如何解释群落系统发育获得的结果，用以验证群落物种组成的中性理论或生态位理论？ 群落系统发育分析有哪些不足以及未来的发展方向如何？

这次讲座较为系统得回顾了上述知识， 结合相应的练习， 希望对本领域的研究人员入门和进行数据分析有帮助。

2014年1月 于西双版纳

HK80：香港所用地理坐标之间的相互转换

目前， 香港所用的地理坐标，主要包括：

  1. WGS84 经纬度
  2. 基于WGS84椭球的UTM坐标，其中东经114度以东为50区， 以西为49区
  3. HK1980坐标， 类似于UTM， 但是该坐标系的原点位于香港大屿山岛的西南， 以保证香港所有的坐标均为正。
  4. 基于HK1980 所用椭球的经纬度
  5. 基于HK1980所用椭球的UTM坐标

香港地政署提供了坐标系的介绍， 以及相应的转换公式， 参数 （见 http://www.geodetic.gov.hk/smo/gsi/data/pdf/explanatorynotes.pdf），同时提供了相应的网上工具（http://www.geodetic.gov.hk/smo/tform/tform.aspx）， 以完成坐标系之间的相互转化。按照其中的公式， 转换后的精度， 在5m以下， 已经足够进行物种分布调查等用途。

但是该工具很难直接进行批处理， 需要逐条粘贴， 运行时较易出错。 为此， 笔者根据说明中的公式和相应的算法， 编写了相应的函数， 并制作了R程序包， 并定名为HK80。 HK80可以完成五种坐标系之间的相互转换， 共二十个函数， 详情请参阅帮助文件。现将该程序包共享， 以方便有需要的人士。 并请用户提出宝贵意见， 以便将可能的错误及时更正。

张金龙 谨识

2014年3月20日

HK80_0.0.1.tar.gz Linux 源代码

HK80_0.0.1.zip Windows版本

HK80-manual.pdf 程序包使用说明

基于数码相机和GPS整理物种坐标的一个方案

需要的硬件: 数码相机, GPS

需要的软件: GeoSetter, XnView, 7zip

  1. 拍摄时： 在拍摄照片前， 相机的时间要和GPS调整成一致， 最好精确到5S以内

  2. 用XnView为照片重命名， 具体为按下F2， 重命名即可，至少提供中文正式名称

  3. 为照片添加GPS坐标。 GPS导出的航迹保存为GPX格式， 这样可以记录每个航迹点的时间。 保存到照片所在的文件夹下，注意路径不能有中文。用GeoSetter匹配每张照片的拍摄时间， 以及GPS航迹上各坐标点的时间， 为每张照片添加GPS坐标。GeoSetter会将信息写入照片的EXIF信息中。

  4. 用GeoSetter将航迹和照片导出为kmz文件，用7zip解压缩kmz文件， 将其中的kml文件用KMLCSV转换为csv即可。

使用knitr处理Rnw文档

Rnw文档是保存了R代码的Latex文档。人们希望Latex中的R代码可以自动运行， 并将运行结果， 复制到Latex代码中。 同时， 给出绘制的图形等。

例如以下的代码， 是一篇Rnw文档的内容， 将其拷贝， 粘贴到记事本中， 另存为test.Rnw文件即可。 注意编码要选ANSI。

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
\documentclass[CJK]{cctart}
\usepackage{verbatim}
\title{中文文章里面如何使用knitr？ }
\author{张金龙}
\begin{document}
\maketitle
在本例子中， 我们将演示在R中， 如何调用 hclust() 并且将结果写入 LaTeX{} 文件中:
<<>>=
x<-c(1,2,6,8,11)
dim(x)<-c(5,1)
d<-dist(x)
d
hc<-hclust(d, "single")
dend<-as.dendrogram(hc)
@
绘制的树状图:， 目前，kintr程序包在处理Rnw文档时， 不允许注释中有中文
begin{center}
<<fig=TRUE,echo=FALSE>>=
par(mfrow = c(2, 2),mar = c(4,3,1,2))
plot(dend)
plot(dend, nodePar=list(pch = c(1,NA),
   cex=0.8, lab.cex=0.8), type = "t",
   center=TRUE)
plot(dend, edgePar=list(col = 1:2,
   lty = 2:3), dLeaf=1,
   edge.root = TRUE)
plot(dend, nodePar=list(pch = 2:1,
   cex=.4*2:1, col=2:3),
   horiz=TRUE)
@
\end{center}
\end{document}

在Rnw文档中， R代码需要放在 <<fig=TRUE,echo=TRUE>>= 和 @ 之间。 其中的代码一定要可以执行。 在后续处理中，R的Sweave函数和knitr程序包中的knitr函数， 都能够生成相应的tex文件， 以及所需要包含的pdf图形。

fig=TRUE 表示， 要包括代码生成的图形；
echo=TRUE 表示， 要包括运行结果所生成的代码。

Sweave函数就是为了处理Rnw文件而生的。 但是Sweave在处理中文的Latex文档时， 需要设定的参数很多， 并且常常出错。幸运的是，我们还有谢益辉编写的knitr程序包。

kintr在处理Rnw文档非常方便。 假设Rnw文档放在

C:/Users/Jinlong/Desktop/Sweave

文件夹之下，名为test.Rmw， 则使用以下R代码:

1
2
3
setwd("C:/Users/USER/Desktop/Sweave") ### 转换到工作目录，
library(kintr) ###导入kintr程序包
knit("test.Rnw")

即可在该目录下， 生成.tex文档， 以及需要包含的文件。

用TexWorks编译成pdf文档即可。

注意：kintr程序包在处理Rnw文档时，R代码的注释中不允许注释中有中文

苍茫的大地－香港麦理浩径第八段

麦理浩径是香港最早开辟的行山路径，全长达100公里， 于1979年10月启用。 沿途经过西贡、马鞍山、 狮子山， 大帽山， 城门等地， 沿途风景优美。

麦理浩径第八段，从海拔从400m的铅矿坳开始，最高达950米。

这是冬天里的麦理浩径： 在碧蓝天空的笼罩下之下， 弯弯曲曲的小路穿行在枯黄的野草中， 点缀在路旁的巨石，不知是谁用心布置的， 感觉真是恰到好处。 下午的阳光， 暖暖得照着， 凉风习习， 举目望见那点点色彩， 就是行山的人们，更让人觉得人的渺小和大自然的辽阔与悠远。

路就在脚下， 前行， 才是希望所在。

信心， 给我们希望； 前行， 给我们风景。

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

用skycalc计算日出日落时刻

给定经纬度， 要计算某年某月某日该地点的日出日落时刻， 考虑的因素包括：

  1. 太阳在天球上的位置 （赤经,赤纬）： 一般要获得计算当日， 以及前一日， 后一日的太阳赤经赤纬， 通过插值的方法， 获得太阳在日出日落时刻的准确位置

  2. 城市所在的时区 （天文学上的习惯是， 在格林尼治子午圈以东为负， 以西为正） ，以便将计算获得的世界时转换为当地时区的时间。

这里介绍用skycalc程序包如何计算一个地点， 给定日期的日出日落时刻。

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
library(skycalc)
### 儒略日转换为年月日时分秒的函数
Julian2Date <-function(jd){
int2 <-
function(v){
  return(floor(v));
}
   r =list();
  D = int2(jd+0.5);
  F= jd+0.5-D #取得日数的整数部份A及小数部分F
  if(D >=2299161){
       c= int2((D-1867216.25)/36524.25);
       D= D +1+c- int2(c/4);
   }
  D  = D +1524;              
  r$Y = int2((D -122.1)/365.25);#年数
  D  = D - int2(365.25*r$Y);  
  r$M = int2(D/30.601); #月数
  D  = D - int2(30.601*r$M);
  r$D = D; #日数
  if(r$M>13){
      r$M = r$M -13;
      r$Y = r$Y -4715;
   }else{
      r$M = r$M -  1;
      r$Y = r$Y -4716;
   }
  #日的小数转为时分秒
  F = F *24;
  r$h=int2(F);
  F = F - r$h;
  F = F *60;
  r$m=int2(F);
  F = F -r$m;
  F = F *60;
  r$s=F;
  return(r);
}

sunRTScalc <-function(year, month, day, longitude, latitude,  zone,
                      bGregorianCalendar =TRUE,
                      type =c("rise", "transit", "set")){
   type <- match.arg(type)
   JD = CAADate_DateToJD(year, month, day, bGregorianCalendar);
   SunDetails1 = CAAElliptical_Calculate_Sun(JD -1)
   Alpha1 = SunDetails1$ApparentGeocentricRA;
   Delta1 = SunDetails1$ApparentGeocentricDeclination;
   SunDetails2 = CAAElliptical_Calculate_Sun(JD)
   Alpha2 = SunDetails2$ApparentGeocentricRA;
   Delta2 = SunDetails2$ApparentGeocentricDeclination;
   SunDetails3 = CAAElliptical_Calculate_Sun(JD +1)
   Alpha3 = SunDetails3$ApparentGeocentricRA;
   Delta3 = SunDetails3$ApparentGeocentricDeclination;
   RiseTransitSetTime =CAARiseTransitSet_Calculate(JD,
   Alpha1, Delta1, Alpha2, Delta2, Alpha3,Delta3,
   longitude, latitude, -0.8333);
   if(type =="rise"){
       rtsJD =(JD +(RiseTransitSetTime$details.Rise/24.00));
   }
   if(type =="set"){
       rtsJD =(JD +(RiseTransitSetTime$details.Set/24.00));
   }
   if(type =="transit"){
       rtsJD =(JD +(RiseTransitSetTime$details.Transit/24.00));
   }
   lclJD = rtsJD -(zone/24.00);
   return(Julian2Date(lclJD))
}
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
### 例一
### 计算美国波士顿 (西经 74.73057,  北纬 39.275787) 2009年8月9日的日出日落时刻
sunRTScalc(year =2009, month =8, day =9, longitude =74.73057,  
          latitude =39.275787, zone =4, bGregorianCalendar=TRUE, type ="rise")
### 因为波士顿的日落已经发生在格林尼治子午圈的下一天 (波士顿当地时间(西四区)2009年8月9日20h01m39m，对应世界时2009年8月10日00h01m39m)， 所里这里在日期上需要加1.

sunRTScalc(year =2009, month =8, day =10, longitude =74.73057,  
          latitude =39.275787, zone =4, bGregorianCalendar=TRUE, type ="set")


### 例二
### 计算北京天安门 (东经 -116.391325,  北纬 39.907356) 2014年1月5日的日出日落时刻
### 日出 因为北京的日出时刻发生在北京时间4h45m43.7s，对应的世界时为1月4日23h36m13.8s， 因此， 在日期上需要减去一天
sunRTScalc(year =2014, month =1, day =4, longitude =-116.391325,  
          latitude =39.907356, zone =-8, bGregorianCalendar =TRUE, type ="rise")
### 日落
sunRTScalc(year =2014, month =1, day =5, longitude =-116.391325,  
          latitude =39.907356, zone =-8, bGregorianCalendar =TRUE, type ="set")
### 中天
sunRTScalc(year =2014, month =1, day =5, longitude =-116.391325,  
          latitude =39.907356, zone =-8, bGregorianCalendar =TRUE, type ="transit")

### 例三
### 计算2014年6月20日北京天安门 (东经 -116.391325,  北纬 39.907356) 的日出日落时刻
### 日出 因为北京的日出时刻发生在北京时间4h45m43.7s，对应的世界时为6月19日22h45m43.7s， 因此， 在日期上需要减去一天
sunRTScalc(year =2014, month =6, day =19, longitude =-116.391325,  
          latitude =39.907356, zone =-8, bGregorianCalendar =TRUE, type ="rise")
### 日落
sunRTScalc(year =2014, month =6, day =20, longitude =-116.391325,  
          latitude =39.907356, zone =-8, bGregorianCalendar =TRUE, type ="set")
### 中天
sunRTScalc(year =2014, month =6, day =20, longitude =-116.391325,  
          latitude =39.907356, zone =-8, bGregorianCalendar =TRUE, type ="transit")

####例四
### 计算2020年1月15日澳大利亚墨尔本(东经 -144.963171, 南纬-37.814247) 的日出日落时刻
### 墨尔本为位于东10区，澳大利亚夏时制开始于上一年10月的第一个星期天， 下一年4月第一个星期天止。1月15日，正处于夏令时， 因此时区取-11
### 日出, 由于澳大利亚日出时， 世界时仍然处在2020年1月14日， 因此这里day取14
sunRTScalc(year =2020, month =1, day =14, longitude =-144.963171,  
          latitude =-37.814247, zone =-11, bGregorianCalendar =TRUE, type ="rise")
### 日落
sunRTScalc(year =2020, month =1, day =15, longitude =-144.963171,  
          latitude =-37.814247, zone =-11, bGregorianCalendar =TRUE, type ="set")
### 中天
sunRTScalc(year =2020, month =1, day =15, longitude =-144.963171,  
          latitude =-37.814247, zone =-11, bGregorianCalendar =TRUE, type ="transit")

AAplus天文算法C++库在R中的实现：skycalc程序包

1991年， 比利时天文学家Jean Meeus的著作 天文算法 出版。 该书分为四十几个章节， 每章介绍一些天象的计算方法， 同时介绍用计算机如何实现。 该书出版后很快吸引了众多的天文爱好者。 当时， 真正实现计算机窗口化的Windows95还没有诞生， 但是从那时起， 利用Jean Meeus提出的算法， 用计算机程序来进行天象计算，就已经成为很多对天文和计算爱好的一部分。

该书从插值方法、迭代方法开始讲起， 继而介绍日期的转换， 节日的确定以及日月和各大行星的位置， 亮度的计算方法和各种坐标系统的相互转换，如地平坐标、赤道坐标、黄道坐标、 银河系坐标、 太阳中心坐标等，至此， 可以方便计算各天体的出没时刻以及晨昏蒙影的时刻。 书中还包括土星环， 土星卫星， 木星卫星， 各行星近日点时刻， 春分， 夏至， 秋分， 冬至的计算等等。该书采用VOS87行星运动理论，星体位置的计算， 考虑到了多方面的摄动影响， 用多项式表示天体位置， 达到了极高的计算精度。因此， 该书从出版至今，被认为是学习天文历算的最佳入门和练习教程之一。

国内外的很多天文爱好者， 一直在致力于用计算机语言实现该书所介绍的各种算法，并因此诞生了多个语言的多个版本。如

其中， 爱尔兰软件工程师 P.J. Naughter 用C++写成的AAplus库， 按照章节， 实现了天文算法的大部分算法， 并且公布了源代码（http://www.naughter.com/aa.html）。 该库大约每一两个月更新一次， 更正其中的错误， 添加一些新的参数等。每个函数， 都有对于参数和返回值的详细说明， 并有检验的程序， 对函数的功能进行检验。所以， 对于对天文历算感兴趣的人来说，该C++库可以说是极好的入门材料。

由于是C++的源代码库， 该库的可移植性应该是很好的， 可以在多个平台上运行。但是C++必须要经过编译和连接之后才能执行， 因此， 每个函数在调用的时候， 都需要用C++写相应的main函数，这为程序调试， 学习带来了很多不便之处。

由于R语言可以直接给出函数的运行结果，在调用时， 非常方便。 函数的运行结果， 还可以绘图， 并进一步编写星历表等。为此，从2013年11月开始， 笔者根据R语言Rcpp程序包提供的R语言和C++接口，将绝大部分AAplus的函数参数，转换为SEXP对象，SEXP对象在AAplus中可以编译执行。函数的返回值， 若为单个数值， 则用wrap函数， 转变为R可以识别的SEXP对象， 若为多个数值， 则返回R的list。 至此， AAplus库的300多个函数，其中363个函数添加了R语言的参数接口。C++的程序写完后，我便着手写R的wrapper函数，在Windows下， C++函数， 都需要先编译成dll动态链接库， 而R中， 用.C函数， 可以调用dll函数。 所有的wrapperfunction写完后， 又为每个R函数准备相应的帮助文档， 以及运行实例。 这样前后大约花了一个半月。

最后程序包取名为skycalc， 取天象计算之意， 目前的版本是0.1。

由于本人水平有限， 对于天文历算的很多问题理解不深，程序包中难免存在着很多错误， 敬请用户批评指正， 以便在后续版本中， 及时改正。

张金龙 谨识

2014年1月

于香港新界大埔

源代码及安装 https://github.com/helixcn/skycalc

三维图中的颜色渲染

在R中绘制三维图常用persp， 但是persp函数并未提供颜色渲染的相应参数。这里提供了一种方法，在persp中， 按照取值不同， 添加相应的颜色。

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
### 给出x，y以及f取值序列， 要绘制的是f的函数面
x <-seq(-10, 10, length=30)
y <- x
f <-function(x, y){ r <-sqrt(x^2+y^2); 10*sin(r)/r }
z <-outer(x, y, f)
z[is.na(z)]<-1
op <-par(bg ="white", mar=c(1,1,1,1))
### 生成颜色序列
col.table1 <- topo.colors(400)
col.table2 <- heat.colors(400)
col.table3 <- terrain.colors(400)
col.table4 <- cm.colors(400)
### 计算行列
ncz <-ncol(z)
nrz <-nrow(z)
### zfacet是用来做cut的， 如果想获得与矩阵相同的grid数目， 需要在行列的开始和结束，都减去1
### 这里的zfacet是 breaks
zfacet <- z[-1, -1]+z[-1, -ncz]+z[-nrz, -1]+ z[-nrz, -ncz]
### 进行cut， 生成400种颜色, 生成每个方格，所属的颜色等级
col.ind <- as.numeric(cut(zfacet, 400))
###col.table1[col.ind] 本句是将col.table1颜色序列， 按照zfacet的值进行排列， 以保证对应的方格， 出现对应等级的颜色。 因为z矩阵的保存方式实际上是一个vector。 col.ind能够才能够和grid的颜色匹配。
par(mfrow =c(2,2), mar ＝ c(1, 1, 1, 1))
### 四种颜色的三维图
persp(x, y, z, theta =30, phi =30, expand =0.5, col= col.table1[col.ind], main ="topo.colors")
persp(x, y, z, theta =30, phi =30, expand =0.5, col= col.table2[col.ind], main ="heat.colors")
persp(x, y, z, theta =30, phi =30, expand =0.5, col= col.table3[col.ind], main ="terrain.colors")
persp(x, y, z, theta =30, phi =30, expand =0.5, col= col.table4[col.ind], main ="cm.colors")
注： lattice程序包的wireframe函数提供了颜色的渲染方法。
wireframe(volcano, drape =TRUE, aspect =c(61/87, 0.4))
wireframe(volcano, shade =TRUE, aspect =c(61/87, 0.4))
wireframe(z, shade =TRUE, aspect =c(61/87, 0.4))
wireframe(z, drape =TRUE, aspect =c(61/87, 0.4))

新西兰Dunedin的杜鹃花海

新西兰南部的小城但尼丁， 10月下旬正好迎来春天。 春光明媚。这里是杜鹃花的海洋， 四处绽放着杜鹃花。 但尼丁植物园更有一个杜鹃花展区。徜徉在杜鹃花海中， 春天的气息沁人心脾。

这里的杜鹃花大多是我国引种的。 在惊叹当地植物园杜鹃花品种的保存以及育种的专业的同时，我不由感叹，我国虽是这些杜鹃花的原产国， 但是并没有将其利用好， 而是任其自然情况下自生自灭。 有些甚至遭到严重的破坏。

自然保育的观念， 还需要进一步深入人心啊

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img

img