基于R语言的风廓线雷达产品数据简单处理

于红明　刘宝君　宋志龙

摘 要：风廓线雷达产品数据在气象上应用广泛，对其产品数据资料的分析、处理和可视化十分重要。本文利用免费的开源软件R语言，对雷达资料产品数据进行批量读取、计算和绘图。展示了R语言在处理风廓线雷达资料的特点：简洁、易学。

关键词：R语言 风廓线雷达产品数据 风廓线图

中图分类号：P49 文献标识码：A 文章编号：1672-3791（2014）08（c）-0039-03

风廓线雷达作为中小尺度天气系统有效的探测工具，能够24小时不间断的提供：水平风向、风速、垂直气流等气象要素随高度的分布，是进行边界层和高空气象探测的重要设备，能对现有的气象观测进行补充。其观测资料在气象上应用广泛：对高影响天气的监测和预警；对灾害天气系统的天气学分析；用于数值预报中的观测资料同化和数据后处理等。因此开展风廓线雷达的处理十分必要[1～4]。

现有的气象数据处理一般使用Fortran，C等高级编程语言。虽然有着处理速度快，计算结果可靠等优点。但是其使用并不便：需要掌握编程语法和算法，不同计算机平台可能需要不同的编译系统等。所以本文针对风廓线雷达资料，使用R语言来对其进行简单处理。R语言的全称：为统计计算和图形展示而设计的一种编程语言和统计环境。目前有一个R核心开发团队对其进行定期维护和更新[5]。其主要的特点在于：（1）完全免费：可以在其官方网站（http：//cran.r-project.org）下载到完整的安装包并免费使用；（2）开源软件：R语言源代码完全公开，任何人均能提供各种数据处理模块，下载相应模块后能迅速地完成数据处理等工作；（3）图形显示：可以利用R语言处理数据后，直接获取各种统计分析图形；（4）多平台使用：R语言可以在Windows、mac、Unix这些操作系统上安装，不需要重复编译，其脚本在不同操作系统间可以任意使用。针对风廓线雷达资料，我们可以利用R语言对其进行简要的数据处理，统计分析，然后利用图形显示系统得到各种分析图片。

本文将在第二部分介绍R语言读取和储存风廓线雷达产品数据资料；文章第三部分对资料进行简要处理，最后会利用R语言的图形显示功能作图。

1 雷达产品数据的读取

测站观测的风廓线雷达资料有两种数据[6]：原始数据和产品数据。其中原始数据为二进制格式，主要包含功率谱数据文件，瞬时径向谱数据文件。产品数据文件为文本格式，包括实时的采样高度上的、半小时平均的采样高度上的、一小时平均的采样高度上的产品数据文件。本文所处理资料为实时的采样高度上的产品数据文件。R语言也能读取处理二进制格式文件，对雷达原始数据资料的处理会在以后的工作中进行。

1.1 单个产品数据的读取

R语言有许多函数能够直接从文本文件中读取数据，比较常用的有：read.table（），read.csv（），read.fwf（）。其中和Fortran比较接近的是read.fwf（），可以指定读取数据的长度和格式。于Fortran不同的是，R语言在读取数据的时候，不用先给定数据类型，程序会直接读取数据，并存储到一个数据框（data.frame）里。例如针对本文要处理的雷达产品数据，可以直接使用read.table（）命令读取产品数据：

raw_data<-read.table（fname，fill=TRUE）

其中raw_data为一个数据框（data.frame），用于储存我们需要处理的雷达产品数据，fname 为要读取的雷达数据文件名字；“header=FALSE”表示该实体数据中没有数据说明头文件；由于本文要处理的雷达产品数据不是规则的表格形式，所以需要使用参数：“fill=TRUE”，来自动填满不是数据表格的部分。产品数据文件的前三行为测站基本参数，最后一行为结束行，中间部分为实际数据，包含：采样高度，水平风向，水平风速，垂直风速，水平方向可信度，垂直方向可信度，Cn2。为了便于资料处理，将实体数据单独储存到一个名为r_data的新数据框中：

r_data<-raw_data[4：（length（raw_data[，1]）-1），]

