[发明专利]基于相似变换模型的栅格数据坐标转换并行方法

权利要求书

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）将大地坐标转换为地心坐标，转换公式为：

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

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

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

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

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

 

其中：( 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₀；

T₁ = Z +b*e²*sin³β₀；(S₁) ²= T₁² +[w – a* e²*cos³β₀]²；sin²φ₁=T₁²/ S₁²；

cosφ₁=( w – a* e²*cos³β₀)/S₁; R_N = a / (1-e²*sin²φ₁)^1/2;

If sin²φ₁≥sin267.5°then H = Z / sinφ₁ +RN (e²-1).

Else cosφ₁=[ W- a* e²*cos³β₀]/ S₁and H=W/cosφ₁- R_N.

L=arctan(sinφ1 cosφ1/)

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

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

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

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

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

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