1、地图投影算法5.1 地图投影的概念 计算机地图制图时，是将分布在地球上的事物符号化，并绘制到计算机屏幕或者纸上。地球上的地物使用经纬度坐标表定其位置，地球的形状是一个扁率很小的近似椭球体。如果把这个椭球体缩小几百万倍，它的形状就接近正球体。地球表面是曲面，计算机屏幕或者纸张是平面的，就象要将一个橘子剥开来摊平，理论上讲橘皮永远是铺不平的。要把地球表面上的物体和现象绘制到平面图纸上，就必须解决球面与平面之间的矛盾。如何使球面转化到平面后少破裂、少重叠、变形小，并使变形符合一定要求，需要一种科学方法，这种方法就是研究从地球面转绘到平面上所采用的数学方法，叫做地图投影。 16世纪航海事业得到很大发展
2、，欧洲人发现了两个新大陆；麦哲伦船队的一只船成功地完成了环球航行。那时，航海家粗略地知道各大陆的位置，并且能借助于罗盘航行，但是他们无法确定以何种航线能达到目的地。为了能在航海图上随时确定航向，比利时地图学家和数学家墨卡托发明的地图投影。 地球面上任意一点的位置是由经纬度来决定的。所谓地图投影就是建立平面上的点（用平面直角坐标或极坐标表示）和地球表面上的点（用纬度和经度表示）之间的函数关系。可用数学公式表达： x = f1 (, ) y = f2 (, ) (5.1) 如果能建立x、y和、之间的函数式，就可以依据地表上的点（，），求出其在平面上的位置（x,y）。 地图投影会存在两中误差，形状变
3、化(也称角度变化)或者面积变化。投影以后能保持形状不变化的投影，称为等角投影 (Conformal mapping)，它的优点除了地物形状保持不变以外，在地图上测量两个地物之间的角度也能和实地保持一致，这非常重要，当在两地间航行必须保持航向的准确；或者另外一个例子是无论长距离发射导弹还是短距离发射炮弹，发射角度必须准确测量出来。因此等角投影是最常被使用的投影。等角投影的缺点是高纬度地区地物的面积会被放大。在墨卡托投影的地图上阿拉斯加和巴西看上去差不多同样大小，但是事实上巴西要大四倍。如图5-6。 投影以后能保持形状不变化的投影，称为等面积投影 (Equivalent mapping)，在有按面
4、积分析需要的应用中很重要，显示出来的地物相对面积比例准确，但是形状会有变化，假设地球上有个园，投影后绘制出来即变成个椭圆了。还有第三种投影，非等角等面积投影，意思是既有形状变化也有面积变化，这类投影既不等角也不等积，长度、角度、面积都有变形。其中有些投影在某个主方向上保持长度比例等于1，称为等距投影。 地图投影的过程是可以想象用一张足够大的纸去包裹地球，将地球上的地物投射到这张纸上。地球表面投影到平面上、圆锥面或者圆柱面上，然后把圆锥面、圆柱面沿母线切开后展成平面。如图5-1所示，根据这张纸包裹的方式，地图投影有可以分成：方位投影、圆锥投影和圆柱投影。根据这张纸与地球相交的方式，地图投影有可以
5、分成切投影和割投影，在切线或者割线上的地物是没有变形的，而距离切线或者割线越远变形越大。 还有不少投影直接用解析法得到。根据所借助的几何面不同可分为伪方位投影、伪圆锥投影、伪圆柱投影等。 总结一些常用投影的适用范围。等角投影常用于航海图、风向图、洋流图等。现在世界各国地形图采用此类投影比较多。等积投影用于绘制经济地区图和某些自然地图。对于大多数数学地图和小比例尺普通地图来说，应优先考虑等积的要求。地理区域，诸如国家、水域和地理分类地区（植被、人口、气候等）相对分布范围，显然是十分重要的内容。 任意投影常用作数学地图，以及要求沿某一主方向保持距离正确的地图。常用作世界地图的投影有墨卡托投影、高尔
6、投影、摩尔威特投影、等差分纬线多圆锥投影、格灵顿投影、桑森投影、乌尔马耶夫投影等。 5.2 地图投影的实现 5.2.1 扩充地图层对象类 继承前一章的地图层对象类TGeoLayer派生出一个新的类TGeoLayerWithProjction，使它具有投影坐标计算和转换的功能。 TGeoLayerWithProjection = class(TGeoLayer) end; TGeoLayerWithProjction中新增域ProjectionType以记录投影的类型。本书共设有12种投影选择。见5.2.2。 投影坐标计算的功能通过方法SetProjectionType实现。它首先释放地图层原先
7、拥有的所有地图对象，重新从磁盘文件中读入原始数据，然后判别ProjectionType，当ProjectionType为0时不作投影处理，否则调用投影计算过程ty完成坐标的转换。ty将根据ProjectionType分别调用11种不同的投影计算过程，这些投影计算过程在下一节中介绍。坐标转换后计算每个地图对象最大矩形范围和地图层的最大矩形范围 procedure TGeoLayerWithProjection.SetProjectionType(prjType: integer); end; 地图投影过程ty根据ProjectionType分别调用不同的投影计算过程。 procedure ty(
8、prjType: integer; xys: PXYArray; nData: integer); var ii: integer; begin case prjType of 1: for ii := 0 to nData - 1 do with xysii do ty1(x, y, x, y); 2: for ii := 0 to nData - 1 do with xysii do ty2(x, y, x, y); 3: for ii := 0 to nData - 1 do with xysii do ty3(x, y, x, y); 4: for ii := 0 to nData -
9、 1 do with xysii do ty4(x, y, x, y); . end; end; ty过程的另外一种写法是： procedure ty(prjType: integer; xys: PXYArray; nData: integer); begin for ii := 0 to nData - 1 do with xysii do case prjType of 1: ty1(x, y, x, y); 2: ty2(x, y, x, y); 3: ty3(x, y, x, y); 4: ty4(x, y, x, y); . end; end; 从书写形式上第二种写法更加简洁，但是
10、分析运行效率，第一种写法效率更高，运算速度更快，因为第一种写法作一次case判断和nData次投影计算过程的调用，而第二种写法作nData次case判断和nData次投影计算过程的调用。 5.2.2 地图投影算法 地图投影计算中，常常用到一些与地球有关的参数，在介绍地图投影算法前先介绍这些常用的参数： 原始经纬度：、 经线弧长：s = 6367558 *-16000*sin(2*)+20*sin(4*)-0.022*sin(6*)+0.00003*sin(8*) 投影系数：n = 6367558/sqrt(1-0.00669*sqr(sin()， 纬线半径：r = n*cos() 极坐标的极径
11、和方位角：、，极坐标向平面坐标X、Y的转换公式为， X =cos() Y =sin() 天顶距和方位角：z 、A 地球半径：eR = 6378245 比例尺：Scale = 100000 下面例举11种常见的地图投影。为了保证精度，投影的过程是先作投影再作比例变化。 1) 正轴等角方位投影（球面投影） 正轴等角方位投影，视点位于地球表面，又称球面投影。这种投影的变形在中心区较小，离中心点越远，变形越大，等变形线是以极点为圆心的同心圆。适于做圆形地区的地图。常用于编制两极地区的地图。如图5-2。 公式： =2Rtg(z/2) =A A = Z=/2- (5.2) X= cos() Y=sin()
12、 算法： procedure ty1(Lo,La:single; var x,y:single); begin . end; 2) 正轴等积方位投影（兰勃特投影） 正轴等积方位投影使的投影图上面积和实际地面相应面积相等，常用于小比例尺的地图中。如图5-3。 公式： =2Rsin(z/2) =A A= Z=/2- (5.3) X = cos() Y=sin() 算法： procedure ty2(Lo,La:single; var x,y:single); . end; 3) 正轴等距方位投影（波斯托投影） 正轴等距方位投影，垂直圈投影后保持长度不变，变形较前两者为小，主要用于两极地区，也用于航
13、行图。当将投影中心选在机场，到各地的距离与方位角均与实地相等。如图5-4。 公式： = RZ =A A= Z=/2- (5.4) X =cos() Y =sin() 算法： procedure ty3(Lo,La:single; var x,y:single); . end; 4) 等角正轴切圆柱投影（墨卡托投影） 等角正轴切圆柱投影由荷兰制图学家墨卡托（Mercator Gerardus）于1569年所创，故又名墨卡托投影。墨卡托投影投影变形是随纬度变化而变化，在同一纬线上各点的长度比、面积比、角度变形均相等，适于做适于纬线方向延伸地区。尤其是赤道附近。而且等角航线投影为直线，主要用于海图制
14、作。如图5-5。 公式: x = Rln(tg(/4+/2) (5.5) y = R 算法： procedure ty4(Lo,La:single; var x,y:single); . end; 5) 等角横切椭圆柱投影（高斯-克吕格投影）（Gauss_Krivger） 横切椭圆柱投影，由德国数学家、物理学家、天文学家高斯于19世纪20年代拟订，后经德国大地测量学家克吕格对投影公式加以补充得到。其经纬线互相垂直，变形最大位于赤道与投影带最外一条经线的交点上，常用于纬度较高地区。如图5-6。 公式： s = 6367558-16000sin(2)+20sin(4)-0.022sin(6)+0.
15、00003sin(8) = e2coscos 6367558 n = - 1-0.00669sin2 x = s+n/2 2sincos+n/24 4sincos3 (5-tg2+9+42)+n/720 6sincos5(61-58tg2+tg4+2702-330tg2) + (5.6) y = ncos+n/6 3cos3(1-tg2+)+n/120 5cos5(5-18tg2+tg4+14-58tg2)+ 算法： procedure ty5(Lo,La:single; var x,y:single); . end; 通用横轴墨卡托投影（U.T.M）在高斯-克吕格投影的基础上乘上0.9996
16、，以减少中低纬地区的变形。 6) 经线为正弦曲线的等积伪圆柱投影（桑生投影）(Sanson) 桑生投影是一种经线为正弦曲线的等面积伪圆柱投影，纬线为间隔相等的平行直线，经线为对称于中央经线的正弦曲线，在每一条纬线上经线间隔相等。赤道与中央经线没有变形，适于赤道附近南北延伸地区。如图5-7。 公式： y = Rcos (5.7) x = R 算法： procedure ty6(Lo,La:single;Var x,y:single); . end; 7) 经线为椭圆的等积伪圆柱投影（摩尔威特投影） 为在椭圆中获取等积性质，高纬度地区必须缩小南北方向上的比例尺，低纬则增大，东西方向基本正确，最小两
17、个变形区位于中纬度地区，摩尔威特投影常用于关心区域为于此时。如图5-8。 公式： x = sqrt(2)Rsin(t) (5.8) y = 2sqrt(2)R/ Cos(t) 2t + Sin(2t) = Sin() 这里为了使算法简单用代替t 算法： procedure ty7(Lo,La:single;Var x,y:single); . end; 8) 爱凯特正弦伪圆柱投影（Eckert） 为改变桑生投影中高纬地区不易表示的缺点而采用“极线”表示两极经线为正弦曲线的投影方式。如图5-9。 公式: X = 2Rt/SQRT(2+) (5.9) 2RCOS2(t/2) Y = - 2+ t + sin(t) =(2+)/2 sin 这里为了使算法简单用代替t 算法： procedure ty8(Lo,La:single; var x,y:single); . end;