其中使用“length（raw_data[，1]）”函数判断raw_data一共有多少行，通过“4：（length（raw_data[，1]）-1）”截取其第四行到倒数第二行到新的数据框r_data中。

值得注意的是，产品数据文件中一般含有缺测值，本文件中用“////////”表示。R在读取其数据的时候将所有的数据的数据类型默认为因子（factor，一种R的数据类型）。而在数据处理计算中，我们使用的数据类型为数值形（numeric），通过如下简单语句即可实现数据类型转换：

r_data<-lapply（r_data，as.character）

r_data<-lapply（r_data，as.numeric）

r_data<-as.data.frame（r_data）

先使用函数“lapply”和“as.character”把每列元素转换为字符型，再使用as.numeric转换为数值型，最后再用as.data.frame把r_data转换成一个数据框。从这些处理过程中不难看出，R语言能十分方便的利用“lapply”函数实现整列（整行）的数据处理，相应的apply类函数还有许多，从而省去了程序中的循环语句编写；而且R的数据类型使用十分灵活，能够方便的将其转换为不同的数据类型，类似的函数还有as.matrix，as.logical等。R在使用as.numeric时，会自动将无法转换成数值的字符，转换成R的缺测值（NA），本数据中的“//////”在经过转换后全部转换为了“NA”。“NA”在R语言中可以参与计算，也可以使用一个简单的函数将数据中的“NA”去除。例如针对本文中的数据：endprint

m_data<-r_data[complete.cases（r_data），]

其中complete.case函数能够获取数据中不含缺测的所有列，进而赋值后的数据框m_data中剔除了r_data的缺测值。

1.2 批量数据文件的读取

通过1.1中几个简单的语句即可完成单个雷达产品数据读取和初始化，而在业务运行中，雷达产品数据肯定是大批量的生成，也需要程序脚本具有批处理功能。

我们将1.1中单个文件的雷达产品数据处理过程，整合成一个R语言函数readradar：

readradar <- function（dir，fname）{… …}

readradar只需要给定路径和文件名，即可读取雷达产品数据资料并去除缺测值。函数最后返还一个数据框，其中包含该文件中所有非缺测产品数据。给定所有需要处理的路径和文件名后，即可完成资料的批处理。

而R语言能够十分便利的获取文件名，因为其能通过脚本语言进入当前运行的操作系统。简单的说，R语言能够在程序内部完成操作系统的文件、文件夹处理、安装包的安装等。例如可以直接使用file.create（create.dir）函数直接在R语言中生成新的文件（文件夹）等。 本文只需要R语言读取要处理的文件名：

fname<-list.files（path="./Qingdao/"

fname中存储了所有“./Qingdao”下的雷达数据文件名。通过循环即可完成雷达数据的批量读取：

for （i in 1：length（fname）） {

m_data<-readradar（dir，fname[i]）

}

别的批量处理过程（资料简要处理，雷达风廓线图等）和读取类似，只需要加入循环即可。

2 雷达资料的简要处理和画图

第二节中给出了R语言对雷达产品数据文件的批读取。R语言最为实用的优点在于其计算和画图功能。

2.1 雷达风廓线的简要计算

为了便于处理资料。我们利用names函数将m_data数据框的每一列分别命名：

names（m_data）<-c（"hgt"，"h_dir"，"h_speed"，"w_speed"，"h_rel"，"w_rel"，"cn2"）

每一列表示对应的资料的采样高度，水平风向，水平风速，垂直风速，水平方向可信度，垂直方向可信度，Cn2。

水平风的u分量计算如下：

m_data$u <- m_data$h_speed * cos（（270 - m_data$h_dir）* pi/180）

其中pi为R语言中自带的圆周率的取值3.141593。而计算出来的u风速可以直接存储到数据框m_data的新的一列中。

同理可以计算出水平风速的v分量：

m_data$v <- m_data$h_speed * sin（（270 - m_data$h_dir）* pi/180）

上述计算过程可以看出，R语言中数据框的计算处理十分方便，不需要使用循环语句计算不同高度上的u/v分量，而且计算结果可以直接作为数据框新的一列存到原数据框中。而R语言除了基本的计算外，还能使用内部函数做各种复杂数学计算，例如矩阵的求逆、线性回归分析、抽样分布、显著性试验等。如果R语言自带的函数已经不能解决当前数学问题，还能上网搜索下载对应函数包。

