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

1、基于R语言的风廓线雷达产品数据简单处理基于R语言的风廓线雷达产品数据简单处理摘要：风廓线雷达产品数据在气象上应用广泛，对其产品数据资料的分析、处理和可视化十分重要。本文利用免费的开源软件R语言，对雷达资料产品数据进行批量读取、计算和绘图。展示了R语言在处理风廓线雷达资料的特点：简洁、易学。关键词：R语言 风廓线雷达产品数据 风廓线图风廓线雷达作为中小尺度天气系统有效的探测工具,能够24小时不间断的提供:水平风向、风速、垂直气流等气象要素随高度的分布,是进行边界层和高空气象探测的重要设备,能对现有的气象观测进行补充。其观测资料在气象上应用广泛:对高影响天气的监测和预警；对灾害天气系统的天气学分析

2、；用于数值预报中的观测资料同化和数据后处理等。因此开展风廓线雷达的处理十分必要14。现有的气象数据处理一般使用Fortran,C等高级编程语言。虽然有着处理速度快,计算结果可靠等优点。但是其使用并不便:需要掌握编程语法和算法,不同计算机平台可能需要不同的编译系统等。所以本文针对风廓线雷达资料,使用R语言来对其进行简单处理。R语言的全称:为统计计算和图形展示而设计的一种编程语言和统计环境。目前有一个R核心开发团队对其进行定期维护和更新5。其主要的特点在于:(1)完全免费:可以在其官方网站(http:/cran.r-project.org)下载到完整的安装包并免费使用；(2)开源软件:R语言源代码

3、完全公开,任何人均能提供各种数据处理模块,下载相应模块后能迅速地完成数据处理等工作；(3)图形显示:可以利用R语言处理数据后,直接获取各种统计分析图形；(4)多平台使用:R语言可以在Windows、mac、Unix这些操作系统上安装,不需要重复编译,其脚本在不同操作系统间可以任意使用。针对风廓线雷达资料,我们可以利用R语言对其进行简要的数据处理,统计分析,然后利用图形显示系统得到各种分析图片。本文将在第二部分介绍R语言读取和储存风廓线雷达产品数据资料；文章第三部分对资料进行简要处理,最后会利用R语言的图形显示功能作图。1 雷达产品数据的读取测站观测的风廓线雷达资料有两种数据6:原始数据和产品数

4、据。其中原始数据为二进制格式,主要包含功率谱数据文件,瞬时径向谱数据文件。产品数据文件为文本格式,包括实时的采样高度上的、半小时平均的采样高度上的、一小时平均的采样高度上的产品数据文件。本文所处理资料为实时的采样高度上的产品数据文件。R语言也能读取处理二进制格式文件,对雷达原始数据资料的处理会在以后的工作中进行。1.1 单个产品数据的读取R语言有许多函数能够直接从文本文件中读取数据,比较常用的有:read.table(),read.csv(),read.fwf()。其中和Fortran比较接近的是read.fwf(),可以指定读取数据的长度和格式。于Fortran不同的是,R语言在读取数据的时

5、候,不用先给定数据类型,程序会直接读取数据,并存储到一个数据框(data.frame)里。例如针对本文要处理的雷达产品数据,可以直接使用read.table()命令读取产品数据:raw_data<-read.table(fname,fill=TRUE)其中raw_data为一个数据框(data.frame),用于储存我们需要处理的雷达产品数据,fname 为要读取的雷达数据文件名字；“headerFALSE”表示该实体数据中没有数据说明头文件；由于本文要处理的雷达产品数据不是规则的表格形式,所以需要使用参数:“fillTRUE”,来自动填满不是数据表格的部分。产品数据文件的前三行为测站基

6、本参数,最后一行为结束行,中间部分为实际数据,包含:采样高度,水平风向,水平风速,垂直风速,水平方向可信度,垂直方向可信度,Cn2。为了便于资料处理,将实体数据单独储存到一个名为r_data的新数据框中:r_data<-raw_data4:(length(raw_data,1)-1),其中使用“length(raw_data,1)”函数判断raw_datpr_data<-lapply(r_data,as.numeric)r_data<-as.data.frame(r_data)先使用函数“lapply”和“as.character”把每列元素转换为字符型,再使用as.nume

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

