CN102591622B - 基于相似变换模型的栅格数据坐标转换并行方法 - Google Patents

Abstract

本发明公开了基于相似变换模型的栅格数据坐标转换并行方法，属于高性能地理计算领域。其步骤为：步骤1：并行初始化；步骤2：用户输入参数；步骤3：调用GDAL数据读写数据函数；步骤4：确定坐标转换类型及转换步骤；步骤5：进行源文件边界采样及坐标变换；步骤6：对目标文件平均分块；步骤7：各进程对数据块的边界进行采样及坐标变换；步骤8：读取数据块数据、源文件中与数据块对应范围内的数据，存入数组；步骤9：对数据块中每一像元的坐标进行由目标参考系到源目标参考系的转换；步骤10：计算像元在各波段的值并赋给数组中该像元点相应元素；步骤11：将数据块数组写入目标文件。本发明可有效提高栅格数据坐标转换数据量和效率。

Description

基于相似变换模型的栅格数据坐标转换并行方法

技术领域

本发明属于高性能地理计算领域，特别是涉及一种基于相似变换的栅格数据坐标转换的栅格数据处理方法。

背景技术

栅格数据是以二维矩阵的形式来表示地物或现象的数据组织方式，每个矩阵单位称为一个栅格单元，栅格的每个数据表示地物或现象的属性数据。作为空间数据最基本的组织方式之一，栅格数据在空间数据动态更新、地理信息系统时空数据分析、时空数据共享等方面都发挥这重要的作用。

随着对地观测技术的迅速发展，为开展区域地理过程的研究，将获取多尺度、多时相的栅格数据。然而，这些栅格数据的坐标参照系往往会有差异，这给地理数据集成、地理计算、地理过程分析等带来了困难，甚至会引发分析结果的错误，因此有必要对这些栅格数据进行坐标转换。坐标参照系是提供系统原点、尺度、定向及其时间演变的一组协议、算法和常数。坐标系是一种在给定维数空间中用坐标来表示点的方法。坐标转换的目的是将具有不同坐标参考系的数据统一到同一坐标参考系中。坐标参考系的差异主要来自于坐标参考系定义的差异，即基准、坐标原点位置、坐标轴定向和尺度定义的差异。

国内外诸多研究机构和学者已开展了大量的坐标转换的研究工作。坐标转换的包括同一基准下和不同基准地心坐标系、大地坐标系、投影坐标系间的转换，其中同一基准下大地坐标系到地心坐标系的转换常采用解析法，地心坐标系到大地坐标的转换多采用近似计算法，数值迭代法和直接法。其中，近似计算法包括Bowring1976年在期刊Survey Review,23卷,181号发表的“TRANSFORMATION FROM SPATIAL TO GEOGRAPHICAL COORDINATES”一文中提出的近似计算公式，Toms,R.1996年在会议Proceedings on Standards for theInteroperability of Distributed Simulation上的“An Improved Algorithm for Geocentric to GeodeticCoordinate Conversion”一文中提出的近似计算公式等；数值迭代法有Laskowski Piotr1991年在期刊Bulletin Géodésique第65,第1集中发表的“Is Newton’s iteration faster than simpleiteration for transformation between geocentric and geodetic coordinates?”一文中提出的Bowring公式的迭代算法，J.Pollard2002年发表在期刊Journal of Geodesy第76卷1号的“Iterativevectormethods for computing geodetic latitude and height from rectangular coordinates”一文中提出的迭代算法等；直接法有H.Vermeille2002年发表在期刊Journal of Geodesy第76卷8号的“Direct transformation from geocentric coordinates to geodetic coordinates”中提出的直接算法等。不同基准的坐标系间的转换最为常用的为相似变换模型，相似变换是指在各个方向上变换的比率必须相同的一种比例变换，即变换前后图形的形状不变（大小可变）。坐标变换中常涉及的是平移、旋转、缩放三种基本的相似操作，七参数模型（7-parameter Helmert模型或Bursa-Wolf模型等）、三参数模型（Geocentric Translations等）都属于相似变换模型，EPSG（EPSG:The European Petroleum Survey Group）在其Surveying and Positioning Guidance NoteNumber7,part2：Coordinate Conversions and Transformations including Formulas对坐标变换中涉及的模型赋予了标识码并对模型原理做了详细描述。