2.2 雷达风廓线图

R语言主要有三种绘图函数：高级、低级和交互式。通过调用高级绘图函数，能在R语言中直接绘制各种统计图；低级绘图函数能够对现有的图进行修改；交互式绘图能够让用户直接利用鼠标修改图形。大量的内置函数让绘图变得十分的简易。例如本文需要的u风场随高度变化图可以通过以下函数实现：

plot（m_data$u，m_data$hgt）

在plot函数中增加各种参数，能够优化所绘图形。例如本文中使用如下命令获取风廓线图：

plot（m_data$u，m_data$hgt，type="b"，

main=title，xlim=c（-20，20），ylim=c（0，5000），

xlab="风速（m/s）"，ylab="采样高度（m）"，pch=16，col=2）

其中tpye=“b”表示曲线为点画线；main为主标题；xlim和ylim表示横纵坐标取值范围；xlab、ylab表示横、纵坐标标题；pch=16表示点为实心圆圈，col=12表示颜色为红色。

通过使用低级绘图命令如points、lines等，可以在刚画的u风场廓线后增加v风场廓线：

points（m_data$v，m_data$hgt）

lines（m_data$v，m_data$hgt）

图例则可以用lengend函数添加：

legend（"topleft"，pch=c（16，17），col=c（2，19），lty=c（1，1），legend=c（"u"，"v"），bty="n"）

R语言还能使用par函数对图片进行设置，例如本文中使用如下函数将2014年4月26日每隔3小时的风廓线显示到同一图片中：

par（mfrow=c（2，4），mar=c（5， 4， 4， 2） + 0.1，oma=c（0.1，0.1，0.4，0.1））

其中mfrow=c（2，4）表示图片分割为两行四列，参数mar和oma指定了图片的间距。

将上述命令和第二节的数据读取函数整合到一个R脚本中，运行脚本即可得到如图1风廓线图。

R语言的绘图函数很多，除了文中绘制风廓线图外还能绘制：直方图、散点图、饼图等。在安装绘图包后还能绘制3D图等。其绘图功能，能满足大部分气象资料统计分析的出图需求。

3 结语和讨论

免费的开源软件R语言，能够用于数据分析和图形显示。文中使用其对风廓线雷达资料进行了批量读取、计算和绘图。

使用read.table函数简洁的实现资料读取。

R语言的操作系统函数，能便利的实现资料的批处理。

R语言的各种绘图函数，能迅速地绘制较为美观的风廓线图。

R语言的统计计算功能十分强大，本文只是简单的使用基础计算，在随后的工作中，可以使用R语言实现，雷达数据的质量控制，统计分析等。

本文只分析了雷达产品数据，R语言也能处理二进制数据资料，以后的工作中可以利用R语言分析其他类型的气象资料。

参考文献

[1] Green，J.L.，Atmospheric measurements by VHF pulsed Doppler radar.IEEE Trans.On Geoscience electronics.GE-17（4）：262-280.1979

[2] 何平.相控阵风廓线雷达[M].气象出版社，2006：104—122.

[3] 王欣，卞林根，彭浩，李剑东.风廓线仪系统探测试验与应用[J].应用气象学报，2005，16（5）：693-698.

[4] 胡明宝.风廓线雷达数据处理和应用研究[D].南京信息工程大学博士论文，2012.

[5] John M.Chambers.Facets of R.The R Journal，1（1）：5-8，2009.

[6] 关于进行风廓线雷达数据传输的通知.气测函[2011]61号.内部资料endprint

m_data<-r_data[complete.cases（r_data），]

其中complete.case函数能够获取数据中不含缺测的所有列，进而赋值后的数据框m_data中剔除了r_data的缺测值。

1.2 批量数据文件的读取

通过1.1中几个简单的语句即可完成单个雷达产品数据读取和初始化，而在业务运行中，雷达产品数据肯定是大批量的生成，也需要程序脚本具有批处理功能。

