一种面向虚拟地球的栅格数据快速建库方法.pdf

摘要
申请专利号：	CN201110366702.8	申请日：	2011.11.18
公开号：	CN102509022A	公开日：	2012.06.20
当前法律状态：	终止	有效性：	无权
法律详情：	未缴年费专利权终止IPC(主分类):G06F 19/00申请日:20111118授权公告日:20150218终止日期:20151118\|\|\|授权\|\|\|实质审查的生效IPC(主分类):G06F 19/00申请日:20111118\|\|\|公开
IPC分类号：	G06F19/00(2011.01)I	主分类号：	G06F19/00
申请人：	武汉大学
发明人：	向隆刚; 曾珍
地址：	430072 湖北省武汉市武昌区珞珈山武汉大学
优先权：
专利代理机构：	武汉科皓知识产权代理事务所(特殊普通合伙) 42222	代理人：	薛玲
PDF下载：	PDF下载

内容摘要

本发明涉及一种面向虚拟地球的栅格数据快速建库方法。该方法提出，如果第0层一个瓦片对应的原始遥感影像数据无法一次完全载入内存，从中间层（内存需求最小层）开始构建金字塔。这样打破了以往逐层构建的局限性，仅需扫描原始数据文件一次，减少了I/O操作次数，消除了处理海量数据时读取数据的瓶颈问题。本发明还设计了瓦片处理线程和瓦片计算线程两类线程，让两类线程同时运行，充分利用了系统资源，根据硬件不同的条件，配置多个瓦片计算线程，使程序读数据和写数据的速度相匹配，实现了读写操作的同步性，节省了读写操作相互等待的时间。提高了金字塔构建的效率，尤其适合于对超大遥感数据文件的处理。

权利要求书

1：一种面向虚拟地球的栅格数据快速建库方法，其特征在于，包括以下步骤：步骤 1，读取原始遥感影像的信息，包括原始遥感影像的波段数 bn、原始遥感影像的尺寸 pw×ph、原始遥感影像的影像分辨率 pd 以及原始遥感影像对应的地理范围（rw, rs, re, rn）；定义瓦片结构体，所述瓦片结构体是用于在内存中存储瓦片数据并对瓦片进行标识的数据结构体；步骤 2，待构建金字塔的最顶层为分辨率最低的层次，标记为第 0 层；根据原始遥感影像对应的地理范围（rw, rs, re, rn）和第 0 层瓦片分辨率 yd[0]，计算出第 0 层瓦片的起止范围，计算方式如下， sr[0] = (rw+180+1e -8)/yd[0] er[0] = (re+180+1e-8)/yd[0] sc[0] = (rs+90+1e-8)/yd[0] ec[0] = (rn+90+1e-8)/yd[0] 其中， sr[0] 为第 0 层的开始行、 er[0] 为第 0 层的结束行、 sc[0] 为第 0 层的开始列、 ec[0] 为第 0 层的结束列、 e 为数学参数；第 0 层瓦片分辨率 yd[0] 根据原始遥感影像的影像分辨率 pd 计算，计算公式如下， l -1 yd[0]=pd×tilesize×2 其中， l 为预设的待构建金字塔的层数， tilesize 为预设的瓦片大小；步骤 3，比较第 0 层一个瓦片对应原始遥感影像所需的内存 memvolume[0] 与预设的可申请的内存，设预设的可申请内存上限为 ms 个字节，当所需内存大小 memvolume[0] 不大于可申请内存上限 ms 时进入步骤 4，否则进入步骤 5 ；所述第 0 层一个瓦片对应原始遥感影像所需的内存大小 memvolume[0] 根据原始遥感影像的波段数 bn 计算，计算公式如下， memvolume[0] =pixelvolume[0]×bn 其中， pixelvolume[0] 为第 0 层中存储一个瓦片覆盖原始遥感影像的其中一个波段的信息所需的内存大小，根据第 0 层瓦片分辨率 yd[0]、原始遥感影像的尺寸 pw×ph 和原始遥感影像对应的地理范围（rw, rs, re, rn）计算，计算公式如下， pixelvolume[0] = (yd[0]×pw/(re-rw)+0.5)×( yd[0]×ph/(rn-rs)+0.5) ×sizeof(type) 其中 type 为存储像素值的数据类型， sizeof（）是获取数据类型在计算机内占有字节数的函数；步骤 4，逐个处理第 0 层的每一个瓦片，处理方式为计算第 0 层中某瓦片在原始遥感影像的像素覆盖范围，并将该范围对应的原始影像数据的数据块读入到内存，记为数据块 i，数据块 i 对应的第 0 层瓦片及其四分后代瓦片的瓦片结构体记为集合 sti，将集合 sti 初始化；遍历集合 sti，当遍历到集合 sti 中的某个瓦片结构体时，在数据块 i 中获取相应的数据，通过插值获得瓦片数据并编码；处理第 0 层的所有瓦片后进入步骤 6 ；步骤 5，设某一层中存储一个瓦片对应的原始影像所需的内存大小为 memvolume[level]，存储从第 0 层到该层上一层的所有瓦片结构体所需的内存大小是 2 tlstuvolume[level]，设 memvolume[level]+tlstuvolume[level] 取值最小的层为内存需求最小层 wl ；进行以下步骤，步骤 5.1，求取内存需求最小层 wl，计算公式如下， wl=max(0, └ (l -e -1)/2 ┘ ) 其中 e = └ log4(m /3) ┘， └┘代表向下取整； m 表示第 0 层瓦片个数，根据步骤 2 所得第 0 层的开始行 sr[0]、第 0 层的结束行 er[0]、第 0 层的开始列 sc[0] 和第 0 层的结束列 ec[0] 计算，计算公式如下， m=(er[0]-sr[0]+1)*(ec[0]-sc[0]+1) 步骤 5.2，将第 0 层到第 wl-1 层的所有瓦片结构体记为集合 ts，将集合 ts 初始化，并加载到内存之中；步骤 5.3，逐个处理第 w1 层的每一个瓦片，处理方式为计算第 w1 层中某瓦片在原始遥感影像的像素覆盖范围，并将该范围对应的原始影像数据的数据块读入内存，记为数据块 i，数据块 i 对应的第 w1 层瓦片及其四分后代瓦片的瓦片结构体记为集合 sti，将集合 sti 初始化；遍历集合 sti 和集合 ts，当遍历到集合 sti 中的某个瓦片结构体时，在数据块 i 中获取相应的数据，通过插值获得瓦片数据并编码；当遍历到集合 ts 中的某个瓦片结构体时，判断瓦片的覆盖范围是否包括数据块 i，如果包括则插值获得该瓦片内相应部分区域的瓦片数据，在获得该瓦片内所有瓦片数据后进行编码；处理第 w1 层的所有瓦片后进入步骤 6 ；步骤 6，输出多分辨率金字塔。
2：如权利要求 1 所述面向虚拟地球的栅格数据快速建库方法，其特征在于：步骤 4 和步骤 5 中，设置瓦片处理线程和瓦片计算线程两类线程，并开辟两块内存区间主缓存和备用缓存来存储数据块；瓦片处理线程将原始遥感影像的数据块读入备用缓存，并在瓦片计算线程遍历完瓦片结构体后，将备用缓存设置为主缓存，并继续将下一个数据块读入备用缓存；瓦片计算线程读取主缓存内的数据块，插值计算瓦片数据。
3：如权利要求 1 或 2 所述面向虚拟地球的栅格数据快速建库方法，其特征在于：瓦片结构体包括变量 count、 status、 level、 row、 col、 pixelSpan 和指针 pTileBmp，其中变量 count 标识瓦片中有效值的个数、变量 status 标识瓦片状态、变量 level 标识瓦片所在的层数、变量 row 标识瓦片的行号、变量 col 标识瓦片的列号、变量 pixelSpan 标识瓦片中每个像素的分辨率，指针 pTileBmp 指向瓦片数据所在内存位置。

说明书