虽然坐标转换方法的研究较为充分，但现有研究多是针对单机硬件环境，立足于桌面型软件系统，侧重于算法精度、适用性的提升，侧重于单机环境下算法处理效率的研究。面对海量栅格数据，尤其是坐标转换这种高密度计算，为达到栅格数据处理对数据量及速度的要求，应进行并行方法的设计与实现，充分发挥并行集群与多核处理器的高性能计算优势。并行计算是指在并行机上，将一个应用分解成多个子任务，分配给不同的处理器，各个处理器之间相互协同，并行地执行子任务，从而达到加速求解速度，或者求解应用问题规模的目的。

发明内容

1.发明要解决的技术问题

本发明针对现有的栅格数据坐标转换串行算法模式和传统的硬件平台难以支撑海量栅格数据快速坐标转换的需求的问题，本发明提供了一种基于相似变换模型的栅格数据坐标转换并行方法，该方法在处理的数据量和效率上可以得到比较满意的结果，并具备良好的通用性。

2.本发明的技术方案如下：

原理：在本发明中，栅格数据的坐标转换总体上采用间接法，即首先根据源文件属性及目标参考系等参数生成一个具有目标坐标参考系的初始目标文件，然后计算该初始文件中的每一像元的值。本发明主要采取数据并行策略，由于计算涉及目标文件中每一像元，同时考虑负载均衡，将数据平均划分并行计算。

基于相似变换模型的栅格数据坐标转换并行方法，其步骤为（步骤流程图见图1）：

步骤1：并行初始化，即完成MPI（Message Passing Interface）的初始化和GDAL（GeospatialData Abstraction Library）的格式驱动的注册。本方法使用MPI作为并行编程工具，需要注明并行开始的位置，并获取并行处理的总进程数和当前运行的进程号。此外，本方法是基于GDAL完成栅格的读写操作，需要注册GDAL的格式驱动。

步骤2：用户以命令行形式输入参数，参数包括源文件名及路径、目标文件名及路径、目标坐标参考系。如源文件坐标参考系缺失或错误，则参数还应包括源坐标参考系。坐标参考系的表达方式可以为以下几种中的一种：EPSG码，PROJ.4定义（PROJ.4：CartographicProjections library），WKT定义（WKT：Well Known Text）。

步骤3：调用GDAL数据读写数据函数，将源文件作为GDALDataset类的对象，获取属性信息。如用户自定义了源坐标参考系，则源文件的坐标参考系为用户自定义源坐标参考系。

步骤4：根据源坐标参考系和目标坐标参考系确定坐标转换类型及转换步骤。具体方法为：如源坐标参考系与目标参考系的基准相同，坐标系不同，则坐标转换的类型为坐标系转换，坐标系类型有地心坐标系、大地坐标系、投影坐标系，三种坐标系两两间有算法实现转换，转换的步骤即为栅格数据行列坐标转为投影坐标，投影坐标与于其他两类坐标进行转换。如源坐标参考系与目标参考系的基准不相同，则为基准转换，源坐标参考系的数据向目标坐标参考系的转换步骤为：源坐标参考系的行列坐标——源坐标参考系的投影坐标——源坐标参考系的大地坐标——源坐标参考系的地心坐标——WGS84地心坐标——目标参考系的地心坐标——目标参考系的大地坐标——目标参考系的投影坐标——目标参考系的行列坐标，目标参考系的数据向源坐标参考系的转换与上述步骤相反。

步骤5：对源文件的边界进行采样，对采样点进行源坐标参考系到目标坐标参考系的坐标变换，坐标转换的步骤已由步骤4确定，坐标转换的原理将在步骤9中详细描述，获得目标文件的地理范围，目标文件的分辨率与源文件的分辨率相同，目标文件的所有像元的值均为0，计算行列坐标到投影坐标的仿射变换参数集，从而建立一个初始的目标文件。