我们将1.1中单个文件的雷达产品数据处理过程，整合成一个R语言函数readradar：

readradar <- function（dir，fname）{… …}

readradar只需要给定路径和文件名，即可读取雷达产品数据资料并去除缺测值。函数最后返还一个数据框，其中包含该文件中所有非缺测产品数据。给定所有需要处理的路径和文件名后，即可完成资料的批处理。

而R语言能够十分便利的获取文件名，因为其能通过脚本语言进入当前运行的操作系统。简单的说，R语言能够在程序内部完成操作系统的文件、文件夹处理、安装包的安装等。例如可以直接使用file.create（create.dir）函数直接在R语言中生成新的文件（文件夹）等。 本文只需要R语言读取要处理的文件名：

fname<-list.files（path="./Qingdao/"

fname中存储了所有“./Qingdao”下的雷达数据文件名。通过循环即可完成雷达数据的批量读取：

for （i in 1：length（fname）） {

m_data<-readradar（dir，fname[i]）

}

别的批量处理过程（资料简要处理，雷达风廓线图等）和读取类似，只需要加入循环即可。

2 雷达资料的简要处理和画图

第二节中给出了R语言对雷达产品数据文件的批读取。R语言最为实用的优点在于其计算和画图功能。

2.1 雷达风廓线的简要计算

为了便于处理资料。我们利用names函数将m_data数据框的每一列分别命名：

names（m_data）<-c（"hgt"，"h_dir"，"h_speed"，"w_speed"，"h_rel"，"w_rel"，"cn2"）

每一列表示对应的资料的采样高度，水平风向，水平风速，垂直风速，水平方向可信度，垂直方向可信度，Cn2。

水平风的u分量计算如下：

m_data$u <- m_data$h_speed * cos（（270 - m_data$h_dir）* pi/180）

其中pi为R语言中自带的圆周率的取值3.141593。而计算出来的u风速可以直接存储到数据框m_data的新的一列中。

同理可以计算出水平风速的v分量：

m_data$v <- m_data$h_speed * sin（（270 - m_data$h_dir）* pi/180）

上述计算过程可以看出，R语言中数据框的计算处理十分方便，不需要使用循环语句计算不同高度上的u/v分量，而且计算结果可以直接作为数据框新的一列存到原数据框中。而R语言除了基本的计算外，还能使用内部函数做各种复杂数学计算，例如矩阵的求逆、线性回归分析、抽样分布、显著性试验等。如果R语言自带的函数已经不能解决当前数学问题，还能上网搜索下载对应函数包。

2.2 雷达风廓线图

R语言主要有三种绘图函数：高级、低级和交互式。通过调用高级绘图函数，能在R语言中直接绘制各种统计图；低级绘图函数能够对现有的图进行修改；交互式绘图能够让用户直接利用鼠标修改图形。大量的内置函数让绘图变得十分的简易。例如本文需要的u风场随高度变化图可以通过以下函数实现：

plot（m_data$u，m_data$hgt）

在plot函数中增加各种参数，能够优化所绘图形。例如本文中使用如下命令获取风廓线图：

plot（m_data$u，m_data$hgt，type="b"，

main=title，xlim=c（-20，20），ylim=c（0，5000），

xlab="风速（m/s）"，ylab="采样高度（m）"，pch=16，col=2）

其中tpye=“b”表示曲线为点画线；main为主标题；xlim和ylim表示横纵坐标取值范围；xlab、ylab表示横、纵坐标标题；pch=16表示点为实心圆圈，col=12表示颜色为红色。

通过使用低级绘图命令如points、lines等，可以在刚画的u风场廓线后增加v风场廓线：

points（m_data$v，m_data$hgt）

lines（m_data$v，m_data$hgt）

图例则可以用lengend函数添加：

legend（"topleft"，pch=c（16，17），col=c（2，19），lty=c（1，1），legend=c（"u"，"v"），bty="n"）

R语言还能使用par函数对图片进行设置，例如本文中使用如下函数将2014年4月26日每隔3小时的风廓线显示到同一图片中：

par（mfrow=c（2，4），mar=c（5， 4， 4， 2） + 0.1，oma=c（0.1，0.1，0.4，0.1））