一种面向虚拟地球的栅格数据快速建库方法
    技术领域本发明涉及空间数据组织技术领域，尤其是涉及一种面向虚拟地球的栅格数据快速建库方法。
     背景技术如今越来越多的空间数据以栅格的形式存储（包括地形数据和影像数据）。栅格数据成为了重要的数据源。随着对地观测技术的快速发展，人类能够快速地获取多时相、多波段、多分辨率的海量遥感影像数据。随着分辨率的不断提高，全球数据的数据量也在以惊人的速度增加，如仅福建省的 30m×30m 分辨率的 24 位 bmp 格式的遥感影像就有 700MB，若为 1m×1m 分辨率，数据量将达到 630GB。如何去存储、组织和管理栅格数据，如何高效地检索、显示用户感兴趣的数据，使栅格数据能够更好地为实际应用服务越来越受人们的关注。
     有关文献：高峰，安培浚 . 国际空间和对地观测技术发展战略新动向 [J] . 遥感技术与应用， 2008，(06) ： 686-696.; 张更新，张昭，朱江 . 全球对地观测卫星现状及其产业发展综述 [J]. 数字通信世界，2009， (10) ： 72-75. ；王恩泉 . 中国版 Google Earth 的空间数据组织与管理研究 [D]. 北京：中国测绘科学研究院， 2007.。
     随着 Google Earth 的发布与流行，虚拟地球技术进入了人们的视野。虚拟地球能够模拟显示大范围高分辨率的海量遥感影像数据，不仅用于各个行业领域内的专业应用研究 , 也能够满足人们足不出户即可感知整个地球。目前比较热门虚拟地球产品除了 Google Earth，还有美国航天航空局（NASA）的 WorldWind、 Bing 3D 和 Skyline 等，国内不少公司与高校开展了相关研究，建立自主版权的虚拟数字地球平台，如武汉大学测绘遥感信息工程国家重点实验室自主研制的虚拟地球系统 GeoGlobe。
     有关文献：苗放，叶成名，刘瑞，孔祥生，李康荣，徐松浦，张远红 . 新一代数字地球平台与 “数字中国” 技术体系架构探讨 [J]. 测绘科学， 2007， (06) ： 157-158. ； http:// news.3snews.net/industry/20081225/6054.shtml。
     虚拟地球通常采用多分辨率瓦片金字塔结构来组织栅格数据，以提供快速的网络调度和三维渲染。对栅格数据建库的过程也就是构建多分辨率金字塔的过程。通过构建多分辨率金字塔，在多尺度金字塔结构索引机制的配合下，可以快速获取指定分辨率和范围的栅格数据，并适合于网络传输和渲染处理，从而保证用户根据视景内的范围，在任一角度浏览的数据只是金字塔的子集。面向多分辨率瓦片金字塔的数据组织与调度是虚拟地球系统得以实现网络三维高效可视化的核心思想。如何对栅格数据快速构建多分辨率的金字塔一直是人们研究的热点。
     有关文献：龚健雅，李斌等 . 当代 GIS 的若干理论与技术 . 武汉：武汉测绘科技大学出版社 . ；吴晨琛 .GeoGlobe 中尺度空间数据集管理机制 [D]. 武汉：武汉大学， 2008.。
     早期对金字塔的研究重点主要是如何提高成像质量和传输。1979 年， Tanimoto 等人提出了均值金字塔。均值金字塔的思想简单而直观，但应用到图像的渐进传输中却能获得很好的效果，引起了人们的关注。
     有关文献：朱程辉，何勇，王金玲 . 基于小波金字塔的快速图像匹配算法 [J] . 微计算机信息，2010， (10) ： 127-128.。
     同年， K.R.Sloan 和 Tanimoto 又提出了 Reduced-sum 金字塔的结构，该结构克服了均值金字塔数据量过大的问题，将金字塔中的像素数减少到与原始图像相同。
     有关文献：刘国权，李守轩 . 基于小波图像金字塔的 SSDA 快速模板匹配算法 [J]，科技广场，2007， (11) ： 134-136.。
     但无论是均值金字塔还是 Reduced-sum 金字塔，其画面质量都不够理想，图像的混叠失真较大。于是在 1983 年， P.J.Burt 等与 J.Nao 等提出了高斯 - 拉普拉斯金字塔结构，该结构在采样的同时利用高斯型低通滤波器对图像进行平滑，从而获得良好的视觉效果。
     有关文献：丁素英 . 基于小波变换的 Laplacian 金字塔图像数据压缩 [J] . 潍坊学院学报，2009， (04) ： 34-36. ；叶勤，张小虎，王栋 . 一种基于区域生长与空间形状约束的高分辨率遥感图像道路提取方法 [J]. 遥感信息， 2010， (02) ： 25-29.。
     而后的学者对如何对大数据量的金字塔进行有效的存储和管理进行了研究。 1989 年，又有学者提出了 Difference 金字塔结构，通过将相邻层的像素相减来去除相关性，从而减少金字塔的信息熵。有关文献：欧阳平，张玉方 . 形态学开闭运算在居民地边缘检测中的应用 [J]. 测绘通报， 2009， (01) ： 40-41.。
     随着计算机技术的快速发展，计算机能提供足够的空间存储金字塔数据，但是随着分辨率的提升，即使覆盖较小区域，一个遥感数据文件的尺寸也是比较大的，可能达到上十 G，甚至上百 G，从这样的超大数据文件（1 个或者多个）构建出多分辨率瓦片金字塔，其瓶颈在于读操作。
     目前，一般采用的金字塔构建方法是逐层逐块构建金字塔，对每层金字塔都获取逐个像素的扫描行，然后将扫描行放人缓存，当缓存已满时采用邻近法替换数据。其缺点是多次扫描原始遥感文件，计算一个顶层瓦片时，涉及到原始遥感文件的较大范围， I/O 效率低下。
     发明内容针对目前在虚拟地球系统中对栅格数据构建多分辨率金字塔的方法无法高效处理大数据量的遥感影像，为了提高构建多分辨率金字塔效率，本发明提供的技术方案是面向虚拟地球的一种针对海量栅格数据构建金字塔的方法，主要是针对 tif、 img 等可以存储大数据量格式的文件。
     本发明的技术方案为一种面向虚拟地球的栅格数据快速建库方法，包括以下步骤：步骤 1，读取原始遥感影像的信息，包括原始遥感影像的波段数 bn、原始遥感影像的尺寸 pw×ph、原始遥感影像的影像分辨率 pd 以及原始遥感影像对应的地理范围（rw, rs, re, rn）；定义瓦片结构体，所述瓦片结构体是用于在内存中存储瓦片数据并对瓦片进行标识的数据结构体，步骤 2，待构建金字塔的最顶层为分辨率最低的层次，标记为第 0 层；根据原始遥感影
     像对应的地理范围（rw, rs, re, rn）和第 0 层瓦片分辨率 yd[0]，计算出第 0 层瓦片的起止范围，计算方式如下， sr[0] = (rw+180+1e-8)/yd[0] er[0] = (re+180+1e-8)/yd[0] sc[0] = (rs+90+1e-8)/yd[0] ec[0] = (rn+90+1e-8)/yd[0] 其中， sr[0] 为第 0 层的开始行、 er[0] 为第 0 层的结束行、 sc[0] 为第 0 层的开始列、 ec[0] 为第 0 层的结束列、 e 为数学参数；第 0 层瓦片分辨率 yd[0] 根据原始遥感影像的影像分辨率 pd 计算，计算公式如下， l -1 yd[0]=pd×tilesize×2 其中， l 为预设的待构建金字塔的层数， tilesize 为预设的瓦片大小；步骤 3，比较第 0 层一个瓦片对应原始遥感影像所需的内存 memvolume[0] 与预设的可申请的内存，预设的可申请内存上限为 ms 个字节，当所需内存大小 memvolume[0] 不大于可申请内存上限 ms 时进入步骤 4，否则进入步骤 5 ；所述第 0 层一个瓦片对应原始遥感影像所需的内存大小 memvolume[0] 根据原始遥感影像的波段数 bn 计算，计算公式如下， memvolume[0] =pixelvolume[0]×bn 其中， pixelvolume[0] 为第 0 层中存储一个瓦片覆盖原始遥感影像的其中一个波段的信息所需的内存大小，根据第 0 层瓦片分辨率 yd[0]、原始遥感影像的尺寸 pw×ph 和原始遥感影像对应的地理范围（rw, rs, re, rn）计算，计算公式如下， pixelvolume[0] = (yd[0]×pw/(re-rw)+0.5)×( yd[0]×ph/(rn-rs)+0.5) ×sizeof(type) 其中 type 为存储像素值的数据类型， sizeof（）是获取数据类型在计算机内占有字节数的函数；步骤 4，逐个处理第 0 层的每一个瓦片，处理方式为计算第 0 层中某瓦片在原始遥感影像的像素覆盖范围，并将该范围对应的原始影像数据的数据块读入到内存，记为数据块 i，数据块 i 对应的第 0 层瓦片及其四分后代瓦片的瓦片结构体记为集合 sti，将集合 sti 初始化；遍历集合 sti，当遍历到集合 sti 中的某个瓦片结构体时，在数据块 i 中获取相应的数据，通过插值获得瓦片数据并编码；处理第 0 层的所有瓦片后进入步骤 6 ；步骤 5，设某一层中存储一个瓦片对应的原始影像所需的内存大小为 memvolume[level]，存储从第 0 层到该层上一层的所有瓦片结构体所需的内存大小是 tlstuvolume[level]，设 memvolume[level]+tlstuvolume[level] 取值最小的层为内存需求最小层 wl ；进行以下步骤，步骤 5.1，求取内存需求最小层 wl，计算公式如下， wl=max(0, └ (l -e -1)/2 ┘ ) 其中 e = └ log4(m /3) ┘， └┘代表向下取整； m 表示第 0 层瓦片个数，根据步骤 2 所得第 0 层的开始行 sr[0]、第 0 层的结束行 er[0]、第 0 层的开始列 sc[0] 和第 0 层的结束列 ec[0] 计算，计算公式如下，m=(er[0]-sr[0]+1)×(ec[0]-sc[0]+1) 步骤 5.2，将第 0 层到第 wl-1 层的所有瓦片结构体记为集合 ts，将集合 ts 初始化，并加载到内存之中；步骤 5.3，逐个处理第 w1 层的每一个瓦片，处理方式为计算第 w1 层中某瓦片在原始遥感影像的像素覆盖范围，并将该范围对应的原始影像数据的数据块读入内存，记为数据块 i，数据块 i 对应的第 w1 层瓦片及其四分后代瓦片的瓦片结构体记为集合 sti，将集合 sti 初始化；遍历集合 sti 和集合 ts，当遍历到集合 sti 中的某个瓦片结构体时，在数据块 i 中获取相应的数据，通过插值获得瓦片数据并编码；当遍历到集合 ts 中的某个瓦片结构体时，判断瓦片的覆盖范围是否包括数据块 i，如果包括则插值获得该瓦片内相应部分区域的瓦片数据，在获得该瓦片内所有瓦片数据后进行编码；处理第 w1 层的所有瓦片后进入步骤 6 ；步骤 6，输出多分辨率金字塔。
     而且，步骤 4 和步骤 5 中，设置瓦片处理线程和瓦片计算线程两类线程，并开辟两块内存区间主缓存和备用缓存来存储数据块；瓦片处理线程将原始遥感影像的数据块读入备用缓存，并在瓦片计算线程遍历完瓦片结构体后，将备用缓存设置为主缓存，并继续将下一个数据块读入备用缓存；瓦片计算线程读取主缓存内的数据块，遍历瓦片结构体，插值计算瓦片数据。而且，瓦片结构体包括变量 count、 status、 level、 row、 col、 pixelSpan 和指针 pTileBmp，其中变量 count 标识瓦片中有效值的个数、变量 status 标识瓦片状态、变量 level 标识瓦片所在的层数、变量 row 标识瓦片的行号、变量 col 标识瓦片的列号、变量 pixelSpan 标识瓦片中每个像素的分辨率，指针 pTileBmp 指向瓦片数据所在内存位置。
     本发明的优势在于从中间层（算法对于内存的需求是最小的层）开始构建金字塔，打破了以往逐层构建的局限性，仅需扫描原始数据文件一次，减少了 I/O 操作次数，消除了处理海量数据时读取数据的瓶颈问题，这对于超大文件来说提高了金字塔构建的效率。本发明还充分利用了现代计算机的多处理器多核，可并行计算的特点；采用瓦片处理线程和瓦片计算线程两类线程同时运行，充分利用了系统资源，通过配置多个瓦片计算线程，使程序读数据和写数据的速度相匹配，实现了读写操作的同步性，节省了读写操作相互等待的时间。
     附图说明
     图 1 为本发明实施例的流程图；图 2 为本发明实施例的双线性内插示意图；图 3 为本发明实施例的瓦片时空统一编码示意图。
     图 4 为本发明实施例的数据处理流程图。
     具体实施方式为便于理解本发明技术方案起见，首先说明全球多分辨率金字塔：全球多分辨率金字塔是基于球面的一种可以无限细分 , 但又不改变形状的地球体拟合格网 , 当细分到一定程度时 , 可以达到模拟地球表面的目的。本发明中采用的是通过等经纬度格网划分的全球多分辨率金字塔模型。在栅格数据金字塔中，按照全球等经纬度格网划分，每一层被裁切成统一形状、大小的数据块，这样的一个数据块称之为瓦片。
     这种全球多分辨率金字塔非常适合用来组织和管理多分辨率的海量卫星影像数据，因为它遵循以下规律： 1）全球地理坐标的经度范围为 [-180o, +180o]，纬度范围为 [-90o， +90o] ； 2）第 k 层金字塔的分辨率为第 (k+l) 层金字塔的 2 倍； 3）在任意层内，横向和纵向金字塔瓦片的数目之比为 2:1 ； 4）金字塔瓦片的编码顺序由左到右，由下到上；根据上述规则，很容易计算出每层瓦片的个数，同样也能快速计算出给定经纬度的点落在哪一个瓦片内。例如，已知一个点的经纬度为，该点在第 k 层金字塔的行号 RowNo 和列号 ColNo 可以通过下面的公式计算出： RowNo= └ ((φ + 90) / (180/2k ) ┘ mod 2k ColNo = └ ((λ + 180) / (180/2k ) ┘ mod 2(k +1) 其中└┘代表向下取整， mod 是模运算。相反的，知道一个第 k 层的瓦片的行列号（X， Y），同样也可以计算出瓦片的经度范围 west ～ east 和纬度范围 south ～ north ： west = ((x mod 2(k+1)) × 180/2k) – 180 east = west + 180/2k south = ((y mod 2k) × 180/2k) – 90 north = south + 180/2k 因此，只需给出待构建金字塔的层次，每层的瓦片划分是已知的。
     本发明流程可采用计算机软件技术实现自动运行。以下结合附图和实施例详细说明本发明技术方案，参见图 1，实施例流程如下：步骤 1、首先进行预处理，包括统一空间参考、无效域去除和读元数据。
     实施例将 WGS-84 坐标系作为统一的空间基准，这样可以保证多源地理空间数据的无缝浏览，也提高了数据处理的效率。具体实施时，利用现有 Erdas、 Envi 等遥感软件对遥感影像进行处理，即可将其投影坐标转换成 WGS-84 坐标系。
     无效域的去除有两种途径，即无效值过滤与多边形裁剪。无效值过滤方法比较简单高效，在需要去除成像不理想的区域时则可使用多边形裁剪。多边形生成可以采用手工勾勒方式，在手工勾勒的基础上还可以进行精细调位。实施例中，遥感影像成像效果较好，仅仅影像边缘有一些无效区域，可以直接设定无效值，采用无效值过滤的方法进行无效域去除。
     实施例读取原始遥感影像（即元数据）的信息，通过调用 Gdal 开源库读取遥感影像的原始数据文件信息实现，包括原始遥感影像的波段数 bn、原始遥感影像的尺寸 pw×ph、原始遥感影像的影像分辨率 pd 以及原始遥感影像对应的地理范围（rw, rs, re, rn）。实施例采用了现有技术中的 Gdal 库 (Geospatial Data Abstraction Library)， Gdal 库是一个在 X/MIT 许可协议下的开源栅格空间数据转换库。它利用抽象数据模型来表达所支持的各种文件格式，能够读取、写入、转换、处理各种栅格数据格式，并且是跨平台的。
     实施例定义瓦片结构体： typedef struct tagMemTile{ int count;// 定义 count 为整形变量 int status;// 定义 status 为整形变量 int level;// 定义 level 为整形变量 int row;// 定义 row 为整形变量 int col;// 定义 col 为整形变量 double pixelSpan;// 定义 pixelSpan 为浮点数 FIBITMAP* pTileBmp;// 定义 pTileBmp 为指向 FIBITMAP 的指针 }MemTile; 其中 tagMemTile 为结构体名称， MemTile 为结构体 tagMemTile 的变量名， count 为瓦片中有效值的个数， level 是瓦片所在的层数， row 是瓦片的行号， col 的瓦片的列号， pixelSpan 是瓦片中一个像素的分辨率， FIBITMAP 是 FreeImage 中用来存储瓦片数据的结构， pTileBmp 是指向这个结构的指针； status 为瓦片状态，可用于标识内存需求最小层的上层瓦片还是下层的四分瓦片以及瓦片是否已经处理完。步骤 2，待构建金字塔的最顶层为分辨率最低的层次，标记为第 0 层；根据原始遥感影像对应的地理范围（rw, rs, re, rn）和第 0 层瓦片分辨率 yd[0]，计算出第 0 层瓦片的起止范围，计算方式如下， sr[0] = (rw+180+1e-8)/yd[0] er[0] = (re+180+1e-8)/yd[0] sc[0] = (rs+90+1e-8)/yd[0] ec[0] = (rn+90+1e-8)/yd[0] 其中， sr[0] 为第 0 层的开始行、 er[0] 为第 0 层的结束行、 sc[0] 为第 0 层的开始列、 ec[0] 为第 0 层的结束列、 e 为数学参数。
     实施例中，第 0 层瓦片分辨率 yd[0] 根据原始遥感影像的影像分辨率 pd 计算，计算公式如下， yd[0]=pd×tilesize×2l -1 其中， l 为预设的待构建金字塔的层数， tilesize 为预设的瓦片大小，即一个瓦片包含的像素数目。
     步骤 3，比较第 0 层一个瓦片对应原始遥感影像所需的内存 memvolume[0] 与预设的可申请的内存，预设的可申请内存上限为 ms 个字节，当所需内存大小 memvolume[0] 不大于可申请内存上限 ms 时进入步骤 4，否则说明第 0 层中一个瓦片的相关数据不可以完全载入内存，进入步骤 5。
     此时第 0 层瓦片的瓦片结构体的所需内存大小可以忽略。因此实施例中，第0层存储一个瓦片对应原始遥感影像所需的内存大小 memvolume[0] 根据原始遥感影像的波段数 bn 计算，计算公式如下， memvolume[0] =pixelvolume[0]×bn 其中， pixelvolume[0] 为第 0 层中存储一个瓦片覆盖原始遥感影像的其中一个波段的信息所需的内存大小，根据第 0 层瓦片分辨率 yd[0]、原始遥感影像的尺寸 pw×ph 和原始遥
     感影像对应的地理范围（rw, rs, re, rn）计算，计算公式如下， pixelvolume[0] = (yd[0]×pw/(re-rw)+0.5)×( yd[0]×ph/(rn-rs)+0.5) ×sizeof(type) 其中 type 为存储像素值的数据类型， sizeof（）是获取数据类型在计算机内占有字节数的函数；具体实施时，内存上限 ms 一般可设到 500M，即 500×1024×1024 个字节。
     步骤 4，逐个处理第 0 层的每一个瓦片，处理方式为计算第 0 层中某瓦片在原始遥感影像的像素覆盖范围，并将该范围对应的原始影像数据的数据块读入到内存，记为数据块 i，数据块 i 对应的第 0 层瓦片及其四分后代瓦片的瓦片结构体记为集合 sti，将集合 sti 初始化；遍历集合 sti，当遍历到集合 sti 中的某个瓦片结构体时，在数据块 i 中获取相应的数据，通过插值获得瓦片数据并编码；处理第 0 层的所有瓦片后进入步骤 6。
     为便于实施参考起见，提供实施例计算第 0 层中某瓦片在原始遥感影像的像素覆盖范围具体实现方法，设瓦片对应的地理范围为（m_RgnWest, m_RgnEast, m_RgnSouth, m_RgnNorth），该瓦片对应原始影像的像素范围为（xw,xe,ys,yn），步骤 1 中读取遥感影像的元数据中遥感影像的左下角地面坐标为（x0,y0） , 遥感影像的右上角地面坐标为（x1,y1），遥感影像像素在行方形上的地面分辨率为 Dr，遥感影像像素在列方形上的地面分辨率为 Dc，计算公式如下： xw=( m_RgnWest –x0 )/ Dc + 0.5f ys=(y1 - m_RgnSouth)/ Dr + 0.5f xe=( m_RgnEast –x0 )/ Dc + 0.5f yn=(y1 - m_RgnNorth)/ Dr + 0.5f 其中， f 表示 0.5 是 float 类型。
     第 0 层为分辨率最低层，第 0 层一个瓦片对应遥感影像的数据量大于其他各层。因此如果第 0 层的一个瓦片可以读入内存，则只需逐个计算第 0 层的每一个瓦片及其所有四分瓦片，得到各瓦片中所有像素值，所得结果写入数据集中。
     步骤 5，在第 0 层一个瓦片对应遥感影像不能完全载入内存时，实施例计算出所需内存容量最小的那一层，记为 wl。即假设某一层中存储一个瓦片对应的原始遥感影像所需的内存大小为 memvolume[level]，存储从第 0 层到该层上一层的所有瓦片结构体所需的内存大小是 tlstuvolume[level]，设 memvolume[level]+tlstuvolume[level] 取值最小的层为内存需求最小层 wl。倘若实施时，发现 memvolume[level]+tlstuvolume[level] 仍大于内存上限 ms，可以设计此时跳出程序，提示用户修改申请的缓存值。基于这种内存需求最小层设定，执行以下步骤：步骤 5.1，求取内存需求最小层 wl。
     从第 wl 层开始计算多分辨率金字塔，算法对于内存的需求是最小的。这样可以根据设置的缓存大小确定从哪一层开始构建，实施例中 wl=max(0, └ (l -e -1)/2 ┘ ) 其中 e = └ log4(m /3) ┘， └┘代表向下取整； m 表示第 0 层瓦片个数，根据步骤 2 所得第 0 层的开始行 sr[0]、第 0 层的结束行 er[0]、第 0 层的开始列 sc[0] 和第 0 层的结束列 ec[0] 计算，计算公式如下，m=(er[0]-sr[0]+1)×(ec[0]-sc[0]+1) 为便于理解本发明技术方案起见，下面给出了实施例中求取内存需求最小层层数方式的具体说明：设原始遥感影像有 b 个波段，输出的金字塔有 l 层，自上而下，第0层（分辨率最低层） (l -1) 有m （m ≥ 1）个瓦片，第 l -1 层（分辨率最高层）则有 m ×4 个瓦片，总共有 (m/3)×(4l -1) 个瓦片。每个瓦片像素大小是 s 。
     如果从第 (l -1) 层开始，即缓存第 (l -1) 层以上所有瓦片结构体，计为 (l -1) (m /3)×(4 -1)，则所需内存为： Memvolume-1[l -1]+tlstuvolume[l -1]=b ×s 2×((m /3)×(4(l -1)-1)) + b ×s 2×40 = b ×s 2((m /3)(4(l -1)-1) + 40) 如果从第 (l -2) 层开始，即缓存从第 (l -2) 层以上的所有瓦片，则所需内存为： 2 (l -2) memvolume[l -2]+tlstuvolume[l -2]=b ×s ((m /3)(4 -1) + 41) 如果从第 k （0 ≤ k ≤ l -1）层开始，即缓存从第 k 层以上的所有瓦片，则所需内存为： 2 k (l -k -1) T k = memvolume[k ]+tlstuvolume[k ]=b ×s ×((m /3)(4 -1) + 4 )。 2 k (l -k -2)
     T k+1 -T k = 3×b ×s ×((m /3)×4 -4 )，可见T k+1 和 Tk 大小由∆ k = n ×4k -4(l -k -2) 来决定，其中， n =m /3。设 4e ≤ n <4(e +1)，由于 m ≥ 1， e ≥ -1，分两种情况来讨论： (e +k ) (l -k -2) (l +k -2 ) (l -k -2) 1） (l -e ) ≤ 2，此时， ≥4 -4 =0，即 ∆ k 恒为正值，因此 k =0 ∆ k ≥ 4 -4 时内存需求最小。
     2） (l -e )>2，此时：设 (l -e -2)=2h ， ≥ 4(e +h )-4 (e+h)=0， ∆ h = n ×4h -4(l -h -2) ∆ h-1 = n ×4(h -1)-4(l -h -1) < 4(e +h )-4(e +h +1)<0，即 T (k+1 ) ≥ T k     设 (l -e -2)=2h +1， ∆ h +1 = n ×4(h +1)-4(l -h -3) ≥ 4(e +h +1)-4 (e+h)>0， ∆ h = n ×4h -4(l -h -2) < 4(e +h +1)-4(e +h +1)=0，即 T (k+2 ) >T (k+1 )< T k，因而 k =(h +1) 时内存需求最小。
     综上所述， k =max(0, └ (l -e -1)/2 ┘ ) 时，算法对内存的需求最小，其中 e = └ log4(m /3) ┘。
     当 m =1，即第 0 层仅有 1 个瓦片，此时 k = max(0, └ l /2 ┘ ) = └ l /2 ┘，即在等于└ l /2 ┘时内存需求最小。
     计算出内存需求最小层 wl= max(0, └ (l -e -1)/2 ┘ ) 后，就可以按对第 0 层到第 wl-1 层、第 wl 层到第 l -1 层分别按两种不同方式计算各层瓦片的瓦片数据。对第 0 层到第 wl-1 层，是将瓦片结构体初始化后加载到内存中，每从原始数据文件读出一个数据块，遍历这些瓦片：如果某瓦片中有像素所需数据位于该数据块内，插值计算出该像素的值。而对于第 wl 层及该层以下的瓦片则和步骤 4 处理方式一致。实施例采用步骤 5.2 和步骤 5.3 实现。
     步骤 5.2，将第 0 层到第 wl-1 层的所有瓦片结构体记为集合 ts，将集合 ts 初始化，并加载到内存之中。
     例如第 0 层到第 wl-1 层某瓦片 pMemTile 的层行列为（k， rr1， cc1），瓦片分辨率为 Dc，瓦片大小为 tilesize，则初始化该瓦片为： pMemTile=new MemTile;// 创建一个新的瓦片结构体 pMemTile->count=0;// 将 count 初始化为 0 pMemTile->level=k;// 对层号进行赋值 pMemTile->row=rr1;// 对行号进行赋值 pMemTile->col=cc1;// 对列号进行赋值 pMemTile->pixelSpan=Dc/tilesize;// 计算瓦片的像素分辨率并赋值给 pixelSpan pMemTile->pTileBmp=FreeImage_ Allocate(tilesize, tilesize, 32);// 为 pTileBmp 分配内存 pMemTile->status=1;// 标识该瓦片为上层瓦片步骤 5.3，逐个处理第 w1 层的每一个瓦片，处理方式为计算第 w1 层中某瓦片在原始遥感影像的像素覆盖范围，并将该范围对应的原始影像数据的数据块读入内存，记为数据块 i，数据块 i 对应的第 w1 层瓦片及其四分后代瓦片的瓦片结构体记为集合 sti，将集合 sti 初始化；遍历集合 sti 和集合 ts，当遍历到集合 sti 中的瓦片结构体时，在数据块 i 中获取相应的数据，通过插值获得瓦片数据并编码；当遍历到集合 ts 中的瓦片结构体时，判断瓦片的覆盖范围是否包括数据块 i，如果包括则插值获得该瓦片内相应部分区域的瓦片数据，在获得该瓦片内所有瓦片数据后进行编码。当然，不包括则不进行任何处理。
     例如第 w1 层到第 l -1 层某瓦片 pMemTile 的层行列为（k， rr1， cc1），瓦片分辨率为 Dc，瓦片大小为 tilesize，则初始化该瓦片为： pMemTile=new MemTile;// 创建一个新的瓦片结构体 pMemTile->count=0;// 将 count 初始化为 0 pMemTile->level=k;// 对层号进行赋值 pMemTile->row=rr1;// 对行号进行赋值 pMemTile->col=cc1;// 对列号进行赋值 pMemTile->pixelSpan=Dc/tilesize;// 计算瓦片的像素分辨率并赋值给 pixelSpan pMemTile->pTileBmp=NULL;// 暂不为 pTileBmp 分配内存 pMemTile->status=-1;// 标识该瓦片为 w1 层或下层瓦片计算第 w1 层中某瓦片在原始遥感影像的像素覆盖范围，与步骤 4 计算第 0 层中某瓦片在原始遥感影像的像素覆盖范围具体实现方法类似：设第 w1 层某瓦片对应的地理范围为（m_RgnWest, m_RgnEast, m_RgnSouth, m_ RgnNorth），该瓦片对应原始影像的像素范围为（xw,xe,ys,yn），步骤 1 中读取遥感影像的元数据中遥感影像的左下角地面坐标为（x0,y0） , 遥感影像的右上角地面坐标为（x1,y1），遥感影像像素在行方形上的地面分辨率为 Dr，遥感影像像素在列方形上的地面分辨率为 Dc，计算公式如下： xw=( m_RgnWest –x0 )/ Dc + 0.5f ys=(y1 - m_RgnSouth)/ Dr + 0.5f xe=( m_RgnEast –x0 )/ Dc + 0.5f yn=(y1 - m_RgnNorth)/ Dr + 0.5f 其中， f 表示 0.5 是 float 类型。
     当遍历完所有第 w1 层的瓦片后，即将原始影像扫描一遍， 0 到 w1-1 层的瓦片数据也就计算完毕（如果瓦片覆盖的范围超出了原始影像的范围，则以无效值处理）。对于集合 sti 中的瓦片结构体，通过插值获得瓦片数据后可以立即编码；对于集合 ts 中的瓦片结构体，因为每个数据块可能只对应瓦片内部分区域，可以将第 w1 层的所有瓦片对应的影像数据扫描一遍，即对原始文件扫描一遍后，再进行编码。为了便于提高瓦片的处理，具体实施时可以设定当 0 到 w1-1 层的某瓦片处理完后，将 status 的值设为 0，之后便不用计算这个瓦片是否是在目前内存中数据块的范围内了。为了提高处理效率，实施例在步骤 4 和步骤 5 中，设置瓦片处理线程和瓦片计算线程两类线程，并开辟两块内存区间主缓存和备用缓存来存储数据块；瓦片处理线程将原始遥感影像的数据块读入备用缓存，并在瓦片计算线程遍历完瓦片结构体后，将备用缓存设置为主缓存，并继续将下一个数据块读入备用缓存；瓦片计算线程读取主缓存内的数据块，遍历瓦片结构体，插值计算瓦片数据。
     根据硬件的条件不同，通过调整瓦片计算线程的数目，使程序读数据和写数据的速度相匹配，实现了读写操作的同步性，节省了读写操作相互等待的时间。
     步骤 4 和 5 中，具体插值获得瓦片数据实现可采用现有技术，为便于实施参考起见，提供实施例的插值处理方法说明如下：首先根据瓦片地理范围，计算其像素范围，将该范围的像素通过 Gdal 开源库读入内存。创建的金字塔瓦片里每个点不可能和全球基础数据文件中的每个像素一一对应，所以必须对其进行重采样。遥感影像采用最近邻采样算法，直接取与的像素值为所需点的像素值作为采样值，即为最邻近点，其影像坐标值为（，）点位置最近像元地形影像采用双线性内插法，以确保数据的精度。如图 2 所示，已知、、、四点的坐标，标。
     是点对于邻近四点的水平和垂直距离，假设点与点间的距离为 1，求的影响求出点：点坐首先进行第一次线性插值，根据、对然后进行第二次线性插值，根据、对的影响求出点：最后对、两点做内插，求出：本实施例是遥感影像，故采用最近邻采样算法从数据文件中计算出瓦片数据。如果在输出影像瓦片中，所有单位均是有效值，则输出 JPEG 格式瓦片，否则输出 PNG 格式瓦片。
     步骤 4 和 5 中，具体对瓦片数据的编码实现可采用现有技术，为便于实施参考起见，提供实施例的编码方法说明如下：瓦片编码功能的目的在于将瓦片的多维标识（层、行和列）映射为一维标识，从而提高瓦片的定位速度。除了位置与尺度两个基本属性之外，对于卫星数据来说，还有时相属性。为此从可扩展角度出发，将时相信息集成到瓦片的编码之中，即在空间编码的基础上进行扩展，加入时相信息，从而形成瓦片的时空一体化编码，即在空间编码的基础上进行扩展，加入时相信息，从而形成瓦片的时空一体化编码，即多时相瓦片编码 = 时间编码 + 空间编码，该编码在一个金字塔内是惟一的，可通过时空索引进行查找和存储。另外，经过扩展后的 Morton 编码在商用数据库中可用标准数据类型进行存储和管理 , 因此是比较理想的编码方案。这样将每一个空间数据块作为数据库中的一条记录 , 即可在商用数据库中进行存储 , 在读取每一块空间数据时 , 只要根据制定的时相和地面范围信息即可获取相应范围的影像块编码 , 通过编码即可获得相应的数据块的数据。本实施例中的瓦片编码用一个 64 比特整型表示，从高位到低位共 4 部分，其内部结构如图 3 所示：包括有效位（共 1 位）、时间编码（共 10 位）、层编码（共 5 位）、行列号 Morton 编码（共 48 位）。
     步骤 6、多分辨率金字塔的输出。
     实施例在完成多分辨率金字塔构建后，输出多分辨率金字塔文件。
     总结本发明技术方案过程，如图 4 所示：首先将影像的坐标系统一到 WGS84 坐标系，然后去除无效域，调用 Gdal 库读取元数据。然后对每个瓦片进行处理：通过计算瓦片地理范围获取对应的像素范围，在该像素范围的数据块中通过插值计算瓦片数据对每一个瓦片进行编码，最后输出多分辨率金字塔文件。
     为说明本发明效果起见，对一组大小不一的遥感影像（1G、 5G、 10G、 15G、 20G、 25G、 30G、 35G）分别采用现有一般构建多分辨率金字塔方法和本发明所提供构建多分辨率金字塔方法，对每个遥感影像都构建了 9 层金字塔。实施本发明所提供构建多分辨率金字塔方法的程序采用了 8 个计算线程，缓存大小设置为 1024M，由下表可看出该算法大大提高了构建金字塔的效率。FileSize(inGB) 目前构建多分辨率金字塔方法本发明构建多分辨率金字塔方法 1 5 8m31s 50m38s 55s 4m23s 10 15 1h51m 2h44m 9m21s 13m54s 20 3h39m 21m12s 25 4h35m 26m42s 30 35 4h58s 5h48m 33m5s 35m44s其中， h 表示小时， m 表示分钟， s 表示秒。
     本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

资源描述

《一种面向虚拟地球的栅格数据快速建库方法.pdf》由会员分享，可在线阅读，更多相关《一种面向虚拟地球的栅格数据快速建库方法.pdf（17页珍藏版）》请在专利查询网上搜索。

1、(10)申请公布号 CN 102509022 A (43)申请公布日 2012.06.20 C N 1 0 2 5 0 9 0 2 2 A *CN102509022A* (21)申请号 201110366702.8 (22)申请日 2011.11.18 G06F 19/00(2011.01) (71)申请人武汉大学地址 430072 湖北省武汉市武昌区珞珈山武汉大学 (72)发明人向隆刚曾珍 (74)专利代理机构武汉科皓知识产权代理事务所(特殊普通合伙) 42222 代理人薛玲 (54) 发明名称一种面向虚拟地球的栅格数据快速建库方法 (57) 摘要本发明涉及一种面向虚拟地球的栅格。

2、数据快速建库方法。该方法提出，如果第0层一个瓦片对应的原始遥感影像数据无法一次完全载入内存，从中间层（内存需求最小层）开始构建金字塔。这样打破了以往逐层构建的局限性，仅需扫描原始数据文件一次，减少了I/O操作次数，消除了处理海量数据时读取数据的瓶颈问题。本发明还设计了瓦片处理线程和瓦片计算线程两类线程，让两类线程同时运行，充分利用了系统资源，根据硬件不同的条件，配置多个瓦片计算线程，使程序读数据和写数据的速度相匹配，实现了读写操作的同步性，节省了读写操作相互等待的时间。提高了金字塔构建的效率，尤其适合于对超大遥感数据文件的处理。 (51)Int.Cl. 权利要求书2页。

3、说明书11页附图3页 (19)中华人民共和国国家知识产权局 (12)发明专利申请权利要求书 2 页说明书 11 页附图 3 页 1/2页 2 1.一种面向虚拟地球的栅格数据快速建库方法，其特征在于，包括以下步骤：步骤1，读取原始遥感影像的信息，包括原始遥感影像的波段数bn、原始遥感影像的尺寸pwph、原始遥感影像的影像分辨率pd以及原始遥感影像对应的地理范围（rw, rs, re, rn）；定义瓦片结构体，所述瓦片结构体是用于在内存中存储瓦片数据并对瓦片进行标识的数据结构体；步骤2，待构建金字塔的最顶层为分辨率最低的层次，标记为第0层；根据原始遥感影像对应的地理范围（rw,。

4、 rs, re, rn）和第0层瓦片分辨率yd0，计算出第0层瓦片的起止范围，计算方式如下， sr0 = (rw+180+1e -8 )/yd0 er0 = (re+180+1e -8 )/yd0 sc0 = (rs+90+1e -8 )/yd0 ec0 = (rn+90+1e -8 )/yd0 其中，sr0为第0层的开始行、er0为第0层的结束行、sc0为第0层的开始列、 ec0为第0层的结束列、e为数学参数；第0层瓦片分辨率yd0根据原始遥感影像的影像分辨率pd计算，计算公式如下， yd0=pdtilesize2 l-1 其中，l为预设的待构建金字塔的层数，tilesize为预设的瓦片。

5、大小；步骤3，比较第0层一个瓦片对应原始遥感影像所需的内存memvolume0与预设的可申请的内存，设预设的可申请内存上限为ms个字节，当所需内存大小memvolume0不大于可申请内存上限ms时进入步骤4，否则进入步骤5；所述第0层一个瓦片对应原始遥感影像所需的内存大小memvolume0根据原始遥感影像的波段数bn计算，计算公式如下， memvolume0 =pixelvolume0bn 其中，pixelvolume0为第0层中存储一个瓦片覆盖原始遥感影像的其中一个波段的信息所需的内存大小，根据第0层瓦片分辨率yd0、原始遥感影像的尺寸pwph和原始遥感影像对应的地理范围（r。

6、w, rs, re, rn）计算，计算公式如下， pixelvolume0 = (yd0pw/(re-rw)+0.5)( yd0ph/(rn-rs)+0.5) sizeof(type) 其中type为存储像素值的数据类型，sizeof（）是获取数据类型在计算机内占有字节数的函数；步骤4，逐个处理第0层的每一个瓦片，处理方式为计算第0层中某瓦片在原始遥感影像的像素覆盖范围，并将该范围对应的原始影像数据的数据块读入到内存，记为数据块i，数据块i对应的第0层瓦片及其四分后代瓦片的瓦片结构体记为集合st i ，将集合st i 初始化；遍历集合st i ，当遍历到集合st i 中的某个瓦片结。

7、构体时，在数据块i中获取相应的数据，通过插值获得瓦片数据并编码；处理第0层的所有瓦片后进入步骤6；步骤5，设某一层中存储一个瓦片对应的原始影像所需的内存大小为 memvolumelevel，存储从第0层到该层上一层的所有瓦片结构体所需的内存大小是权利要求书CN 102509022 A 2/2页 3 tlstuvolumelevel，设memvolumelevel+tlstuvolumelevel取值最小的层为内存需求最小层wl；进行以下步骤，步骤5.1，求取内存需求最小层wl，计算公式如下， wl=max(0, (l-e-1)/2) 其中e = log 4 (m/3) ，代。

8、表向下取整；m表示第0层瓦片个数，根据步骤2所得第0层的开始行sr0、第0层的结束行er0、第0层的开始列sc0和第0层的结束列ec0计算，计算公式如下， m=(er0-sr0+1)*(ec0-sc0+1) 步骤5.2，将第0层到第wl-1层的所有瓦片结构体记为集合ts，将集合ts初始化，并加载到内存之中；步骤5.3，逐个处理第w1层的每一个瓦片，处理方式为计算第w1层中某瓦片在原始遥感影像的像素覆盖范围，并将该范围对应的原始影像数据的数据块读入内存，记为数据块 i，数据块i对应的第w1层瓦片及其四分后代瓦片的瓦片结构体记为集合st i ，将集合st i 初始化；遍历集合st i 。

9、和集合ts，当遍历到集合st i 中的某个瓦片结构体时，在数据块i中获取相应的数据，通过插值获得瓦片数据并编码；当遍历到集合ts中的某个瓦片结构体时，判断瓦片的覆盖范围是否包括数据块i，如果包括则插值获得该瓦片内相应部分区域的瓦片数据，在获得该瓦片内所有瓦片数据后进行编码；处理第w1层的所有瓦片后进入步骤6；步骤6，输出多分辨率金字塔。 2.如权利要求1所述面向虚拟地球的栅格数据快速建库方法，其特征在于：步骤4和步骤5中，设置瓦片处理线程和瓦片计算线程两类线程，并开辟两块内存区间主缓存和备用缓存来存储数据块；瓦片处理线程将原始遥感影像的数据块读入备用缓存，并在瓦片计算线程遍。

10、历完瓦片结构体后，将备用缓存设置为主缓存，并继续将下一个数据块读入备用缓存；瓦片计算线程读取主缓存内的数据块，插值计算瓦片数据。 3.如权利要求1或2所述面向虚拟地球的栅格数据快速建库方法，其特征在于：瓦片结构体包括变量count、status、level、row、col、pixelSpan和指针pTileBmp，其中变量 count标识瓦片中有效值的个数、变量status标识瓦片状态、变量level标识瓦片所在的层数、变量row标识瓦片的行号、变量col标识瓦片的列号、变量pixelSpan标识瓦片中每个像素的分辨率，指针pTileBmp指向瓦片数据所在内存位置。权利要求。

11、书CN 102509022 A 1/11页 4 一种面向虚拟地球的栅格数据快速建库方法技术领域 0001 本发明涉及空间数据组织技术领域，尤其是涉及一种面向虚拟地球的栅格数据快速建库方法。背景技术 0002 如今越来越多的空间数据以栅格的形式存储（包括地形数据和影像数据）。栅格数据成为了重要的数据源。随着对地观测技术的快速发展，人类能够快速地获取多时相、多波段、多分辨率的海量遥感影像数据。随着分辨率的不断提高，全球数据的数据量也在以惊人的速度增加，如仅福建省的30m30m分辨率的24位bmp格式的遥感影像就有700MB，若为 1m1m分辨率，数据量将达到630GB。如何去存储、组织。

12、和管理栅格数据，如何高效地检索、显示用户感兴趣的数据，使栅格数据能够更好地为实际应用服务越来越受人们的关注。 0003 有关文献：高峰，安培浚.国际空间和对地观测技术发展战略新动向J .遥感技术与应用，2008， (06)：686-696.;张更新，张昭，朱江.全球对地观测卫星现状及其产业发展综述J.数字通信世界， 2009，(10) ：72-75.；王恩泉.中国版Google Earth的空间数据组织与管理研究D.北京：中国测绘科学研究院，2007.。 0004 随着Google Earth的发布与流行，虚拟地球技术进入了人们的视野。虚拟地球能够模拟显示大范围高分辨率的海量遥感影像。

13、数据，不仅用于各个行业领域内的专业应用研究, 也能够满足人们足不出户即可感知整个地球。目前比较热门虚拟地球产品除了Google Earth，还有美国航天航空局（NASA）的WorldWind、Bing 3D和Skyline等，国内不少公司与高校开展了相关研究，建立自主版权的虚拟数字地球平台，如武汉大学测绘遥感信息工程国家重点实验室自主研制的虚拟地球系统GeoGlobe。 0005 有关文献：苗放，叶成名，刘瑞，孔祥生，李康荣，徐松浦，张远红.新一代数字地球平台与“数字中国”技术体系架构探讨J.测绘科学，2007，(06)：157-158.；http:/ 0006 虚拟地球通常采用多。

14、分辨率瓦片金字塔结构来组织栅格数据，以提供快速的网络调度和三维渲染。对栅格数据建库的过程也就是构建多分辨率金字塔的过程。通过构建多分辨率金字塔，在多尺度金字塔结构索引机制的配合下，可以快速获取指定分辨率和范围的栅格数据，并适合于网络传输和渲染处理，从而保证用户根据视景内的范围，在任一角度浏览的数据只是金字塔的子集。面向多分辨率瓦片金字塔的数据组织与调度是虚拟地球系统得以实现网络三维高效可视化的核心思想。如何对栅格数据快速构建多分辨率的金字塔一直是人们研究的热点。 0007 有关文献：龚健雅，李斌等.当代GIS的若干理论与技术.武汉：武汉测绘科技大学出版社.；吴晨琛.GeoGlob。

15、e中尺度空间数据集管理机制D.武汉：武汉大学，2008.。 0008 早期对金字塔的研究重点主要是如何提高成像质量和传输。1979年，Tanimoto等人提出了均值金字塔。均值金字塔的思想简单而直观，但应用到图像的渐进传输中却能获得很好的效果，引起了人们的关注。说明书CN 102509022 A 2/11页 5 0009 有关文献：朱程辉，何勇，王金玲.基于小波金字塔的快速图像匹配算法J .微计算机信息， 2010，(10)：127-128.。 0010 同年，K.R.Sloan和Tanimoto又提出了Reduced-sum金字塔的结构，该结构克服了均值金字塔数据量过大的问题。

16、，将金字塔中的像素数减少到与原始图像相同。 0011 有关文献：刘国权，李守轩.基于小波图像金字塔的SSDA快速模板匹配算法J，科技广场， 2007，(11)：134-136.。 0012 但无论是均值金字塔还是Reduced-sum金字塔，其画面质量都不够理想，图像的混叠失真较大。于是在1983年，P.J.Burt等与J.Nao等提出了高斯-拉普拉斯金字塔结构，该结构在采样的同时利用高斯型低通滤波器对图像进行平滑，从而获得良好的视觉效果。 0013 有关文献：丁素英.基于小波变换的Laplacian金字塔图像数据压缩J .潍坊学院学报， 2009，(04)：34-36.；叶勤，张小。

17、虎，王栋.一种基于区域生长与空间形状约束的高分辨率遥感图像道路提取方法J.遥感信息，2010，(02)：25-29.。 0014 而后的学者对如何对大数据量的金字塔进行有效的存储和管理进行了研究。1989 年，又有学者提出了Difference金字塔结构，通过将相邻层的像素相减来去除相关性，从而减少金字塔的信息熵。 0015 有关文献：欧阳平，张玉方.形态学开闭运算在居民地边缘检测中的应用J.测绘通报，2009，(01) ：40-41.。 0016 随着计算机技术的快速发展，计算机能提供足够的空间存储金字塔数据，但是随着分辨率的提升，即使覆盖较小区域，一个遥感数据文件的尺寸也是比较大的。

18、，可能达到上十G，甚至上百G，从这样的超大数据文件（1个或者多个）构建出多分辨率瓦片金字塔，其瓶颈在于读操作。 0017 目前，一般采用的金字塔构建方法是逐层逐块构建金字塔，对每层金字塔都获取逐个像素的扫描行，然后将扫描行放人缓存，当缓存已满时采用邻近法替换数据。其缺点是多次扫描原始遥感文件，计算一个顶层瓦片时，涉及到原始遥感文件的较大范围，I/O效率低下。发明内容 0018 针对目前在虚拟地球系统中对栅格数据构建多分辨率金字塔的方法无法高效处理大数据量的遥感影像，为了提高构建多分辨率金字塔效率，本发明提供的技术方案是面向虚拟地球的一种针对海量栅格数据构建金字塔的方法，主要是针。

19、对tif、img等可以存储大数据量格式的文件。 0019 本发明的技术方案为一种面向虚拟地球的栅格数据快速建库方法，包括以下步骤：步骤1，读取原始遥感影像的信息，包括原始遥感影像的波段数bn、原始遥感影像的尺寸pwph、原始遥感影像的影像分辨率pd以及原始遥感影像对应的地理范围（rw, rs, re, rn）；定义瓦片结构体，所述瓦片结构体是用于在内存中存储瓦片数据并对瓦片进行标识的数据结构体，步骤2，待构建金字塔的最顶层为分辨率最低的层次，标记为第0层；根据原始遥感影说明书CN 102509022 A 3/11页 6 像对应的地理范围（rw, rs, re, rn）和第0层。

20、瓦片分辨率yd0，计算出第0层瓦片的起止范围，计算方式如下， sr0 = (rw+180+1e -8 )/yd0 er0 = (re+180+1e -8 )/yd0 sc0 = (rs+90+1e -8 )/yd0 ec0 = (rn+90+1e -8 )/yd0 其中，sr0为第0层的开始行、er0为第0层的结束行、sc0为第0层的开始列、 ec0为第0层的结束列、e为数学参数；第0层瓦片分辨率yd0根据原始遥感影像的影像分辨率pd计算，计算公式如下， yd0=pdtilesize2 l-1 其中，l为预设的待构建金字塔的层数，tilesize为预设的瓦片大小；步骤3，比较第0层一个瓦。

21、片对应原始遥感影像所需的内存memvolume0与预设的可申请的内存，预设的可申请内存上限为ms个字节，当所需内存大小memvolume0不大于可申请内存上限ms时进入步骤4，否则进入步骤5；所述第0层一个瓦片对应原始遥感影像所需的内存大小memvolume0根据原始遥感影像的波段数bn计算，计算公式如下， memvolume0 =pixelvolume0bn 其中，pixelvolume0为第0层中存储一个瓦片覆盖原始遥感影像的其中一个波段的信息所需的内存大小，根据第0层瓦片分辨率yd0、原始遥感影像的尺寸pwph和原始遥感影像对应的地理范围（rw, rs, re, rn）计算，。

22、计算公式如下， pixelvolume0 = (yd0pw/(re-rw)+0.5)( yd0ph/(rn-rs)+0.5) sizeof(type) 其中type为存储像素值的数据类型，sizeof（）是获取数据类型在计算机内占有字节数的函数；步骤4，逐个处理第0层的每一个瓦片，处理方式为计算第0层中某瓦片在原始遥感影像的像素覆盖范围，并将该范围对应的原始影像数据的数据块读入到内存，记为数据块i，数据块i对应的第0层瓦片及其四分后代瓦片的瓦片结构体记为集合st i ，将集合st i 初始化；遍历集合st i ，当遍历到集合st i 中的某个瓦片结构体时，在数据块i中获取相应的数据。

23、，通过插值获得瓦片数据并编码；处理第0层的所有瓦片后进入步骤6；步骤5，设某一层中存储一个瓦片对应的原始影像所需的内存大小为 memvolumelevel，存储从第0层到该层上一层的所有瓦片结构体所需的内存大小是 tlstuvolumelevel，设memvolumelevel+tlstuvolumelevel取值最小的层为内存需求最小层wl；进行以下步骤，步骤5.1，求取内存需求最小层wl，计算公式如下， wl=max(0, (l-e-1)/2) 其中e = log 4 (m/3) ，代表向下取整；m表示第0层瓦片个数，根据步骤2所得第0层的开始行sr0、第0层的结束行er0、第0。

24、层的开始列sc0和第0层的结束列ec0计算，计算公式如下，说明书CN 102509022 A 4/11页 7 m=(er0-sr0+1)(ec0-sc0+1) 步骤5.2，将第0层到第wl-1层的所有瓦片结构体记为集合ts，将集合ts初始化，并加载到内存之中；步骤5.3，逐个处理第w1层的每一个瓦片，处理方式为计算第w1层中某瓦片在原始遥感影像的像素覆盖范围，并将该范围对应的原始影像数据的数据块读入内存，记为数据块 i，数据块i对应的第w1层瓦片及其四分后代瓦片的瓦片结构体记为集合st i ，将集合st i 初始化；遍历集合st i 和集合ts，当遍历到集合st i 中的某个。

25、瓦片结构体时，在数据块i中获取相应的数据，通过插值获得瓦片数据并编码；当遍历到集合ts中的某个瓦片结构体时，判断瓦片的覆盖范围是否包括数据块i，如果包括则插值获得该瓦片内相应部分区域的瓦片数据，在获得该瓦片内所有瓦片数据后进行编码；处理第w1层的所有瓦片后进入步骤6；步骤6，输出多分辨率金字塔。 0020 而且，步骤4和步骤5中，设置瓦片处理线程和瓦片计算线程两类线程，并开辟两块内存区间主缓存和备用缓存来存储数据块；瓦片处理线程将原始遥感影像的数据块读入备用缓存，并在瓦片计算线程遍历完瓦片结构体后，将备用缓存设置为主缓存，并继续将下一个数据块读入备用缓存；瓦片计算线程读取主。

26、缓存内的数据块，遍历瓦片结构体，插值计算瓦片数据。 0021 而且，瓦片结构体包括变量count、status、level、row、col、pixelSpan和指针 pTileBmp，其中变量count标识瓦片中有效值的个数、变量status标识瓦片状态、变量 level标识瓦片所在的层数、变量row标识瓦片的行号、变量col标识瓦片的列号、变量 pixelSpan标识瓦片中每个像素的分辨率，指针pTileBmp指向瓦片数据所在内存位置。 0022 本发明的优势在于从中间层（算法对于内存的需求是最小的层）开始构建金字塔，打破了以往逐层构建的局限性，仅需扫描原始数据文件一次，减少了I/O操作次。

27、数，消除了处理海量数据时读取数据的瓶颈问题，这对于超大文件来说提高了金字塔构建的效率。本发明还充分利用了现代计算机的多处理器多核，可并行计算的特点；采用瓦片处理线程和瓦片计算线程两类线程同时运行，充分利用了系统资源，通过配置多个瓦片计算线程，使程序读数据和写数据的速度相匹配，实现了读写操作的同步性，节省了读写操作相互等待的时间。附图说明 0023 图1为本发明实施例的流程图；图2为本发明实施例的双线性内插示意图；图3为本发明实施例的瓦片时空统一编码示意图。 0024 图4 为本发明实施例的数据处理流程图。 0025 具体实施方式为便于理解本发明技术方案起见，首先说明全球多分辨。

28、率金字塔：全球多分辨率金字塔是基于球面的一种可以无限细分,但又不改变形状的地球体拟合格网,当细分到一定程说明书CN 102509022 A 5/11页 8 度时,可以达到模拟地球表面的目的。本发明中采用的是通过等经纬度格网划分的全球多分辨率金字塔模型。在栅格数据金字塔中，按照全球等经纬度格网划分，每一层被裁切成统一形状、大小的数据块，这样的一个数据块称之为瓦片。 0026 这种全球多分辨率金字塔非常适合用来组织和管理多分辨率的海量卫星影像数据，因为它遵循以下规律： 1）全球地理坐标的经度范围为-180 o , +180 o ，纬度范围为-90 o ，+90 o ； 2）第k层金字。

29、塔的分辨率为第(k+l)层金字塔的2倍； 3）在任意层内，横向和纵向金字塔瓦片的数目之比为2:1； 4）金字塔瓦片的编码顺序由左到右，由下到上；根据上述规则，很容易计算出每层瓦片的个数，同样也能快速计算出给定经纬度的点落在哪一个瓦片内。例如，已知一个点的经纬度为，该点在第k层金字塔的行号RowNo和列号ColNo可以通过下面的公式计算出： RowNo= (+ 90) / (180/2 k ) mod 2 k ColNo = (+ 180) / (180/2 k ) mod 2 (k+1) 其中代表向下取整，mod是模运算。 0027 相反的，知道一个第k层的瓦片的行列号（X，Y），同样也。

30、可以计算出瓦片的经度范围westeast和纬度范围southnorth： west = (x mod 2 (k+1) ) 180/2 k ) 180 east = west + 180/2 k south = (y mod 2 k ) 180/2 k ) 90 north = south + 180/2 k 因此，只需给出待构建金字塔的层次，每层的瓦片划分是已知的。 0028 本发明流程可采用计算机软件技术实现自动运行。以下结合附图和实施例详细说明本发明技术方案，参见图1，实施例流程如下：步骤1、首先进行预处理，包括统一空间参考、无效域去除和读元数据。 0029 实施例将WGS-84坐标。

31、系作为统一的空间基准，这样可以保证多源地理空间数据的无缝浏览，也提高了数据处理的效率。具体实施时，利用现有Erdas、Envi等遥感软件对遥感影像进行处理，即可将其投影坐标转换成WGS-84坐标系。 0030 无效域的去除有两种途径，即无效值过滤与多边形裁剪。无效值过滤方法比较简单高效，在需要去除成像不理想的区域时则可使用多边形裁剪。多边形生成可以采用手工勾勒方式，在手工勾勒的基础上还可以进行精细调位。实施例中，遥感影像成像效果较好，仅仅影像边缘有一些无效区域，可以直接设定无效值，采用无效值过滤的方法进行无效域去除。 0031 实施例读取原始遥感影像（即元数据）的信息，通过调用Gd。

32、al开源库读取遥感影像的原始数据文件信息实现，包括原始遥感影像的波段数bn、原始遥感影像的尺寸pwph、原始遥感影像的影像分辨率pd以及原始遥感影像对应的地理范围（rw, rs, re, rn）。实施例采用了现有技术中的Gdal库 (Geospatial Data Abstraction Library)，Gdal库是一个在X/MIT许可协议下的开源栅格空间数据转换库。它利用抽象数据模型来表达所支持的各种文件格式，能够读取、写入、转换、处理各种栅格数据格式，并且是跨平台的。说明书CN 102509022 A 6/11页 9 0032 实施例定义瓦片结构体： typedef st。

33、ruct tagMemTile int count;/定义count为整形变量 int status;/定义status为整形变量 int level;/ 定义level为整形变量 int row;/ 定义row为整形变量 int col;/ 定义col为整形变量 double pixelSpan;/ 定义pixelSpan为浮点数 FIBITMAP* pTileBmp;/ 定义pTileBmp为指向FIBITMAP的指针 MemTile; 其中tagMemTile 为结构体名称，MemTile 为结构体tagMemTile的变量名，count为瓦片中有效值的个数， level是瓦片所在的层。

34、数，row是瓦片的行号，col的瓦片的列号， pixelSpan是瓦片中一个像素的分辨率，FIBITMAP是FreeImage中用来存储瓦片数据的结构，pTileBmp是指向这个结构的指针；status为瓦片状态，可用于标识内存需求最小层的上层瓦片还是下层的四分瓦片以及瓦片是否已经处理完。 0033 步骤2，待构建金字塔的最顶层为分辨率最低的层次，标记为第0层；根据原始遥感影像对应的地理范围（rw, rs, re, rn）和第0层瓦片分辨率yd0，计算出第0层瓦片的起止范围，计算方式如下， sr0 = (rw+180+1e -8 )/yd0 er0 = (re+180+1e -8 )/。

35、yd0 sc0 = (rs+90+1e -8 )/yd0 ec0 = (rn+90+1e -8 )/yd0 其中，sr0为第0层的开始行、er0为第0层的结束行、sc0为第0层的开始列、 ec0为第0层的结束列、e为数学参数。 0034 实施例中，第0层瓦片分辨率yd0根据原始遥感影像的影像分辨率pd计算，计算公式如下， yd0=pdtilesize2 l-1 其中，l为预设的待构建金字塔的层数，tilesize为预设的瓦片大小，即一个瓦片包含的像素数目。 0035 步骤3，比较第0层一个瓦片对应原始遥感影像所需的内存memvolume0与预设的可申请的内存，预设的可申请内存上限为ms个。

36、字节，当所需内存大小memvolume0不大于可申请内存上限ms时进入步骤4，否则说明第0层中一个瓦片的相关数据不可以完全载入内存，进入步骤5。 0036 此时第0层瓦片的瓦片结构体的所需内存大小可以忽略。因此实施例中，第0层存储一个瓦片对应原始遥感影像所需的内存大小memvolume0根据原始遥感影像的波段数bn计算，计算公式如下， memvolume0 =pixelvolume0bn 其中，pixelvolume0为第0层中存储一个瓦片覆盖原始遥感影像的其中一个波段的信息所需的内存大小，根据第0层瓦片分辨率yd0、原始遥感影像的尺寸pwph和原始遥说明书CN 1025090。

37、22 A 7/11页 10 感影像对应的地理范围（rw, rs, re, rn）计算，计算公式如下， pixelvolume0 = (yd0pw/(re-rw)+0.5)( yd0ph/(rn-rs)+0.5) sizeof(type) 其中type为存储像素值的数据类型，sizeof（）是获取数据类型在计算机内占有字节数的函数；具体实施时，内存上限ms一般可设到500M，即50010241024个字节。 0037 步骤4，逐个处理第0层的每一个瓦片，处理方式为计算第0层中某瓦片在原始遥感影像的像素覆盖范围，并将该范围对应的原始影像数据的数据块读入到内存，记为数据块i，数据块i对应的第。

38、0层瓦片及其四分后代瓦片的瓦片结构体记为集合st i ，将集合st i 初始化；遍历集合st i ，当遍历到集合st i 中的某个瓦片结构体时，在数据块i中获取相应的数据，通过插值获得瓦片数据并编码；处理第0层的所有瓦片后进入步骤6。 0038 为便于实施参考起见，提供实施例计算第0层中某瓦片在原始遥感影像的像素覆盖范围具体实现方法，设瓦片对应的地理范围为（m_RgnWest, m_RgnEast, m_RgnSouth, m_RgnNorth），该瓦片对应原始影像的像素范围为（xw,xe,ys,yn），步骤1中读取遥感影像的元数据中遥感影像的左下角地面坐标为（x0,y0）, 遥感影像。

39、的右上角地面坐标为（x1,y1），遥感影像像素在行方形上的地面分辨率为Dr，遥感影像像素在列方形上的地面分辨率为Dc，计算公式如下： xw=( m_RgnWest x0 )/ Dc + 0.5f ys=(y1 - m_RgnSouth)/ Dr + 0.5f xe=( m_RgnEast x0 )/ Dc + 0.5f yn=(y1 - m_RgnNorth)/ Dr + 0.5f 其中，f表示0.5是float类型。 0039 第0层为分辨率最低层，第0层一个瓦片对应遥感影像的数据量大于其他各层。因此如果第0层的一个瓦片可以读入内存，则只需逐个计算第0层的每一个瓦片及其所有四分瓦片，。

40、得到各瓦片中所有像素值，所得结果写入数据集中。 0040 步骤5，在第0层一个瓦片对应遥感影像不能完全载入内存时，实施例计算出所需内存容量最小的那一层，记为wl。即假设某一层中存储一个瓦片对应的原始遥感影像所需的内存大小为memvolumelevel，存储从第0层到该层上一层的所有瓦片结构体所需的内存大小是tlstuvolumelevel，设memvolumelevel+tlstuvolumelevel取值最小的层为内存需求最小层wl。倘若实施时，发现memvolumelevel+tlstuvolumelevel仍大于内存上限ms，可以设计此时跳出程序，提示用户修改申请的缓存值。基于。

41、这种内存需求最小层设定，执行以下步骤：步骤5.1，求取内存需求最小层wl。 0041 从第wl层开始计算多分辨率金字塔，算法对于内存的需求是最小的。这样可以根据设置的缓存大小确定从哪一层开始构建，实施例中wl=max(0, (l-e-1)/2) 其中e = log 4 (m/3) ，代表向下取整；m表示第0层瓦片个数，根据步骤2所得第0层的开始行sr0、第0层的结束行er0、第0层的开始列sc0和第0层的结束列ec0计算，计算公式如下，说明书CN 102509022 A 10 8/11页 11 m=(er0-sr0+1)(ec0-sc0+1) 为便于理解本发明技术方案起见，下面。

42、给出了实施例中求取内存需求最小层层数方式的具体说明：设原始遥感影像有b个波段，输出的金字塔有l层，自上而下，第0层（分辨率最低层）有m（m1）个瓦片，第l-1层（分辨率最高层）则有m4 (l-1) 个瓦片，总共有(m/3)(4 l -1) 个瓦片。每个瓦片像素大小是s。 0042 如果从第(l-1)层开始，即缓存第(l-1)层以上所有瓦片结构体，计为 (m/3)(4 (l-1) -1)，则所需内存为： Memvolume-1l-1+tlstuvolumel-1=bs 2 (m/3)(4 (l-1) -1) + bs 2 4 0 = bs 2 (m/3)(4 (l-1) -1) + 4 0。

43、 ) 如果从第(l-2)层开始，即缓存从第(l-2)层以上的所有瓦片，则所需内存为： memvolumel-2+tlstuvolumel-2=bs 2 (m/3)(4 (l-2) -1) + 4 1 ) 如果从第k（0kl-1）层开始，即缓存从第k层以上的所有瓦片，则所需内存为： T k = memvolumek+tlstuvolumek=bs 2 (m/3)(4 k -1) + 4 (l-k-1) )。 0043 T k+1 -T k = 3bs 2 (m/3)4 k -4 (l-k-2) )，可见T k+1 和 T k 大小由 k = n4 k -4 (l-k-2) 来决定，其中，n=m/。