步骤6：根据通信域中的进程数，对目标文件平均分块，各进程根据其进程编号，获取相应数据块及其在目标文件中的位置信息和数据块的大小；

步骤7：各进程对数据块的边界进行采样，对采样点进行目标坐标参考系到源坐标参考系的坐标变换，确定该数据块在源文件中对应的范围，坐标转换的步骤已由步骤4确定，坐标转换原理将在步骤9中详细描述。

步骤8：读取数据块数据，存入一个长度为m的一维数组（m=数据块行数×数据块列数×波段数），读取源文件中与数据块对应范围内的数据，存入一个长度为n的一维数组（n=范围内行数×范围内列数×波段数），数组的元素为文件中某一像元在某一波段的值；

步骤9：对数据块中每一像元的坐标进行由目标参考系到源目标参考系的转换，转换步骤已在步骤4中确定，具体转换原理为（步骤流程图见图2）：

1）行列坐标转换为投影坐标：根据文件的行列坐标到投影坐标的仿射变换参数集进行仿射变换，公式为：

Xp=padfTransform[0]+P*padfTransform[1]+L*padfTransform[2];

Yp=padfTransform[3]+P*padfTransform[4]+L*padfTransform[5];

其中：padfTransform为含6个元素的一维数组，存储行列坐标到投影坐标的仿射变换所需的6个参数；P为像元在栅格文件中所在行号，L为像元在栅格文件中所在列号，行列号均由零开始；Xp,Yp为像元的投影坐标；

2）投影坐标转换为大地坐标，是投影的逆过程；

3）将大地坐标转换为地心坐标，转换公式为：