其中mfrow=c（2，4）表示图片分割为两行四列，参数mar和oma指定了图片的间距。

将上述命令和第二节的数据读取函数整合到一个R脚本中，运行脚本即可得到如图1风廓线图。

R语言的绘图函数很多，除了文中绘制风廓线图外还能绘制：直方图、散点图、饼图等。在安装绘图包后还能绘制3D图等。其绘图功能，能满足大部分气象资料统计分析的出图需求。

3 结语和讨论

免费的开源软件R语言，能够用于数据分析和图形显示。文中使用其对风廓线雷达资料进行了批量读取、计算和绘图。

使用read.table函数简洁的实现资料读取。

R语言的操作系统函数，能便利的实现资料的批处理。

R语言的各种绘图函数，能迅速地绘制较为美观的风廓线图。

R语言的统计计算功能十分强大，本文只是简单的使用基础计算，在随后的工作中，可以使用R语言实现，雷达数据的质量控制，统计分析等。

本文只分析了雷达产品数据，R语言也能处理二进制数据资料，以后的工作中可以利用R语言分析其他类型的气象资料。

参考文献

[1] Green，J.L.，Atmospheric measurements by VHF pulsed Doppler radar.IEEE Trans.On Geoscience electronics.GE-17（4）：262-280.1979

[2] 何平.相控阵风廓线雷达[M].气象出版社，2006：104—122.

[3] 王欣，卞林根，彭浩，李剑东.风廓线仪系统探测试验与应用[J].应用气象学报，2005，16（5）：693-698.

[4] 胡明宝.风廓线雷达数据处理和应用研究[D].南京信息工程大学博士论文，2012.

[5] John M.Chambers.Facets of R.The R Journal，1（1）：5-8，2009.

[6] 关于进行风廓线雷达数据传输的通知.气测函[2011]61号.内部资料endprint

m_data<-r_data[complete.cases（r_data），]

其中complete.case函数能够获取数据中不含缺测的所有列，进而赋值后的数据框m_data中剔除了r_data的缺测值。

1.2 批量数据文件的读取

通过1.1中几个简单的语句即可完成单个雷达产品数据读取和初始化，而在业务运行中，雷达产品数据肯定是大批量的生成，也需要程序脚本具有批处理功能。

我们将1.1中单个文件的雷达产品数据处理过程，整合成一个R语言函数readradar：

readradar <- function（dir，fname）{… …}

readradar只需要给定路径和文件名，即可读取雷达产品数据资料并去除缺测值。函数最后返还一个数据框，其中包含该文件中所有非缺测产品数据。给定所有需要处理的路径和文件名后，即可完成资料的批处理。

而R语言能够十分便利的获取文件名，因为其能通过脚本语言进入当前运行的操作系统。简单的说，R语言能够在程序内部完成操作系统的文件、文件夹处理、安装包的安装等。例如可以直接使用file.create（create.dir）函数直接在R语言中生成新的文件（文件夹）等。 本文只需要R语言读取要处理的文件名：

fname<-list.files（path="./Qingdao/"

fname中存储了所有“./Qingdao”下的雷达数据文件名。通过循环即可完成雷达数据的批量读取：

for （i in 1：length（fname）） {

m_data<-readradar（dir，fname[i]）

}

别的批量处理过程（资料简要处理，雷达风廓线图等）和读取类似，只需要加入循环即可。

2 雷达资料的简要处理和画图

第二节中给出了R语言对雷达产品数据文件的批读取。R语言最为实用的优点在于其计算和画图功能。

2.1 雷达风廓线的简要计算

为了便于处理资料。我们利用names函数将m_data数据框的每一列分别命名：

names（m_data）<-c（"hgt"，"h_dir"，"h_speed"，"w_speed"，"h_rel"，"w_rel"，"cn2"）

每一列表示对应的资料的采样高度，水平风向，水平风速，垂直风速，水平方向可信度，垂直方向可信度，Cn2。

水平风的u分量计算如下：

m_data$u <- m_data$h_speed * cos（（270 - m_data$h_dir）* pi/180）

其中pi为R语言中自带的圆周率的取值3.141593。而计算出来的u风速可以直接存储到数据框m_data的新的一列中。

