一种地理栅格数据任选多边形区域的直方图实时生成方法
【技术领域】
[0001] 本发明属于地理信息处理技术领域，具体设及一种地理信息系统中地理栅格数据 任选多边形区域的直方图实时生成方法。
【背景技术】
[0002] 直方图是一个统计学概念，它通过对选定地域某一方面特征（如海拔、颜色等)进 行统计来表征该地域的全局特征。直方图广泛应用于许多领域比如图像处理、数据挖掘W 及数据库应用。
[0003] 随着地理空间数据采集技术的发展，人类获取地理空间数据的效率不断提高，用 于分析处理的地理空间数据的规模越来越大。与此同时，随着地理数据分析需求的多样化 发展，用户往往需要通过交互的方式快速得到自己选定空间区域的直方图数据。目前绝大 多数GIS软件如ArcGIS、SuperMap、MapGIS等均具有生成直方图的功能，但是都不支持任选 多边形区域的直方图生成，而且其直方图生成耗时比较长，远不能满足与用户实时交互的 需要。
[0004] 当前国内外研究直方图快速生成算法主要有两个方向。
[0005] 一个是利用CPlKCentral Processing Unit)进行直方图生成。随着技术的进步， CPU处理速度越来越快，计算的性能可W随着CPU性能和数量的增加而不断提高，相对而言 I/O性能的增长跟不上计算性能的增长。在地理栅格数据处理领域，当处理大规模数据时， 低性能的I/O将成为影响整体性能的瓶颈，运严重制约了生成直方图运算效率的提高。因 此，利用CPU进行直方图生成往往是通过减少数据的访问量、提高访问速度来提高生成速 度。
[0006] 另一个是利用GPU(Gra地ics Processing Unit)进行直方图生成。GPU作为一种专 为计算密集型、高度并行化应用而设计的高性能计算平台，其运算能力和存储器带宽上相 对于CPU有明显的优势，通过计算统一设备构架CUDA(Compute Unified Device Architecture),GPU可W在单指令多数据SIMD(Single Instruction Multiple Data)编程 模型下发挥其强大的计算能力，因此，高效的直方图生成可W依托运些并行处理器平台。基 于GPU的直方图生成算法的并行化难点在于如何减少数据分布的不规则所产生的写冲突。 目前CUDA的软件开发包已经支持32位色图和256色图两种基础的直方图生成算法，两者的 实现主要依靠于对直方图进行多次复制，即为每个线程生成一个私有直方图副本W减少写 冲突。
【发明内容】
[0007] 本发明的目的是提供一种地理栅格数据任选多边形区域的直方图实时生成方法。 利用本发明为用户提供服务的流程如图1所示，用户根据需求在客户端屏幕自由框选任意 多边形区域，客户端将用户输入的请求发送到服务器，服务器经过计算后再将求得的直方 图结果返回到客户端，最后由客户端将结果反馈给用户，而且为了满足实时交互的需求，从 用户输入完毕到得到直方图结果，耗时不应超过0.5秒。
[0008] 为实现上述目的，本发明技术解决方案如下：
[0009] -种地理栅格数据任选多边形区域的直方图实时生成方法，包括W下步骤：
[0010] 1)进行预处理。首先将原始栅格数据分割成一定尺寸的瓦片（可将瓦片大小取为 256 X 256)，其分割示意图如图3所示，图中阴影部分表示原始栅格数据，一个方格代表一个 瓦片。然后分别计算栅格数据中各瓦片的直方图并将其W文件方式存储到服务器。计算一 块瓦片的直方图的过程伪代码如下：
[0012] 下面结合伪代码对该过程进行简要说明。其中tile_sizeX和tile_sizeY为瓦片的 长和宽，tile_data为瓦片中需要统计为直方图的属性数值，X和y表示当前读取点在瓦片中 的位置坐标，interval对应直方图的间隔，tile_hist为该瓦片的直方图结果向量，i表示向 量tile_Mst的第i项。计算一块瓦片的直方图的过程为依此遍历瓦片中的每一个点，对于 其中任意一个点(X，y)，判断该点对应的瓦片属性值ti 1 e_data(X，y)与直方图间隔的关系， 如果存在i X interval < tile_data(x，y)<(i+l) X interval，则为结果直方图向量tile_ hist的第i项加1。
[0013] 2)判定各瓦片与用户选定多边形区域的关系，分别确定处于选定区域内部的瓦片 和处于选定区域边界的瓦片。如图4所示，原始栅格数据被分成多个瓦片，虚线表示选定区 域的边界。那么接下来需要判定哪些瓦片在选定区域内部哪些瓦片与选定区域边界相交。 首先计算能够将用户选定多边形区域完全包容进去的最小矩形即最小矩形外包框，将其表 示为(別，71，^2,72)，其中^1，71)和^2,72)分别为最小矩形外包框的左下角和右上角的 顶点坐标。基于客户端显示的缩放级别可W确定分辨率R，那么客户端显示的一个像素点就 代表RXR个栅格单元。从而可利用下取整函数floor得到纵向填充的范围跨度是从floor (71/1?巧化1〇〇1'(72/10+1，横向填充的范围跨度是从^〇〇1'(^/10至相1〇〇1'^2/10 + 1。然后利 用扫描线填充算法（参见Donald D.Hearn,M.Pauline Baker.Computer Graphics with 0penGL(3th Edition)[M].USA:Prentice 化11，2010.P187-P201.)在分辨率R下对选定区 域进行填充。假设一个瓦片有NXN个栅格单元，那么为每个瓦片建立一个(N/R) X (N/R)的 填充矩阵flag_tile便可W表示瓦片的所有点。在填充过程中，对于在选定区域内的点，将 flag_tile中与该点对应的位置置为1，对于选定区域外的点，flag_tile中与该点对应的位 置置为0。同时记录各填充矩阵中被置为1的点的个数。对于一个填充矩阵，如果矩阵中的每 一个点都被置为1，那么可W判断运个填充矩阵对应的瓦片完全在选定区域内。那些未全部 被置为1的填充矩阵对应的瓦片处在选定区域边界。运些填充矩阵将被用于确定原始栅格 数据中的点的位置来计算直方图数据。
[0014] 3)对于步骤2)确定的处于选定区域内部的瓦片，读入步骤I)中生成的文件，读取 对应瓦片的直方图数据并将各瓦片的直方图数据汇总求和得直方图向量inner_hist，其计 算公式如下：
[0015] inner_hist(i) = 2tiie_xeinner_tiietile_hist_of_tile_x(i)
[0016] 其中inne;r_tile为处于选定区域内部的瓦片集合，i表示innerjiist的第i项， 1：；[16_11131：_0；1^_1：；[16_义为瓦片1：；[16_义对应的直方图向量。
[0017] 4)对处于选定区域边界的瓦片，采用并行方式读取原始栅格数据并统计直方图。
[0018] 首先根据步骤2)确定的处于选定区域边界的瓦片，并根据瓦片对应的填充矩阵确 定瓦片处于选定区域内的部分。将每个处于选定区域边界的瓦片作为一个子任务。如图5所 示，程序的主进程PO将子任务分配给不同的工作进程，将处在选定区域边界的瓦片对应的 填充矩阵flag_tile发送给不同的工作进程，各工作进程分别使用开源栅格空间数据转换 库GDAL(参见GDAUGeospatial Data Abstraction Library .http://www.gdal .org/)在原 始栅格数据中提取瓦片并根据瓦片对应的填充矩阵来确定瓦片处于选定区域内的部分，然 后统计运部分的直方图。
[0019] 对每一个处在选定区域边界的瓦片进行直方图统计的过程伪代码如下：
[0021]下面结合伪代码对该过程进行简要说明。其中tile_sizeX和tile_sizeY为瓦片的 长和宽，R为分辨率因子，flag_tile为该瓦片对应的填充矩阵，X和Y用于表示flag_tile中 的位置坐标，tile_da化为瓦片中需要统计为直方图的属性数值，X和y表示当前读取点在瓦 片中的位置坐标，interva 1对应直方图的间隔，bound_ti 1 e_hist为该瓦片的直方图结果向 量，i表示向量6〇11扣1_^16_1113*的第1项。对每一个处在选定区域边界的瓦片进行直方图统 计的步骤为:遍历该瓦片对应填充矩阵flag_tile中的每一个点，对于其中任意一个点(X， Y)如果flag_tile(X，Y) = l，则表示与该点对应的瓦片上的点处在选定区域内部，即瓦片上 横向范围从X X R到(X+1) X R-1纵向范围从Y X R到(Y+1) X R-1的所有点均处于选定区域内 部。然后，遍历瓦片上运部分区域内的点，对于任意一个点(x，y)，判断该点对应的瓦片属性 值tile_data(x，y)与直方图间隔的关系，如果存在i X interval < tile_data(x，y)<(i+l) X interval，则为结果直方图向量bound_tile_hist的第i项加1。
[0022] 然后将各工作进程求得的结果进行求和得直方图向量boundjiist，计算公式如 下，其中bound_tile为处于选定区域边界的瓦片集合，i表示第i项，bound_tile_hist_of_ tile_x为边界处瓦片tile_x对应的直方图向量。
[0023] bound_h i S t (i) = 2 ti ie_xeb〇imd_ti iebound_t i 1 e_h i S t_o f _t i 1 e_x (i)
[0024] 5)将步骤3)得到的直方图向量inne;r_hist和步骤4)中得到的直方图向量bound_ hist求和得到选定区域的直方图histogram,计算公式如下，其中i表示向量histogram的第 i项：
[0025] histogram(i)=inner_hist(i)+bound_hist(i)
[0026] 本发明的有益效果是：
[0027] 通过对原始栅格数据进行一定的预处理，减少了计算直方图时对原始栅格数据的 访问量，并结合客户端屏幕分辨率对生成方法进行了优化，同时利用高性能集群进行并行 化，最终可W满足实时获取用户选定任意多边形区域直方图的需求。
【附图说明】
[0028] 图1是本发明为用户提供服务的流程示意图
[0029] 图2是本发明的流程示意图
[0030] 图3是栅格数据分割为瓦片示意图
[0031] 图4是选定区域内瓦片的判定示意图
[0032] 图5是使用GDAL在并行I/O模式下处理地理空间栅格数据的示意图
[0033] 图6是本发明方法与直接从原始栅格进行统计的直方图生成方法的任务耗时对比 (两种方法均采用256个工作进程进行并行化）
【具体实施方式】
[0034] 下面结合附图和具体实例对本发明做进一步描述。
[0035] 图1为本发明为用户提供服务的流程示意图。用户根据需求在客户端屏幕自由框 选任意多边形区域，客户端将用户输入的请求发送到服务器，服务器经过计算后再将直方 图结果返回到客户端，最后由客户端将直方图结果展示给用户。
[