$\begin{array}{c}N=\frac{a}{\sqrt{1-e^2{\sin }^{2}B}}\\ e^2=\frac{a^2-b^2}{a^2}\end{array}$ $\left\{\begin{array}{c}X=(N+H)\cos B\cos L\\ Y=(N+H)\cos B\sin L\\ Z=[N(1-e^2)+H]\sin B\end{array}\right.$

其中：a为参考椭球的长半轴，b为参考椭球的短半轴，N，e为中间参数，B为大地经度，L为大地纬度，H为大地高，(B,L,H)为大地坐标，(X,Y,Z)为地心坐标；

4）不同基准的地心坐标的转换，转换的参数由用户输入的目标坐标参照系或源坐标参照系的表达式获得，如果用户以EPSG码形式输入目标坐标参考系或源坐标参考系，则转换参数将从EPSG相关的文件gcs.csv、pcs.csv中获得，转换模型为相似变换模型，根据参数个数采用相应的三参数模型或七参数模型：

三参数模型：三参数模型假设转换过程只进行了平移操作，公式为：

$\left[\begin{array}{c}{X}_{\mathrm{new}}\\ Y_{\mathrm{new}}\\ Z_{\mathrm{new}}\end{array}\right]=\left[\begin{array}{c}{X}_{\mathrm{original}}\\ Y_{\mathrm{original}}\\ Z_{\mathrm{original}}\end{array}\right]=\left[\begin{array}{c}\mathrm{\ΔX}\\ \mathrm{\ΔY}\\ \mathrm{\ΔZ}\end{array}\right]$

其中：(X_new,Y_new,Z_new)为转换后地心坐标，(X_original,Y_original,Z_original)为转换前地心坐标，△X，△Y，△Z为平移参数；

七参数模型：七参数模型假设转换过程只进行了旋转、缩放、平移操作，公式为：

$\left[\begin{array}{c}{X}_{\mathrm{new}}\\ Y_{\mathrm{new}}\\ Z_{\mathrm{new}}\end{array}\right]=(1+s)\left[\begin{array}{ccc}1& r_z& -r_y\\ -r_z& 1& r_x\\ r_y& -r_x& 1\end{array}\right]\left[\begin{array}{c}{X}_{\mathrm{original}}\\ Y_{\mathrm{original}}\\ Z_{\mathrm{original}}\end{array}\right]+\left[\begin{array}{c}\mathrm{\ΔX}\\ \mathrm{\ΔY}\\ \mathrm{\ΔZ}\end{array}\right]$

其中：(X_new,Y_new,Z_new)为转换后地心坐标，(X_original,Y_original,Z_original)为转换前地心坐标，△X、△Y、△Z为平移参数，r_x、r_y、r_z为旋转参数，s为尺度参数；

5）地心坐标转换为大地坐标，转换采用Toms改进的地心坐标到大地坐标转换算法，该算法是Ralph Toms1996年发表的论文'An Improved Algorithm for Geocentric to GeodeticCoordinate Conversion'中提出，公式为：

B=arctan(Y/X);

W=(X²+Y²)^1/2；aD/c≈1；T₀=Z*(aD/c)；S₀=(T₀ ²+w²)^1/2；sinβ₀=T₀/S₀；cosβ₀=W/S₀；

其中：a为参考椭球的长半轴，b为参考椭球的短半轴，e为参考椭球的偏心率，（B,L,H）为大地坐标，(X,Y,Z)为地心坐标，其余参数均为中间参数；

6）将大地坐标进行投影转换为投影坐标；

7）投影坐标转换为行列坐标，根据文件行列坐标到投影坐标的仿射变换参数及仿射变换公式，反解出投影坐标到行列坐标的变换公式。

步骤10：根据像元在源参考坐标系的坐标，获取该坐标周围像元，并在数组n中找到周围像元在各个波段的值，采用最邻近法，计算像元在各波段的值，并将这些值赋给数组m中该像元相应的元素。

步骤11：利用GDAL读写接口，将数据块数组m写入目标文件。

步骤4中利用并行计算集群与多核处理器架构，采用数据并行的策略，在对目标文件的处理中，将目标文件进行分块，与传统的串行算法相比，处理数据的速度和规模，都有所提高。

3.有益效果

相比于现有技术，本发明的优势在于：

（1）本发明的方法将传统的栅格数据坐标转换算法与新型的硬件建构相结合，提高了计算的速度和规模；

（2）本发明的方法利用并行计算集群与多核处理器架构，采用数据并行的策略，在对目标文件的处理中，将目标文件进行分块，考虑坐标转换针对所有目标文件的像元，与传统的串行算法相比，处理数据的速度和数据量，都有所提高；

（3）本发明结合了效率较高的算法和新型的硬件架构，实践证明，该方法能够有效提高栅格数据坐标转换的数据量和效率，得到比较满意的结果，且能够直接应用于地理信息数据集成与处理。

附图说明

附图1为海量栅格数据坐标转换算法流程图；

附图2为像元坐标转换流程图；

附图3为坐标转换源文件；

附图4为坐标转换目标文件；

附图5为不同通信域进程数时加速比对比图。

具体实施方式

下面结合实施例对本发明做进一步描述。

实施例1

本实例采用附图3所示的栅格数据作为源文件，该文件为遥感影像数据，该数据是由Landsat7号卫星和Landsat5号卫星的主题成像仪(TM)对同一地区的24幅第四波段遥感影像拼接而成的。栅格数据包含31250×34472个30m×30m的像元，数据量为1G。

本实例在Centos Linux6.0操作系统中，采用C++编程语言编程，基于开源地理数据格式转换类库GDAL及PROJ.4实现，使用MPI作为并行编程工具。其中MPI的实现产品选择Openmpi1.4.1，GDAL选择1.8.1版本，PROJ.4选择4.7.0版本。

程序运行选择IBM System x3500-M3X系列服务器环境。服务器硬件配置为CPU2颗，规格为Intel Xeon Quad Core E5620(主频2.4GHz，12MB Cache，四核）；内存为8GB（2根4GB内存条，规格为DDR31333MHz LP RDIMM）；硬盘为2TB（4个500GB硬盘，规格为7.2K6Gbps NL SAS2.5-inch SFF Slim-HS HD），网络为集成的双口千兆以太网。

本实例具体实施步骤如下：

步骤1：MPI并行初始化，获取进程数numprocs和当前进程编号procnum，并注册GDAL格式驱动；

步骤2：用户输入命令行参数：

mpirun–np8hpgc_warp–t_srs EPSG:2436～\data\srcimg.GIF～\data\dstimg.GIF

其中：‘mpirun’为调用MPI应用程序的声明参数，‘-np’、‘–t_srs’为命令行的引导符，‘-np8’表示MPI通信域的进程数为8个，‘–t_srs’引导的‘EPSG:2436’以EPSG码的形式表达目标坐标参考系，为Beijing1954/3-degree Gauss-Kruger CM117E。‘～\data\srcimg.GIF’源文件，‘～\data\dstimg.GIF’为目标文件；

步骤3：调用GDAL数据读写函数，将源文件作为GDALDataset类的对象，获取其行列数、波段数、坐标参考系、分辨率及地理范围等属性数据。本例中，源文件含31250行、34472列像元，单波段，分辨率为30m，坐标参考系为WGS_1984_UTM_Zone_49N，投影坐标的仿射变换参数为（-527463.87,30,0,2687201.17,0，-30）；

步骤4：源坐标参考系与目标坐标参考系基准不同，因此为基准转换。

步骤5：建立初始目标文件。目标文件与源文件分辨率相同，为30m，单波段，坐标参考系为Beijing1954/3-degree Gauss-Kruger CM117E，共有32465×35568个像元（见附图4），对源文件边界采样点进行坐标转换后得到其行列坐标到投影坐标的仿射变换参数为（-1187926.54，30,0,2747037.86,0，-30），所有像元值均为0，使用GDALCreate()函数建立初始目标文件。

步骤6：获取参与计算进程个数，对初始目标文件按行平均分为numprocs块，各进程获取数据块在目标文件中相对于左上角即第0行0列的偏移量为（0，iY），数据块有DstXSize列，nYChunkSize行，核心代码如下：

int procnum;//进程编号,从零开始

int numprocs;//指定通信域的进程数

MPI_Init(&argc,&argv);//MPI初始化函数

MPI_Comm_rank(MPI_COMM_WORLD,&procnum);//获取进程编号

MPI_Comm_size(MPI_COMM_WORLD,&numprocs);//获取进程总数

int nYChunkSize;//数据块的行数

int iY;//数据块在目标文件中的行号；

int DstXSize;//目标文件的列数

int pulx,puly,plrx,plry;

nYChunkSize=ceil(GDALGetRasterYSize(hDstDS)/(double)numprocs);

iY=procnum*nYChunkSize;

if(nYChunkSize₊iY>GDALGetRasterYSize(hDstDS))

    nYChunkSize=GDALGetRasterYSize(hDstDS)-iY;

DstXSize=GDALGetRasterXSize(hDstDS);

步骤7：沿数据块的上、下、左、右边界各取20个采样点，进行坐标转换，转换核心代码将在步骤9中列出，取采样点转换后坐标的最值，得到坐标块在源文件中对应的范围。

步骤8：读入数据。使用函数GDALDatasetRasterIO()读入数据块及源文件中相应范围内的数据，存入相应大小的数组中；

步骤9：坐标转换，核心代码如下：

1）目标文件的行列坐标到投影坐标的仿射变换：

2）投影变换通过GDAL调用PROJ.4中的pj_inv()函数实现；

3）大地坐标转换为地心坐标：

//Latitude纬度，Longitude经度，Height大地高,X、Y为地心坐标

    Sin_Lat=sin(Latitude);

    Cos_Lat=cos(Latitude);

    Sin2_Lat=Sin_Lat*Sin_Lat;

    Rn=gi->Geocent_a/(sqrt(1.0e0-gi->Geocent_e2*Sin2_Lat));

    *X=(Rn+Height)*Cos_Lat*cos(Longitude);

    *Y=(Rn+Height)*Cos_Lat*sin(Longitude);

    *Z=((Rn*(1-gi->Geocent_e2))+Height)*Sin_Lat;

4）三参数模型：

5）七参数模型：

6）地心坐标转换为大地坐标：

步骤10：调用函数GDALDatasetRasterIO()将数据写入目标文件。

在不同进程数条件下，测试方法的加速比（见附图5），结果表明，随着进程数的增加，加速比先变大，随后边小，进程数为8时，加速比达到最大，为3.32。可见，基于相似变换模型的栅格数据坐标转换并行方法在效率上优于传统的坐标转换算法，并具备一定的通用性。

Claims (2)

1.基于相似变换模型的栅格数据坐标转换并行方法，其步骤为：

步骤1：并行初始化，即完成MPI的初始化和GDAL的格式驱动的注册，本方法使用MPI作为并行编程工具，需要注明并行开始的位置，并获取并行处理的总进程数和当前运行的进程号；此外，本方法是基于GDAL完成栅格的读写操作，需要注册GDAL的格式驱动；

步骤2：用户以命令行形式输入参数，参数包括源文件名及路径、目标文件名及路径、目标坐标参考系；如源文件坐标参考系缺失或错误，则参数还应包括源坐标参考系；坐标参考系的表达方式可以为以下几种中的一种：EPSG码，PROJ.4定义，WKT定义；

步骤3：调用GDAL数据读写数据函数，将源文件作为GDALDataset类的对象，获取属性信息；如用户自定义了源坐标参考系，则源文件的坐标参考系为用户自定义源坐标参考系；

步骤4：根据源坐标参考系和目标坐标参考系确定坐标转换类型及转换步骤；具体方法为：如源坐标参考系与目标参考系的基准相同，坐标系不同，则坐标转换的类型为坐标系转换，坐标系类型有地心坐标系、大地坐标系、投影坐标系，三种坐标系两两间由算法实现转换，转换的步骤即为栅格数据行列坐标转为投影坐标，投影坐标与其他两类坐标进行转换；如源坐标参考系与目标参考系的基准不相同，则为基准转换，源坐标参考系的数据向目标坐标参考系的转换步骤为：源坐标参考系的行列坐标——源坐标参考系的投影坐标——源坐标参考系的大地坐标——源坐标参考系的地心坐标——WGS84地心坐标——目标参考系的地心坐标——目标参考系的大地坐标——目标参考系的投影坐标——目标参考系的行列坐标，目标参考系的数据向源坐标参考系的转换与上述步骤相反；

步骤5：对源文件的边界进行采样，对采样点进行源坐标参考系到目标坐标参考系的坐标变换，坐标转换的步骤已由步骤4确定，坐标转换的原理将在步骤9中详细描述，获得目标文件的地理范围，目标文件的分辨率与源文件的分辨率相同，目标文件的所有像元的值均为0，计算行列坐标到投影坐标的仿射变换参数集，从而建立一个初始的目标文件；

步骤6：根据通信域中的进程数，对目标文件平均分块，各进程根据其进程编号，获取相应数据块及其在目标文件中的位置信息和数据块的大小；

步骤7：各进程对数据块的边界进行采样，对采样点进行目标坐标参考系到源坐标参考系的坐标变换，确定该数据块在源文件中对应的范围，坐标转换的步骤已由步骤4确定，坐标转换原理将在步骤9中详细描述；

步骤8：读取数据块数据，存入一个长度为m的一维数组，m=数据块行数×数据块列数×波段数，读取源文件中与数据块对应范围内的数据，存入一个长度为n的一维数组，n=范围内行数×范围内列数×波段数，数组的元素为文件中某一像元在某一波段的值；

步骤9：对数据块中每一像元的坐标进行由目标参考系到源目标参考系的转换，转换步骤已在步骤4中确定，具体转换原理为：

1）行列坐标转换为投影坐标：根据文件的行列坐标到投影坐标的仿射变换参数集进行仿射变换，公式为：

Xp=padfTransform[0]+P*padfTransform[1]+L*padfTransform[2];

Yp=padfTransform[3]+P*padfTransform[4]+L*padfTransform[5];

其中：padfTransform为含6个元素的一维数组，存储行列坐标到投影坐标的仿射变换所需的6个参数；P为像元在栅格文件中所在行号，L为像元在栅格文件中所在列号，行列号均由零开始；Xp,Yp为像元的投影坐标；

2）投影坐标转换为大地坐标，是投影的逆过程；

3）将大地坐标转换为地心坐标，转换公式为：

$\begin{array}{c}N=\frac{a}{\sqrt{1-e^2{\sin }^{2}B}}\\ e^2=\frac{a^2-b^2}{a^2}\end{array}$ $\left\{\begin{array}{c}X=(N+H)\cos B\cos L\\ Y=(N+H)\cos B\sin L\\ Z=[N(1-e^2)+H]\sin B\end{array}\right.$

其中：a为参考椭球的长半轴，b为参考椭球的短半轴，N，e为中间参数，B为大地经度，L为大地纬度，H为大地高，(B,L,H)为大地坐标，(X,Y,Z)为地心坐标；

4）不同基准的地心坐标的转换，转换的参数由用户输入的目标坐标参照系或源坐标参照系的表达式获得，如果用户以EPSG码形式输入目标坐标参考系或源坐标参考系，则转换参数将从EPSG相关的文件gcs.csv、pcs.csv中获得，转换模型为相似变换模型，根据参数个数采用相应的三参数模型或七参数模型：

三参数模型：三参数模型假设转换过程只进行了平移操作，公式为：

$\left[\begin{array}{c}{X}_{\mathrm{new}}\\ Y_{\mathrm{new}}\\ Z_{\mathrm{new}}\end{array}\right]=\left[\begin{array}{c}{X}_{\mathrm{original}}\\ Y_{\mathrm{original}}\\ Z_{\mathrm{original}}\end{array}\right]=\left[\begin{array}{c}\mathrm{\ΔX}\\ \mathrm{\ΔY}\\ \mathrm{\ΔZ}\end{array}\right]$

其中：(X_new,Y_new,Z_new)为转换后地心坐标，(X_original,Y_original,Z_original)为转换前地心坐标，△X，△Y，△Z为平移参数；

七参数模型：七参数模型假设转换过程只进行了旋转、缩放、平移操作，公式为：

$\left[\begin{array}{c}{X}_{\mathrm{new}}\\ Y_{\mathrm{new}}\\ Z_{\mathrm{new}}\end{array}\right]=(1+s)\left[\begin{array}{ccc}1& r_z& -r_y\\ -r_z& 1& r_x\\ r_y& -r_x& 1\end{array}\right]\left[\begin{array}{c}{X}_{\mathrm{original}}\\ Y_{\mathrm{original}}\\ Z_{\mathrm{original}}\end{array}\right]+\left[\begin{array}{c}\mathrm{\ΔX}\\ \mathrm{\ΔY}\\ \mathrm{\ΔZ}\end{array}\right]$

其中：(X_new,Y_new,Z_new)为转换后地心坐标，(X_original,Y_original,Z_original)为转换前地心坐标，△X、△Y、△Z为平移参数，r_x、r_y、r_z为旋转参数，s为尺度参数；

5）地心坐标转换为大地坐标，转换采用Toms改进的地心坐标到大地坐标转换算法，公式为：

B=arctan(Y/X);

W=(X²+Y²)^1/2；aD/c≈1；T₀=Z*(aD/c)；S₀=(T₀ ²+w²)^1/2；sinβ₀=T₀/S₀；cosβ₀=W/S₀；

其中：a为参考椭球的长半轴，b为参考椭球的短半轴，e为参考椭球的偏心率，（B,L,H）为大地坐标，(X,Y,Z)为地心坐标，其余参数均为中间参数；

6）将大地坐标进行投影转换为投影坐标；

7）投影坐标转换为行列坐标，根据文件行列坐标到投影坐标的仿射变换参数及仿射变换公式，反解出投影坐标到行列坐标的变换公式；

步骤10：根据像元在源参考坐标系的坐标，获取该坐标周围像元，并在数组n中找到周围像元在各个波段的值，采用最邻近法，计算像元在各波段的值，并将这些值赋给数组m中该像元相应的元素；

步骤11：利用GDAL读写接口，将数据块数组m写入目标文件。

2.根据权利要求1所述的基于相似变换模型的栅格数据坐标转换并行方法，其特征在于步骤4中利用并行计算集群与多核处理器架构，采用数据并行的策略，在对目标文件的处理中，将目标文件进行分块。

Images