同理可以计算出水平风速的v分量：

m_data$v <- m_data$h_speed * sin（（270 - m_data$h_dir）* pi/180）

上述计算过程可以看出，R语言中数据框的计算处理十分方便，不需要使用循环语句计算不同高度上的u/v分量，而且计算结果可以直接作为数据框新的一列存到原数据框中。而R语言除了基本的计算外，还能使用内部函数做各种复杂数学计算，例如矩阵的求逆、线性回归分析、抽样分布、显著性试验等。如果R语言自带的函数已经不能解决当前数学问题，还能上网搜索下载对应函数包。

2.2 雷达风廓线图

R语言主要有三种绘图函数：高级、低级和交互式。通过调用高级绘图函数，能在R语言中直接绘制各种统计图；低级绘图函数能够对现有的图进行修改；交互式绘图能够让用户直接利用鼠标修改图形。大量的内置函数让绘图变得十分的简易。例如本文需要的u风场随高度变化图可以通过以下函数实现：

plot（m_data$u，m_data$hgt）

在plot函数中增加各种参数，能够优化所绘图形。例如本文中使用如下命令获取风廓线图：

plot（m_data$u，m_data$hgt，type="b"，

main=title，xlim=c（-20，20），ylim=c（0，5000），

xlab="风速（m/s）"，ylab="采样高度（m）"，pch=16，col=2）

其中tpye=“b”表示曲线为点画线；main为主标题；xlim和ylim表示横纵坐标取值范围；xlab、ylab表示横、纵坐标标题；pch=16表示点为实心圆圈，col=12表示颜色为红色。

通过使用低级绘图命令如points、lines等，可以在刚画的u风场廓线后增加v风场廓线：

points（m_data$v，m_data$hgt）

lines（m_data$v，m_data$hgt）

图例则可以用lengend函数添加：

legend（"topleft"，pch=c（16，17），col=c（2，19），lty=c（1，1），legend=c（"u"，"v"），bty="n"）

R语言还能使用par函数对图片进行设置，例如本文中使用如下函数将2014年4月26日每隔3小时的风廓线显示到同一图片中：

par（mfrow=c（2，4），mar=c（5， 4， 4， 2） + 0.1，oma=c（0.1，0.1，0.4，0.1））

其中mfrow=c（2，4）表示图片分割为两行四列，参数mar和oma指定了图片的间距。

将上述命令和第二节的数据读取函数整合到一个R脚本中，运行脚本即可得到如图1风廓线图。

R语言的绘图函数很多，除了文中绘制风廓线图外还能绘制：直方图、散点图、饼图等。在安装绘图包后还能绘制3D图等。其绘图功能，能满足大部分气象资料统计分析的出图需求。

3 结语和讨论

免费的开源软件R语言，能够用于数据分析和图形显示。文中使用其对风廓线雷达资料进行了批量读取、计算和绘图。

使用read.table函数简洁的实现资料读取。

R语言的操作系统函数，能便利的实现资料的批处理。

R语言的各种绘图函数，能迅速地绘制较为美观的风廓线图。

R语言的统计计算功能十分强大，本文只是简单的使用基础计算，在随后的工作中，可以使用R语言实现，雷达数据的质量控制，统计分析等。

本文只分析了雷达产品数据，R语言也能处理二进制数据资料，以后的工作中可以利用R语言分析其他类型的气象资料。

参考文献

[1] Green，J.L.，Atmospheric measurements by VHF pulsed Doppler radar.IEEE Trans.On Geoscience electronics.GE-17（4）：262-280.1979

[2] 何平.相控阵风廓线雷达[M].气象出版社，2006：104—122.

[3] 王欣，卞林根，彭浩，李剑东.风廓线仪系统探测试验与应用[J].应用气象学报，2005，16（5）：693-698.

[4] 胡明宝.风廓线雷达数据处理和应用研究[D].南京信息工程大学博士论文，2012.

[5] John M.Chambers.Facets of R.The R Journal，1（1）：5-8，2009.

[6] 关于进行风廓线雷达数据传输的通知.气测函[2011]61号.内部资料endprint