本发明属于雷达数据分析处理技术领域，具体涉及基于GIS的X波段面雨量雷达数据提取方法。
背景技术：
X波段面雨量雷达测雨系统能够解决目前离散雨量站网雨量监测数据点、面关系不能完全匹配问题，而当前水文专业领域尚缺乏较为完善的天气雷达监测体系。因此，如何将雷达数据融合到现有的水文监测平台中，形成天地一体的水文资源监测体系成为急需解决的问题。其中，雷达数据如何进行科学、完整的提取是关键难题。
技术实现要素：
为解决上述问题，本发明公开了一种雨量雷达数据提取方法，利用GIS技术将极坐标投影的雷达数据转化为格网化的空间数据，并根据给出的点、面提取对应的雨量值，为下一步准确水文分析提供数据来源。
为了达到上述目的，本发明提供如下技术方案：
基于GIS的X波段面雨量雷达数据提取方法包括如下步骤：
步骤一，数据的存储
以雷达所在位置为圆心，按极坐标方式存储，共360个径向，从0度至359度排列，每径向数据存储若干雨强值，每个雨强值占用两个字节；
步骤二，极坐标数据转换为栅格数据
步骤二-1，生成栅格数据
(1)生成空栅格
将雷达所在位置的经纬度通过投影转换为坐标，调用函数生成空栅格数据；
(2)极坐标数据写入到空栅格
①定义一个二维数组，数组中存储有极坐标方位角angel，库数n，利用二重循环将雷达极坐标数据写入到二维数组中；
②从用来存储极坐标数据的空栅格抽象获得像素块，用于临时存储雨量数据；
③通过下式计算栅格点到雷达中心的距离：
式中dis为栅格点到雷达中心距离，x，y为栅格点坐标；
④通过下式计算栅格点到雷达中心的方位角：
式中azim为方位角，θ＝arctan(y/x)；
⑤根据栅格点到雷达中心的距离dis、方位角azim在极坐标中找到相应的雨量值，从而获得该栅格点的雨量值，从像素块左上角开始，先X方向在Y方向，依次遍历，赋予栅格点雨量值，将像素块重新写回空栅格，完成雷达极坐标数据到栅格数据的转换。
进一步的，还包括以下步骤：
步骤三，雷达雨量图像展示
步骤三-1，对栅格图进行分级渲染；
步骤三-2，利用色阶渲染形成雨量颜色分级；
步骤三-3，对雨量图进行透明度设置。
进一步的，还包括步骤四，雷达雨量数据提取：用点、面批量化方式提取雷达雨量值，根据站点和/或行政边界和/或小流域提取。
与现有技术相比，本发明具有如下优点和有益效果：
本发明提供的基于GIS的X波段面雨量雷达数据提取方法，将极坐标投影的雷达数据转化为格网化的空间数据，实现了雷达极坐标数据到面雨量数据的转换、展示，并能够根据给出的点、面提取对应的雨量值，为下一步准确水文分析提供数据来源，为雷达雨量在常规水文降雨监测，城市洪涝减灾、中小流域山洪泥石流灾害预警、水资源管理中的应用打下基础。
附图说明
图1为雷达雨量图
具体实施方式
以下将结合具体实施例对本发明提供的技术方案进行详细说明，应理解下述具体实施方式仅用于说明本发明而不用于限制本发明的范围。
本发明提出的基于GIS的X波段面雨量雷达数据提取方法包括如下步骤：
步骤一，数据的存储
X波段测雨雷达数据有严谨的数据结构和时空排列顺序，其包含文件头和雨强数据两部分。文件头数据长度为256个字节，包含数据标识、数据版本号、雷达站址和探测参数等信息，保留字节全部填充0。
表1列出了文件头数据的数据结构。
表1文件头数据结构
数据的存储方式如下：记录方式以雷达所在位置为圆心，按极坐标方式存储，共360个径向，从0度至359度排列，每径向数据存储600(本发明中总库数为600个)个雨强值，每个雨强值占用两个字节。数据存储结构如表2所示，其中n为库数。
表2径向数据存储结构
步骤二，极坐标数据转换为栅格数据
本发明将雨量雷达极坐标数据巧妙地转为1200*1200的二维矩阵。
步骤二-1，生成栅格数据
坐标系统是GIS数据重要的数学基础，用于表示地理要素、图像和观测结果的参照系统，坐标系统的定义能够保证雷达数据在软件中正确的显示其位置、方向和距离。本发明中采用西安1980投影。
(1)生成空栅格
雷达所在位置的经纬度在数据头文件中已知，通过投影转换后设为(centx,centy)，单位为米。利用C#语言调用ArcEngine中IRasterWorkspace2接口的CreateRasterDataset函数可以生成空栅格。该函数的主要参数如表3所示：
表3生成空栅格数据主要参数说明
(2)极坐标数据写入到空栅格
①定义一个二维数组rain[angle,n]，angel为极坐标方位角，n为库数，利用二重循环将雷达极坐标数据写入到rain[angle,n]中。
②从用来存储极坐标数据的空栅格抽象获得PixelBlock块(像素块，实际为二维灰度矩阵)，其主要作用是临时存储雨量数据。
③计算栅格点到雷达中心的距离，如公式2-1所示：
式中dis为栅格点到雷达中心距离，x，y为栅格点坐标。
④计算栅格点到雷达中心的方位角，如公式2-2所示：
式中azim为方位角，θ＝arctan(y/x)。
⑤栅格点到雷达中心的距离dis、方位角azim对应维数组rain[angle,n]中angle，n两个数值，从而获得该栅格点的雨量值。从PixelBlock左上角开始，先X方向在Y方向，依次遍历，赋予栅格点雨量值。将PixelBlock重新写回空栅格，完成雷达极坐标数据到栅格数据的转换。
具体的说，得到栅格点到雷达中心的方位角azim后，找到极坐标中与方位角相同的结束角，在该结束角下存储有600个雨量值，由于距离dis＝库数*step，因此库数＝dis/step，根据栅格点到雷达中心的距离dis能够得到相应的库数，也就能找到相应库数下存储的雨量值，即获得该栅格点的雨量值。
在将雷达数据转换为栅格数据后，还包括雷达雨量图像展示的步骤。
步骤三，雷达雨量图像展示
步骤三-1，栅格雨量图渲染
为了方便的从雷达面雨量图上区分雨量的量值、位置，需要对栅格图进行分级渲染。ArcEngine提供IRasterRenderer接口对栅格图层进行分级渲染。
步骤三-2，雨量色阶
利用色阶渲染形成雨量颜色分级，以便用户识别不同颜色对应的雨量值。
步骤三-3，图层透明度设置
避免雨量图覆盖水文GIS地图中的边界、流域信息，还需要对雨量图进行一定的透明度设置。
通过以上步骤，最终雷达雨量图在水文GIS地图中的展示如图1所示。
步骤四，雷达雨量数据提取
为了对雨量雷达的雨量值进行验证，需要用已有的遥测站数据进行比较；小流域的面雨量提取对中小河流洪水预报的准确性起到关键作用。用点、面批量化方式提取雷达雨量值，可以根据站点和/或行政边界和/或小流域提取。
本发明方案所公开的技术手段不仅限于上述实施方式所公开的技术手段，还包括由以上技术特征任意组合所组成的技术方案。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也视为本发明的保护范围。