一种并行方式栅格影像切片方法
【技术领域】
[0001] 本发明属于地理空间信息处理技术领域，涉及一种并行方式栅格影像切片方法。 技术背景
[0002] 随着卫星传感器技术以及无人机航拍技术的快速发展，遥感影像的空间和时间分 辨率都有大幅度地提高，单幅遥感影像文件的数据量也急剧增加。当前主流GIS软件以及互 联网地图应用在WebGIS(网络地理信息系统)解决方案中都广泛采用地图切片（Tile，又称 瓦片）的发布策略，这种方法通过预先将原始影像按照一定的切分规则切割成一张张大小 相同的瓦片，以加载小数据量瓦片的方式达到在网络带宽受限的条件下实现影像的高效显 示。但是目前的商业GIS软件其地图切片的生成以及发布成本较高，操作较复杂，以行业内 知名的ArcGIS系列软件为例，要将栅格影像进行切片不仅要安装ArcGIS Desktop，而且还 要安装庞大的ArcGIS Server，并且操作十分复杂繁琐，极大增加时间和物力成本。最重要 一点是随着单幅遥感影像文件的分辨率以及数据量的急剧增加，其相应的切片数量会呈现 几何级数式的增长，传统的串行算法以及商业GIS软件是通过单机预先切好瓦片，再对外提 供，这种传统切片技术计算资源利用低下，并且没有错误恢复机制，一旦切片过程中出现故 障，整个切片工作得推倒重来，无法在原有的进度上继续进行。因此在硬件性能高速发展的 情况下，如何利用高性能计算技术，方便快捷并且高效快速地进行栅格影像切片，已经成为 WebG IS地图应用中快速可视化方面必须要解决的重要问题。
[0003] 栅格影像切片是指将指定地理范围的栅格影像，在某一级别比例尺下，基于一定 的切割准则切割成若干行和列的固定尺寸的正方形栅格图片的过程，这些规整的影像切片 又称为瓦片。每个瓦片在该级别比例尺下都有一个独立的编码，该编码由瓦片的行列号以 及比例尺的级别构成。当浏览器前端在影像浏览时，根据当前的地理范围动态地计算所需 显示的影像切片的行列号，通过行列号以及当前比例尺级别得到瓦片的编码，而后通过编 码向服务器端请求相应瓦片进行显示，从而达到快速响应目的。
[0004] 目前栅格影像并行切片方法主要有三种思路:一种是CPU(Central Processing Unit，中央处理器)+GPU(Graphic Processing Unit，图形处理单元）的方式，利用GPU的计 算能力进行并行加速，这种方法并行能力受限于GPU硬件的能力，并行程度有限，而且会提 高系统架构成本;另一种是利用分布式集群系统，将切片任务划分为多个子任务，各子任务 在多个分布式节点上同时进行，这种方法可以比较方便将原有的串行方法并行化，但是当 数据规模比较大的时候，前期的数据分布式存储以及后期的结果合并都比较耗时;再一种 是基于多线程并行技术，各线程将任务划分成多个子任务，每个子任务同时进行，这种方法 可以充分利用本机计算资源，但是由于线程间并行属于细粒度并行，各线程共享父进程的 内存空间，容易造成系统的不稳定，并且这种方法可扩展性不是很好；目前基于MPI (Message Passing Interface，消息传递接口）的多进程方式进行栅格影像的并行切片研 究较少，这种方法利用共享外存的高性能集群，可扩展性强，稳定性强，可以充分利用多机 计算资源，由于集群之间共享外存，所以数据无需提前分布式存储，可以极大地提高影像切 片的效率。
【发明内容】
[0005] 本发明的目的在于提供了一种基于新思路的栅格数据并行切片方法。本发明利用 多进程技术对任务进行划分，每个进程独立进行其相应切片工作，当某个进程出现问题时， 其余进程仍可正常完成切片工作，互不影响。本发明切出的结果瓦片兼容目前绝大多数互 联网地图应用，可以直接投入生产，并且本发明支持断点，一旦发生断电等现象，算法可在 原有工作的基础上继续进行切片，避免了重复工作的问题。
[0006] 针对上述技术问题，本发明提供了一种并行方式栅格影像切片方法。步骤如下：
[0007] 第一步:获取原始栅格影像，设置目标瓦片的级别和进程总数。
[0008] 设定目标瓦片的级别level以及MPI参与的进程总数n，以及瓦片输出路径。MPI会 为每一个进程分配一个进程号i (0 < i〈n，i取整数），后续可通过进程号来达到对每个进程 单独控制的目的。
[0009] 第二步:读取原始栅格影像的元数据信息。
[0010] 设置其中一个进程为主进程，主进程读取原始栅格影像的元数据信息，包括影像 的长宽、波段数，投影坐标系、数据的无效值；
[0011]第三步:将原始栅格影像变换至WebMercator投影。
[0012]当前主流WebGIS软件以及互联网地图应用的切分瓦片坐标系统普遍采用Google 提出的WebMercator投影坐标系，因此在进行切片前首先需将原始栅格影像的坐标系统转 换至WebMercator目标投影系，然后基于WebMercator投影下的坐标进行切分。为了加快投 影变换的效率，本发明采用以下策略:主进程根据之前读取的原始影像投影坐标系，如果原 始栅格影像不是WebMercator投影系(在地理栅格影像中，所有投影坐标系都对应一个由开 放地理空间联盟制定的WKT编码，利用GDAL库的相关函数从原始影像中读取该影像的WTK形 式投影信息，通过比较原始影像的WKT编码和WebMercator投影的WTK就可以判断该影像的 投影系是否为WebMercator投影）。主进程利用⑶AL类库的GDALAutoCreateWarpedVRT函数， 重采样方法采用最邻近内插法，将原始栅格影像变换至WebMercator投影坐标系下，得到投 影变换结果影像，以vrt格式文件进行保存，后续的所有切片操作都是基于投影变换结果影 像进行，并将原始栅格影像的无效值、影像宽度、波段数等元数据信息写入vrt文件中。vrt 是一个xml格式的文件，本身不存储像素的色彩信息，它通过与原始栅格影像相关联，通过 描述信息记录自身像素的坐标（自身像素坐标就是指像素点在投影变换结果影像像素坐标 系下的坐标)与相关联栅格影像的映射关系，以及投影变换结果影像的长宽、波段数等元数 据信息，当通过vrt文件读取栅格数据时，首先基于映射关系，找到其在原始栅格影像文件 中的位置，然后从原始栅格影像文件相应位置中把栅格数据读出。由于vrt文件只记录一些 描述性信息，包括投影变换结果影像的长宽、波段数、像素点的数据类型、投影坐标系，以及 与原始栅格影像的映射关系。所以vrt文件数据量小，因此vrt文件用于作为中间文件的投 影变换结果影像具有很大的优势，可以减少大量10操作。
[0013] 第四步:计算投影变换结果影像的瓦片最大级别和最小级别，并更新目标瓦片的 级别。
[0014] 假设以像素为单位每个瓦片的大小为tilesizeXtilesize，投影变换结果影像长 为Xsize、宽为Ysize，以米(m)为单位投影变换结果影像分辨率为resit。以一张世界地图而 言，Resolution(h)表示其第h级别的瓦片分辨率，Resolution(h)可以通过以下公式计算得 出1^8〇1111:;[011(11) = (2\31\1?)/(1：;[168丨26\211)，1?为地球半径等于6378137111，11为整数，31为 圆周率。投影变换结果影像的瓦片最大级别tmaxz为Resolution (h)最接近res It但不大于 resltB^^hfE； j|?T^minRes= (resit Xmax(Xsize ,Ysize))/tilesize 数，贝丨】瓦片最小级别tminz为1^801111:;[011(11)最接近111；[111^8但不大于111；[111^8时的11值。如果之 前设定的目标瓦片级别level>tmaxz则将tmaxz的值赋给level，如果level〈tminz则将 tminz的值赋给level。而后在瓦片输出目录下创建名为level的文件夹。
[0015] 第五步：计算投影变换结果影像的瓦片行列号范围，并按照路径为/目标瓦片级 另IJ/瓦片列号/瓦片行号.png的路径建立空瓦片文件；
[0016] 计算当前level级别下投影变换结果影像地理范围内瓦片的行列号范围，具体计 算方法如下:假设投影变换结果影像地理范围为[ominX，ominY，omaxX，omaxY]，其中该范围 为一个矩形四边形，omi nX，omaxY为左上角点坐标，omaxX，omi nY为右下角点坐标;瓦片行列 号范围为[tminX，tminY，tmaxX，tmaxY]，其中该范围为一个矩形四边形，tminX，tmaxY为左 上角行列号，tmaxX，tminY为右下角行列号，则其之间满足以下映射关系：
[0017] tminX = ceil( (ominX+originShift)/ (Resolution Xtilesize))-1；
[0018] tminY = cei1((ominY+originShift)/(Resolution Xtilesize))-1；
[0019] tmaxX = ceil( (omaxX+originShift)/ (Resolution Xtilesize))-1；
[0020] tmaxY = cei1((omaxY+