专利名称：基于相似变换模型的栅格数据坐标转换并行方法
技术领域：
本发明属于高性能地理计算领域，特别是涉及一种基于相似变换的栅格数据坐标转换的栅格数据处理方法。
背景技术：
栅格数据是以二维矩阵的形式来表示地物或现象的数据组织方式，每个矩阵单位称为一个栅格单元，栅格的每个数据表示地物或现象的属性数据。作为空间数据最基本的组织方式之一，栅格数据在空间数据动态更新、地理信息系统时空数据分析、时空数据共享等方面都发挥这重要的作用。随着对地观测技术的迅速发展，为开展区域地理过程的研究，将获取多尺度、多时相的栅格数据。然而，这些栅格数据的坐标参照系往往会有差异，这给地理数据集成、地理计算、地理过程分析等带来了困难，甚至会引发分析结果的错误，因此有必要对这些栅格数据进行坐标转换。坐标参照系是提供系统原点、尺度、定向及其时间演变的一组协议、算法和常数。坐标系是一种在给定维数空间中用坐标来表示点的方法。坐标转换的目的是将具有不同坐标参考系的数据统一到同一坐标参考系中。坐标参考系的差异主要来自于坐标参考系定义的差异，即基准、坐标原点位置、坐标轴定向和尺度定义的差异。国内外诸多研究机构和学者已开展了大量的坐标转换的研究工作。坐标转换的包括同一基准下和不同基准地心坐标系、大地坐标系、投影坐标系间的转换，其中同一基准下大地坐标系到地心坐标系的转换常采用解析法，地心坐标系到大地坐标的转换多采用近似计算法，数值迭代法和直接法。其中，近似计算法包括BoWringl976年在期刊 Survey Review, 23 卷，181 号发表的“TRANSFORMATION FROM SPATIAL TO GEOGRAPHICAL COORDINATES” 一文中提出的近似计算公式，Toms, R. 1996年在会议Proceedings on Standards for the Interoperability of Distributed Simulation 上的 “An Improved Algorithm for Geocentric to Geodetic Coordinate Conversion，，一文中提出的近似计算公式等；数值迭代法有Laskowski Piotr 1991年在期刊Bulletin Geodesique 第 65,第 I 集中发表的“Is Newton’ s iteration faster than simple iteration for transformation between geocentric and geodetic coordinates，，一文中提出的 Bowring公式的迭代算法，J. Pollard2002年发表在期刊Journal of Geodesy第76卷I 号的 “Iterative vectormethods for computing geodetic latitude and height from rectangular coordinates” 一文中提出的迭代算法等；直接法有H. Vermeille2002年发表在期干丨J Journal of Geodesy 第 76 卷 8 号的“Direct transformation from geocentric coordinates to geodetic coordinates”中提出的直接算法等。不同基准的坐标系间的转换最为常用的为相似变换模型，相似变换是指在各个方向上变换的比率必须相同的一种比例变换，即变换前后图形的形状不变(大小可变)。坐标变换中常涉及的是平移、旋转、缩放三种基本的相似操作，七参数模型(7-parameter Helmert模型或Bursa-Wolf模型等)、三参数模型(Geocentric Translations 等)都属于相似变换模型,EPSG(EPSG: The EuropeanPetroleum Survey Group)在其 Surveying and Positioning Guidance Note Number 7, part 2-Coordinate Conversions and Transformations including Formulas 对坐标变换中涉及的模型赋予了标识码并对模型原理做了详细描述。虽然坐标转换方法的研究较为充分，但现有研究多是针对单机硬件环境，立足于桌面型软件系统，侧重于算法精度、适用性的提升，侧重于单机环境下算法处理效率的研究。面对海量栅格数据，尤其是坐标转换这种高密度计算，为达到栅格数据处理对数据量及速度的要求，应进行并行方法的设计与实现，充分发挥并行集群与多核处理器的高性能计算优势。并行计算是指在并行机上，将一个应用分解成多个子任务，分配给不同的处理器， 各个处理器之间相互协同，并行地执行子任务，从而达到加速求解速度，或者求解应用问题规模的目的。
发明内容
I.发明要解决的技术问题
本发明针对现有的栅格数据坐标转换串行算法模式和传统的硬件平台难以支撑海量栅格数据快速坐标转换的需求的问题，本发明提供了一种基于相似变换模型的栅格数据坐标转换并行方法，该方法在处理的数据量和效率上可以得到比较满意的结果，并具备良好的通用性。2.本发明的技术方案如下
原理在本发明中，栅格数据的坐标转换总体上采用间接法，即首先根据源文件属性及目标参考系等参数生成一个具有目标坐标参考系的初始目标文件，然后计算该初始文件中的每一像元的值。本发明主要采取数据并行策略，由于计算涉及目标文件中每一像元，同时考虑负载均衡，将数据平均划分并行计算。基于相似变换模型的栅格数据坐标转换并行方法，其步骤为(步骤流程图见图I) 步骤I :并行初始化，即完成MPI (Message Passing Interface)的初始化和GDAL
(Geospatial Data Abstraction Library)的格式驱动的注册。本方法使用MPI作为并行编程工具，需要注明并行开始的位置，并获取并行处理的总进程数和当前运行的进程号。此外，本方法是基于GDAL完成栅格的读写操作，需要注册GDAL的格式驱动。步骤2 :用户以命令行形式输入参数，参数包括源文件名及路径、目标文件名及路径、目标坐标参考系。如源文件坐标参考系缺失或错误，则参数还应包括源坐标参考系。坐标参考系的表达方式可以为以下几种中的一种EPSG码，PRO J. 4定义(PRO J. 4 : Cartographic Projections library), WKT 定义(WKT :We 11 Known Text)。步骤3 :调用GDAL数据读写数据函数，将源文件作为GDALDataset类的对象，获取属性信息。如用户自定义了源坐标参考系，则源文件的坐标参考系为用户自定义源坐标参考系。步骤4 :根据源坐标参考系和目标坐标参考系确定坐标转换类型及转换步骤。具体方法为如源坐标参考系与目标参考系的基准相同，坐标系不同，则坐标转换的类型为坐标系转换，坐标系类型有地心坐标系、大地坐标系、投影坐标系，三种坐标系两两间有算法实现转换，转换的步骤即为栅格数据行列坐标转为投影坐标，投影坐标与于其他两类坐标进行转换。如源坐标参考系与目标参考系的基准不相同，则为基准转换，源坐标参考系的数据向目标坐标参考系的转换步骤为源坐标参考系的行列坐标一源坐标参考系的投影坐
标-源坐标参考系的大地坐标-源坐标参考系的地心坐标-WGS84地心坐标-目
标参考系的地心坐标-目标参考系的大地坐标-目标参考系的投影坐标-目标参
考系的行列坐标，目标参考系的数据向源坐标参考系的转换与上述步骤相反。步骤5 :对源文件的边界进行采样，对采样点进行源坐标参考系到目标坐标参考系的坐标变换，坐标转换的步骤已由步骤4确定，坐标转换的原理将在步骤9中详细描述， 获得目标文件的地理范围，目标文件的分辨率与源文件的分辨率相同，目标文件的所有像元的值均为0，计算行列坐标到投影坐标的仿射变换参数集，从而建立一个初始的目标文件。步骤6 :根据通信域中的进程数，对目标文件平均分块，各进程根据其进程编号， 获取相应数据块及其在目标文件中的位置信息和数据块的大小；
步骤7 :各进程对数据块的边界进行采样，对采样点进行目标坐标参考系到源坐标参考系的坐标变换，确定该数据块在源文件中对应的范围，坐标转换的步骤已由步骤4确定， 坐标转换原理将在步骤9中详细描述。步骤8 :读取数据块数据，存入一个长度为m的一维数组(m=数据块行数X数据块列数X波段数)，读取源文件中与数据块对应范围内的数据，存入一个长度为η的一维数组(η=范围内行数X范围内列数X波段数),数组的元素为文件中某一像元在某一波段的值；
步骤9 :对数据块中每一像元的坐标进行由目标参考系到源目标参考系的转换，转换步骤已在步骤4中确定，具体转换原理为(步骤流程图见图2)
1)行列坐标转换为投影坐标根据文件的行列坐标到投影坐标的仿射变换参数集进行仿射变换，公式为
Xp = padfTransform[O] + P^padfTransform[I] + L^padfTransform[2];
Yp = padfTransform[3] + P^padfTransform[4] + L^padfTransform[5];
其中padfTranSform为含6个元素的一维数组，存储行列坐标到投影坐标的仿射变换所需的6个参数；P为像元在栅格文件中所在行号，L为像元在栅格文件中所在列号，行列号均由零开始；Xp, Yp为像元的投影坐标；
2)投影坐标转换为大地坐标，是投影的逆过程；
3)将大地坐标转换为地心坐标,转换公式为
N:
■e sin
X = (AT+//) cos J cos ￡
Y= (N^-H)CosBsm L Z = [iV(l-e2) + //]sin B
其中a为参考椭球的长半轴，6为参考椭球的短半轴，见e为中间参数，沒为大地经度， L为大地纬度，"为大地高，(B, L, H)为大地坐标，(X, Y, Z)为地心坐标；
4)不同基准的地心坐标的转换，转换的参数由用户输入的目标坐标参照系或源坐标参照系的表达式获得，如果用户以EPSG码形式输入目标坐标参考系或源坐标参考系，则转换参数将从EPSG相关的文件gcs. csv、pcs. csv中获得，转换模型为相似变换模型，根据参数
7个数采用相应的三参数模型或七参数模型
三参数模型三参数模型假设转换过程只进行了平移操作，公式为
权利要求
1.基于相似变换模型的栅格数据坐标转换并行方法，其步骤为步骤I :并行初始化，即完成MPI的初始化和GDAL的格式驱动的注册，本方法使用MPI 作为并行编程工具，需要注明并行开始的位置，并获取并行处理的总进程数和当前运行的进程号；此外，本方法是基于GDAL完成栅格的读写操作，需要注册GDAL的格式驱动；步骤2 :用户以命令行形式输入参数，参数包括源文件名及路径、目标文件名及路径、 目标坐标参考系；如源文件坐标参考系缺失或错误，则参数还应包括源坐标参考系；坐标参考系的表达方式可以为以下几种中的一种EPSG码，PR0J.4定义，WKT定义；步骤3 :调用GDAL数据读写数据函数，将源文件作为GDALDataset类的对象，获取属性信息；如用户自定义了源坐标参考系，则源文件的坐标参考系为用户自定义源坐标参考步骤4 :根据源坐标参考系和目标坐标参考系确定坐标转换类型及转换步骤；具体方法为如源坐标参考系与目标参考系的基准相同，坐标系不同，则坐标转换的类型为坐标系转换，坐标系类型有地心坐标系、大地坐标系、投影坐标系，三种坐标系两两间由算法实现转换，转换的步骤即为栅格数据行列坐标转为投影坐标，投影坐标与其他两类坐标进行转换；如源坐标参考系与目标参考系的基准不相同，则为基准转换，源坐标参考系的数据向目标坐标参考系的转换步骤为源坐标参考系的行列坐标-源坐标参考系的投影坐标-源坐标参考系的大地坐标-源坐标参考系的地心坐标-WGS84地心坐标-目标参考系的地心坐标一目标参考系的大地坐标一目标参考系的投影坐标一目标参考系的行列坐标，目标参考系的数据向源坐标参考系的转换与上述步骤相反；步骤5 :对源文件的边界进行采样，对采样点进行源坐标参考系到目标坐标参考系的坐标变换，坐标转换的步骤已由步骤4确定，坐标转换的原理将在步骤9中详细描述，获得目标文件的地理范围，目标文件的分辨率与源文件的分辨率相同，目标文件的所有像元的值均为0，计算行列坐标到投影坐标的仿射变换参数集，从而建立一个初始的目标文件； 步骤6 :根据通信域中的进程数，对目标文件平均分块，各进程根据其进程编号，获取相应数据块及其在目标文件中的位置信息和数据块的大小；步骤7 :各进程对数据块的边界进行采样，对采样点进行目标坐标参考系到源坐标参考系的坐标变换，确定该数据块在源文件中对应的范围，坐标转换的步骤已由步骤4确定， 坐标转换原理将在步骤9中详细描述；步骤8 :读取数据块数据，存入一个长度为m的一维数组，m=数据块行数X数据块列数X波段数，读取源文件中与数据块对应范围内的数据，存入一个长度为η的一维数组，η= 范围内行数X范围内列数X波段数，数组的元素为文件中某一像元在某一波段的值； 步骤9 :对数据块中每一像元的坐标进行由目标参考系到源目标参考系的转换，转换步骤已在步骤4中确定，具体转换原理为I)行列坐标转换为投影坐标根据文件的行列坐标到投影坐标的仿射变换参数集进行仿射变换，公式为Xp = padfTransform[O] + P^padfTransform[I] + L^padfTransform[2];Yp = padfTransform[3] + P^padfTransform[4] + L^padfTransform[5];其中padfTranSform为含6个元素的一维数组，存储行列坐标到投影坐标的仿射变换所需的6个参数；P为像元在栅格文件中所在行号，L为像元在栅格文件中所在列号，行列号均由零开始；Xp, Yp为像元的投影坐标；2)投影坐标转换为大地坐标，是投影的逆过程；3)将大地坐标转换为地心坐标,转换公式为
2.根据权利要求I所述的基于相似变换模型的栅格数据坐标转换并行方法，其特征在于步骤4中利用并行计算集群与多核处理器架构，采用数据并行的策略，在对目标文件的处理中，将目标文件进行分块。
全文摘要
本发明公开了基于相似变换模型的栅格数据坐标转换并行方法，属于高性能地理计算领域。其步骤为步骤1并行初始化；步骤2用户输入参数；步骤3调用GDAL数据读写数据函数；步骤4确定坐标转换类型及转换步骤；步骤5进行源文件边界采样及坐标变换；步骤6对目标文件平均分块；步骤7各进程对数据块的边界进行采样及坐标变换；步骤8读取数据块数据、源文件中与数据块对应范围内的数据，存入数组；步骤9对数据块中每一像元的坐标进行由目标参考系到源目标参考系的转换；步骤10计算像元在各波段的值并赋给数组中该像元点相应元素；步骤11将数据块数组写入目标文件。本发明可有效提高栅格数据坐标转换数据量和效率。
文档编号G06F19/00GK102591622SQ201110441700
公开日2012年7月18日 申请日期2011年12月27日 优先权日2011年12月20日
发明者刘永学, 李满春, 李飞雪, 蒲英霞, 陈冲, 陈振杰, 靳志宾, 黄涛 申请人:南京大学