44、3。设4 e n2，此时：设(l-e-2)=2h， h = n4 h -4 (l-h-2) 4 (e+h) -4 (e+h) =0， h-1 = n4 (h-1) -4 (l-h-1) 0， h = n4 h -4 (l-h-2) T (k+1) count=0;/将count 初始化为0 pMemTile-level=k;/对层号进行赋值 pMemTile-row=rr1;/对行号进行说明书CN 102509022 A 11 9/11页 12 赋值 pMemTile-col=cc1;/对列号进行赋值 pMemTile-pixelSpan=Dc/tilesize;/ 计算瓦片的像素分辨。

45、率并赋值给pixelSpan pMemTile-pTileBmp=FreeImage_ Allocate(tilesize, tilesize, 32);/为pTileBmp分配内存 pMemTile-status=1;/ 标识该瓦片为上层瓦片步骤5.3，逐个处理第w1层的每一个瓦片，处理方式为计算第w1层中某瓦片在原始遥感影像的像素覆盖范围，并将该范围对应的原始影像数据的数据块读入内存，记为数据块 i，数据块i对应的第w1层瓦片及其四分后代瓦片的瓦片结构体记为集合st i ，将集合st i 初始化；遍历集合st i 和集合ts，当遍历到集合st i 中的瓦片结构体时，在数据块i中获取。

