像平面坐标到经纬度坐标的转换

GEOSPATIALINFORMATION坐标转换是几乎每个遥感、GIS、航测等软件都具备的一个基本功能[1]，在功能代码实现过程中会遇到大规模的矩阵运算问题，例如一幅图像为1024*1024甚至更多个像元，一般在用一些开发软件进行编码的时候显得计算量巨大，从而导致计算运行速度慢的问题，也是大型商用软件最忌讳的问题。

IDL（InteractiveDataLanguage）交互式数据语言是由Kodak公司的全资子公司RSI开发并投向市场的、面向矩阵的且用于数据可视化研究与应用开发的计算机语言。其行业涉及到航天、军事、医学、地

球科学、天文学、商用开发等各个领域。IDL的长期成功来源于它的灵活多样的操作模式。它可以快速地分析数据和实现可视化。它还拥有强大的程序环境以及成熟的跨平台终端应用商业软件产品，例如高级遥感影像处理系统（ENVI）就是由IDL编写的。

大型的数组矩阵运算是IDL的强势，在IDL中可以为任何IDL数据类型创建数组。紧凑的数组语法能保证数组运行时不用循环的操作。另外，数组操作的优点充分体现在运行速度上，将数组作为整体进行数组操作，其速度远远大于对数组元素的循环操作。一

收稿日期：2009-06-02

项目来源：国家863计划资助项目（2006AA12Z113）；国防科技工业民用专项科研技术研究资助项目（KJSX0601）；

地理空间信息

GEOSPATIALINFORMATION个成为高效的IDL程序员的诀窍就是灵活有效地利用数组[2]。

选择IDL做为开发工具去实现坐标转换这一涉及大量矩阵数组的运算不仅能大大加快运行速度，而且代码通俗易读。

1算法分析设计

1.1算法分析

目前，大多数的坐标转换都是用四参数或者七参数方法实现的，传统的方法适用于大多数情况。而本文所要做的像平面坐标到经纬度坐标的转换不需要用这种传统的方法，而是用一种简单的数学思维去实现此功能一种方法。[3-4]

在进行算法设计之前，应该考虑如下几个问题：1）地球是不规则的球体，两极扁，赤道鼓。不能简单地做为圆球体来计算；[5]

2）每纬度在地理位置上的距离是一致的（约111km），但是每经度之间的距离是不同的（越往两极，东西方向上的距离就越小），不能直接转换为直角坐标系；

3）给出的经纬度坐标发生变化后，需要动态的缩放所显示地图的大小。1.2算法流程设计

算法结构图如图1所示,其中输入参数有lat0,lon0,m,n,col,row,divisor:lat0代表已知点的纬度；lon0代表已知点的经度；m代表已知点在像平面坐标系中的行坐标；n代表代表已知点在像平面坐标系中的列坐标；col代表图像的像元列数；row代表图像的像元行数；divisor代表图像的分辨率。返回的是一个col*row的结构体数组，结构体中每个字段都是col*row的数组。

图1算法基本流程

2算法的详细设计

1）创建假定的经纬度坐标系统。分别创建一个代表纬度的col*row数组和一个代表经度的row*col数

Feb.,2010Vol.8,No.1

组，在这里所建立的这两个数组构成了一个假定的像平面坐标系。而这个假定的坐标系统就是我们以后计算的核心数据。

建立的两个数组如图2所示，假定的是col=100,row=100,即假定图像为100*100的。这两个数组看起来虽然简单，但却蕴含了北半球的经纬度的走向规律，以后作为运算的一部分。创建经纬度数组的代码如下：(此代码只能实现行列数相同的情况，如果图像行列数不相同，可都定义为行和列数较大的那个数来建立方阵，然后根据需要的col列row行截取出来经纬度数组，此代码略。

)图2建立的索引数组（左为纬度数组，右为经度数组）

IDL>arrlat=indgen(col);建立一个有col个元素的col*1的索引数组

IDL>arrlat=reform(arrlat,1,col);变为1*col的索引列数组

IDL>arrlat=rebin(arrlat,col,row,/sample);扩展为col*row的索引数组

IDL>arrlon=transpose(arrlat);转置上面的数组，经度数组创建完毕

IDL>arrlat=rotate(arrlat,2);旋转数组，纬度数组创建完毕

由此两个数组确定了一个假定的经纬度坐标系，如图3所示。

图3

假定的经纬度坐标系

2）计算所有像元点到横轴和纵轴方向上的地理距离。

IDL>Dx=(arrlat-m)*divisor

Feb.,2010Vol.8,No.1

地理空间信息

GEOSPATIALINFORMATION可以明显地看出，当矩阵运算不是很大的时候，数组整体运算同For循环运算比较优势不是很大，但是当矩阵运算越来越大的时候，数组整体运算的优势性就明显表现出来了，速度要比For循环操作快了10

66

地理空间信息

GEOSPATIALINFORMATIONFeb.,2010Vol.8,No.1

表2

精度的测试

转换后的经纬度坐标(35°22′07.23″N，119°

34′22.21″E)(35°22′46.39″N，119°

33′56.48″E)(35°22′24.21″N，119°

33′17.27″E)(39°22′34.63″N，119°

32′59.56″E)

多倍甚至更多，且随着运算量的增加而递增。3.2精度的验证

特征点位GoogleEarth的经纬度坐标P1P2P3P4

(35°22′06.81″N，119°

34′21.76″E)(35°22′45.58″N，119°

33′55.27″E)(35°22′23.50″N，119°

33′16.84″E)(35°22′33.95″N，119°

32′59.01″E)

图4

上为经过坐标转换的航空摄影正射图像

下为有GoogleEarth的QuickBrid影像

如图4所示，准备两幅有相同区域的遥感影像，一幅是经过坐标转换的航空摄影正射影像（中国科学院遥感应用研究所于2008年12月在烟台航空摄影所得），另外一幅有正确的经纬度坐标的QuickBrid影像。从每幅影像中均匀地挑出几个特征明显的点进行比较，比对结果如表2所示。