8、中的数据:m_data<-r_datacomplete.cases(r_data),其中complete.case函数能够获取数据中不含缺测的所有列,进而赋值后的数据框m_data中剔除了r_data的缺测值。1.2 批量数据文件的读取通过1.1中几个简单的语句即可完成单个雷达产品数据读取和初始化,而在业务运行中,雷达产品数据肯定是大批量的生成,也需要程序脚本具有批处理功能。我们将1.1中单个文件的雷达产品数据处理过程,整合成一个R语言函数readradar:readradar <- function(dir,fname) readradar只需要给定路径和文件名,即可读取雷达产品

9、数据资料并去除缺测值。函数最后返还一个数据框,其中包含该文件中所有非缺测产品数据。给定所有需要处理的路径和文件名后,即可完成资料的批处理。而R语言能够十分便利的获取文件名,因为其能通过脚本语言进入当前运行的操作系统。简单的说,R语言能够在程序内部完成操作系统的文件、文件夹处理、安装包的安装等。例如可以直接使用file.create(create.dir)函数直接在R语言中生成新的文件(文件夹)等。 本文只需要R语言读取要处理的文件名:fname<-list.files(path=“./Qingdao/”fname中存储了所有“./Qingdao”下的雷达数据文件名。通过循环即可完成雷达数

10、据的批量读取:for (i in 1:length(fname) m_data<-readradar(dir,fnamei)别的批量处理过程(资料简要处理,雷达风廓线图等)和读取类似,只需要加入循环即可。2 雷达资料的简要处理和画图第二节中给出了R语言对雷达产品数据文件的批读取。R语言最为实用的优点在于其计算和画图功能。2.1 雷达风廓线的简要计算为了便于处理资料。我们利用names函数将m_data数据框的每一列分别命名:names(m_data)<-c(“hgt”,”h_dir”,”h_speed”,”w_speed”,”h_rel”,”w_rel”,”cn2”)每一列表示对应

11、的资料的采样高度,水平风向,水平风速,垂直风速,水平方向可信度,垂直方向可信度,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语言中数据框的计算处理十分方便,不需要使用循环语句计算不同

12、高度上的u/v分量,而且计算结果可以直接作为数据框新的一列存到原数据框中。而R语言除了基本的计算外,还能使用内部函数做各种复杂数学计算,例如矩阵的求逆、线性回归分析、抽样分布、显著性试验等。如果R语言自带的函数已经不能解决当前数学问题,还能上网搜索下载对应函数包。2.2 雷达风廓线图R语言主要有三种绘图函数:高级、低级和交互式。通过调用高级绘图函数,能在R语言中直接绘制各种统计图；低级绘图函数能够对现有的图进行修改；交互式绘图能够让用户直接利用鼠标修改图形。大量的内置函数让绘图变得十分的简易。例如本文需要的u风场随高度变化图可以通过以下函数实现:plot(m_data$u,m_data$hgt

13、)在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风场廓

14、线: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)表示图片分割为两行四列

15、,参数mar和oma指定了图片的间距。将上述命令和第二节的数据读取函数整合到一个R脚本中,运行脚本即可得到如图1风廓线图。R语言的绘图函数很多,除了文中绘制风廓线图外还能绘制:直方图、散点图、饼图等。在安装绘图包后还能绘制3D图等。其绘图功能,能满足大部分气象资料统计分析的出图需求。3 结语和讨论免费的开源软件R语言,能够用于数据分析和图形显示。文中使用其对风廓线雷达资料进行了批量读取、计算和绘图。使用read.table函数简洁的实现资料读取。R语言的操作系统函数,能便利的实现资料的批处理。R语言的各种绘图函数,能迅速地绘制较为美观的风廓线图。R语言的统计计算功能十分强大,本文只是简单的使用

16、基础计算,在随后的工作中,可以使用R语言实现,雷达数据的质量控制,统计分析等。本文只分析了雷达产品数据,R语言也能处理二进制数据资料,以后的工作中可以利用R语言分析其他类型的气象资料。参考文献1 Green,J.L.,Atmospheric measurements by VHF pulsed Doppler radar.IEEE Trans.On Geoscience electronics.GE-17(4):262-280.19792 何平.相控阵风廓线雷达M.气象出版社,2006:104122.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 关于进行风廓线雷达数据传输的通知.气测函201161号.内部资料