46、相应的数据，通过插值获得瓦片数据并编码；当遍历到集合ts中的瓦片结构体时，判断瓦片的覆盖范围是否包括数据块i，如果包括则插值获得该瓦片内相应部分区域的瓦片数据，在获得该瓦片内所有瓦片数据后进行编码。当然，不包括则不进行任何处理。 0051 例如第w1层到第l -1层某瓦片pMemTile的层行列为（k，rr1，cc1），瓦片分辨率为Dc，瓦片大小为tilesize，则初始化该瓦片为： pMemTile=new MemTile;/创建一个新的瓦片结构体 pMemTile-count=0;/将count 初始化为0 pMemTile-level=k;/对层号进行赋值 pMemTile-r。

47、ow=rr1;/对行号进行赋值 pMemTile-col=cc1;/对列号进行赋值 pMemTile-pixelSpan=Dc/tilesize;/计算瓦片的像素分辨率并赋值给pixelSpan pMemTile-pTileBmp=NULL;/暂不为pTileBmp分配内存 pMemTile-status=-1;/标识该瓦片为w1层或下层瓦片计算第w1层中某瓦片在原始遥感影像的像素覆盖范围，与步骤4计算第0层中某瓦片在原始遥感影像的像素覆盖范围具体实现方法类似：设第w1层某瓦片对应的地理范围为（m_RgnWest, m_RgnEast, m_RgnSouth, m_ RgnNort。

48、h），该瓦片对应原始影像的像素范围为（xw,xe,ys,yn），步骤1中读取遥感影像的元数据中遥感影像的左下角地面坐标为（x0,y0）, 遥感影像的右上角地面坐标为（x1,y1），遥感影像像素在行方形上的地面分辨率为Dr，遥感影像像素在列方形上的地面分辨率为Dc，计算公式如下： xw=( m_RgnWest x0 )/ Dc + 0.5f ys=(y1 - m_RgnSouth)/ Dr + 0.5f xe=( m_RgnEast x0 )/ Dc + 0.5f yn=(y1 - m_RgnNorth)/ Dr + 0.5f 其中，f表示0.5是float类型。 0052 当遍历完所有第w1层的瓦片后，即将原始影像扫描。

展开阅读全文