用C语言实现地图坐标投影的转换

毕业设（论文水利与生态工程学院

系（院）

测绘工程专业毕业设计（论文）题目

斯-word档可自由复制编辑

高斯-克格影的程序实现TheconversionofGauss-Krugerlanguage毕业设计

word档可自由复制编辑

论文是对用C语言实现地图投影转换的详细论述，前两章是对C语言和地图投影的相关知识进行简单介绍，并将研究地图投影转换的目的、意义和研究现状进行阐述，让我们对相关知识有了深刻认识与了解；后两章是具体实现的程序代码，程序运行实现的截图以及具体数据进行计算对程序进行验证，说明其可行性与正确性，得出最后的结论。关键词：地图投影C语高克吕格投影高斯正反算word档可自由复制编辑

AbstractdissertationistotheIssueofprojectionachievedwithinThefirsttwointroduceustheofprojectionbriefly,andshowthepurpose,situationofstudyofprojectionbywhichgiveusaoftwodescribecode,andthecasesdataintosurecorrectness,drawaMapClanguage;Gauss-Krugerprojection;fortheandword档可自由复制编辑

摘

目录I

..........................................................................Ⅱ第章

引

1.1地图投影简介

1.2究地图投影转换的目的和意义1.3图投影研究概括及趋势综述

11.4究思路和技术方法

1.5C言概述

3第章高斯克格影高投）

2.1斯投影的基本概念

2.2斯投影坐标正反算公式

7第章c言序码

3.1图投影变换的计算程序流程图

133.2C语言程序源代码

163.3程序运行截图

19第章算

结语

22参文

25致谢

26word档可自由复制编辑

第一章

引言1.1地图投简地图投影,Map，按照一定的数学法则，把参考椭球面上的点、线投影到可展面上的方法。地图投影就是指建立地球表面（或其他星球表面或天球面）上的点与投影平面（即地图平面）上点之间的一一对应关系的方法。即建立之间的数学转换公式。它将作为一个不可展平的曲面即地球表面投影到一个平面的基本方法，

保证了空间信息在区域上的联系与完整。这个投影过程将产生投影变形，而且不同的投影方法具有不同性质和大小的投影变形。1.2研地投转的的意地图投影变换是指从一种地图投影点的坐标变换为另一种地图投影点的坐标。图投影总类繁多，形式各异。在实际应用中，地理信息系统的数据来源各不相同，很难保证这些数据的投影方式一致。为了对这些不同来源的地理数据进行有效的处理分析，就必须通过地图投影变换把这些投影方式不同的数据定位在同一种坐标系下，使其具有相同的地图投影。同时，在地图编制过程中，也经常会遇到地图资料与新编地图之间投影不一致的情况，必须进过投影转换将该地图资料转绘到新编地图投影坐标网格中。人类一切经济活动都离不开地理空间，各类专业信息都必须以地形基础信息为空间载体，所以必须研究地图数据库中数字化地图数据处理、空间信息定位和变换，以

地满足各类专业信息系统建设的需要。综上所述，采用

C语言作为软件开发平台，实现各种常见地图投影之间的任意变换，也具有很强的现实意义。1.3地投研概及势述地图投影及其变换模型已渗透到地理信息系统的各个方面，无论是数据的获取、预处理、信息的存贮记录，还是数据的处理、应用和输出都需要有一个空间定位的框架即共同的地理坐标和平面坐标系统。由此可见地图投影在地理信息系统中的重要作用。另外，在地图编制过程中也常遇到原始资料与新编地图投影不一致的情况，此时也需要通过投影变换来解决此类问题。传统的投影变换方法均有工序繁、速度慢、word档可自由复制编辑

精度低等缺点，已不能适应制图自动化的要求。随着计算机技术和数学手段的高速进步，地图的编制以实现了完全自动化。投影变换计算自动化要求建立两种不同投影方式之间点与点的变换关系式，

通过各种关系式编程实现各种投影之间的任意变换，不再需要人工参与变换，样既提高了变换精度又提高了工作效率。目前国内外很多地图制图和地理信息系统软件，例如国外的

ArcgisMapInfo等，国内的MapGis等，都已不同程度地提供了地图投影变换功能模块，因此，采用C语言工具软件实现各种常见的地图投影之间的任意变换，为地图投影变换方法研究提出一套行之有效的方法。1.4究路技方研究思路：1首先弄懂地图投影转换原理。地图投影的实质是用数学的方法实现球面上点到平面上点的转换。

地图投影的分类很多，每种投影都有一套坐标计算公式。“高-克格投影UTM投影兰勃特等角投影”等投影之间的转换公式应熟练掌握。地图投影变换方法主要有以下两类：1）传统地图投影变换方法，包括网格转绘法蓝图拼贴法和纠正仪法等。这些方法都有工序繁杂，速度较慢度较低的缺点，难于男足制图自动化的要求。）数字地图投影，包括解析变换法、数字变换法和数字解析变换法。数字投影变换方法的主要思想就是通过一定的数学手段，

寻找两个不同投影性质地图之间点位一一对应关系式。变换法是根据原始地图的投影方程式反解原始地图投影点的经纬度，再代入新编图的地图投影方程式，得到两种投影的平面直角坐标关系；如建立两种不同地图资料间相应点坐标的直接关系，

即在两种地图上量出相应点的平面直角坐标值，带入逼近多项式，分别组成线性方程组，求出系数值，再代入逼近多项式，即可求出新编图投影点的坐标。第二种方法具有普遍性，适用于地图投影转换。假定原始地图资料上点的坐标为（坐标为（X,Y则对应投影点间坐标转换的基本关系式为

xy新编地图上对应点X=F1(x,y)Y=F2(x,y)皆为制图区域内单值、连续的函数。对于任意两种不同投影性质的地图资料来说，其投影关系式分别为word档可自由复制编辑

由式反解得B=B(x,y)把式代入（式，便有上式即为地图投影转换的数学通式，使用于任何两种投影间变换函数表达式随投影间转换类型和转换方法的不同而不同。2、选择某种投影方式，构建投影变换公式。3、用编程语言实现地图投影转换。1.5C语概

1.4.4）f3f4的具体C语言是一种计算机程序设。它既具有高级语言的特点，又具有汇编语言的特点。它于推出，1978后已先后被移植到大、中、小及微型机上。它既可以作为工作系统设计语言，编写系统应用程序，也可以作为应用程序设计语言，编写不依赖计算机硬件的应用程序。它的应用范围广泛，具备很强的数据处理能力，在软件开发上、各类科研都需要用到C语言，适于编写系统软件，三维，二维图形和动画等。之所以选择使用C语言进行地图投影的转换是因为它有以下特点：简洁紧凑、灵活方便、运算符丰富、数据类型丰富是结构式语言、语法限制不太严格，程序设计自由度大、允许直接访问物理地址，对硬件进行操作、生成目标代码质量高，程序执行效率高、适用范围大，可移植性好等。word档可自由复制编辑

第二章高斯克格投影（高斯投影）2.1高投的本念高斯投影又称横轴椭圆柱等角投影，是德国测量学家高斯于～1830年首先提出的。实际上，直到年，由德国另一位测量学家克吕格推导出实用的坐标投影公式后，这种投影才得到推广，所以该投影又称高斯-克吕格投影。想象有一椭圆柱面横套在地球椭球体外面，并与某一条子午线（称中央子午线或轴子午线）相切，椭圆柱的中心轴通过椭球体中心，然后用一定的投影方法将中央子午线两侧各一定经差范围内的地区投影到椭圆柱面上，再将此柱面展开即成为投影面。图2.1高斯-克吕格投影的几何念（2）分带投影我国规定按经6

0

度进行分带投影为大比例尺测图和工程测量采

0

带投影。特殊情况下工程测量控制网也用1.5带或任意带。高斯投带：0子午线起，每隔经0自西向东分带，依次编号…。我

中央子午线的经度，69

0

起每

至35

0

，共计12带，带号用n示，中央子午线L的经度用表示，word档可自由复制编辑

n

带号和经度的关系）

图2.2高斯-克吕格投影分带示图高斯投

0

带是

带的基础上分成的其中央子午线一部分0

带中央子午线重合一部分则6带分界子午线重合带号用n表示30带中央子午线用L表示关系是：L

。（3）高斯平面直角坐标系在投影面上，中央子午线和赤道的投影都是直线，并且以中央子午线和赤道的交点作为坐标原点，以中央子午线的投影为纵坐标轴，以赤道的投影为横坐标轴，这样便形成了高斯平面直角坐标系在我国

坐标均为正，坐标的最大（在赤道上约为330KM。为避免出现负的横坐标，可在横坐标上加。此外还应在坐标前面冠以带号，这种坐标称为国家统一坐标。如某Y=19123456.789m，该点位19带内，其相对于中央子午线而言的横坐标是：首先去掉带号，再减去500KM最后得y=-376543.211m。（4）高斯投影的特性与优点高斯-克吕格特性（1）等角投影——投影前后的角度相等，但长度和面积有变形；（2）等距投影——投影前后的长度相等，但角度和面积有变形；（3）等积投影——投影前后的面积相等，但角度和长度有变形。高斯平面投影的特点①中央子午线无变形；②无角度变形，图形保持相似；③离中央子午线越远，变形越大。由此可见，在测量中，如果中央子午线输错了，投影的中央子午线就会编离实地坐标系正word档可自由复制编辑

确的中央子午线，变形就越大，最终的结果就使用测量的误差更大。3度带与6度带的划分：1．我国采用6度分带和度分带：1∶2.5及1∶5的地形图采用6度分带投影，即经差为，从零度子午线开始，自西向东每个经差6度为一投影带，全球共分60个带，用1，2，3，，5，……表示．即东经0～6为第一带，其中央经线的经度为东经，东经6～12度为第二带，其中经线的经度为9度1∶1万的地形图采用3度分带，从东1.5度的经线开始，每隔3度为一带，用1，2，3……表示，全球共划分120个投影带，即东经1.5～4.5度为第1带，其中央经线经度为东经3度，东经4.5～7.5度为第带，其中央经线的经度为东经6度。我省处于东经113度至东经120度之间，共跨第、39、此三个带，其中东经度以西便为第38，其中央经线为东经114度；东经～118.5为带，其中央经线为东经；东经118.5以东到山海关为带，其中央经线为东经。地形图上公里网横坐标前2就是带号，例如：1∶5万地形图上的横坐标为20345486，其中20即为带号，345486为横坐标值。2．当地中央经线经度的计算六度带中央经线经度的计算：当地中央经线经度＝6°×当地带号－3°，例如：地形图上的横坐标为20345所处的六度带的中央经线经度为＝117°（适用于1∶2．5万和1∶5万地形图三度带中央经线经度的计算：中央经线经度＝3°×当地带号（适用于∶1万地形图3、怎样计算当地的中央子午线？当地中央子午线取决于当地的直角坐标系统，首先确定您的直角坐标系统是带还是度带投影公式推算:6度带中央子午线计算公式：当地经度/6=N；中央子午线N（带号）如果没有除尽时，即N有余数时，

中央子午线L=6*N-33度带中央子午线计算公式：当地经度/3=N，中央子午线L=3*N我国的经度范围西开始于73°东至135°，可分成六度带十一（13号带—号带中央经线依次（75°81°…123°129°、135°三度带二十二24号—45号带央经线依次…132°六度带可用于中小比例尺（如1：250000）测图，三度带可用于大比例尺（如1：10000）word档可自由复制编辑

2020测图，城建坐标多采用三度带的高斯投影4、如何判断投影坐标是3度带坐标还是6度带坐标如(4231898,21655933)其中21为带号，同样所定义的东伪偏移值也需要加上带号，如21带的东伪偏移值为21500000米。如你的工作区经度在120度至126范围则该坐标系为6度带坐标系，该带的中央经度为123度。如(2949320,36353822)其中36为带号，已知该地点位于贵阳市附近，而从地图上我们看到贵阳大概的经度是东经108度左右因此可以所以该坐标系为3度带坐标系，该带的中央经度为108度。而不可能为6度带：2.2斯影标反公（1）高斯投影坐标正算公式：lx,y高斯投影须满足以下三个条件：①中央子午线投影后为直线；②中央子午线投影后长度不变；③投影具有正形性质，即正形投影条件。由第一条件知中央子午线东西两侧的投影必然对称于中央子午线(8-10)式中l的偶函数，yl奇函数；

l

，即

l

/

，如展开l级数，收敛。lll0246

mll13

3l55

（2-1）式中

,

是待定系数，它们都是纬度B的函数。由第三个条件知：,(2-1)式分别l

和q求偏导数并代入上式mml1

2

l5

4

dml4dq

l

4

2ll24

3

ml6

5

1l3ldqdqdq

l

5

上两式两边相等，其必要充分条件是同次l前的系数应相等，即word档可自由复制编辑

dmmdq1m21m3dq

(2-2)(2-2)是一种递推公式，只要确定了

就可依次确定其余各系数。由第二条件知:位于中央子午线上的点，投影后的纵坐标x应等于投影前从赤道量至该点的子午线弧长X，即(2-1)式第一式中，当

l

时有：xXm

(2-4)顾及(对于中央子午线)

cosrdqM得：dBcm0NBdqdq

(2-5)1dm1dBNmsinBcosB22dBdq

(2-6)依次求得

,m,m,3

6

并代入(2-1)式，得到高斯投影正算公式N

2

NsimBcos3B(5224

l

4Ncos(61t2)lword档可自由复制编辑

33

cosBcos2

)l

cosBt4120

)l

（2）高斯投影坐标反算公式x,yl投影方程：B(x)lx)满足以下三个条件：

(2-7)①x坐标轴投影后为中央子午线是投影的对称轴；②x坐标轴投影后长度不变；③投影具有正形性质，即正形投影条件。①由x求底点纬度(足纬度

B

f

,对应的有底点处的等量纬度

f

，求x,y与l的关系式，fq(yl,y由于y和椭球半径相比较小(1/16.37)可将q必是y的偶函数，l必是y的奇函数。

q,l

展开为y的幂级数又由于是对称投影，2y0lyy13

(2-8)n,,n,由第三条件知：

是待定系数，它们都是x的函数

，

，(2-9)(3-8)式分别对x和y求偏导数并代入上式word档可自由复制编辑

241234f241234fdn02dxdx

nn15

nnyny246

dn1yyy5dx

上式相等必要充分条件，是同次幂y前的系数相等，dn1dndn1n0,n1232

第二条件，当时，点在中央子午线上，x=X，对应的点称为底点，其纬度为底点纬度

B

f

，也就是x=X时的子午线弧长所对应的纬度，设所对应的等量纬度为

f

。也就是在底点展开为y的幂级数。由(2-8)nq

f依次求得其它各系数dndqdBndXff

1NB

f

11Ncosrfff

(2-10)1dBn122ff

tf2N2f

f

(2-11)…………将

,,n04

6

代入(2-8)式得(2-12)

f

t2yfN4cosBf

f

t22ffff24cosBffword档可自由复制编辑

ffdq226ffdq226dq3dqdqdq

f

3fN3f

f

(2-13)n,,将②求B与f

5代入(2-8)得(2-14)式。(最后表达式)y的关系。由(8-7)式

MNcos

dB

知：f(f()ffB()(q)fff按台劳级数q展开f

(2-14)(2-15)1B13Bf(q)dqdq2dq

3

(2-16)Bf

1B

f

f

1B

f

f

（2-17）由此可求出各阶导数：f

f

f

(2-18)2BcosB2fffff

)

(2-19)BB22fffffffffword档可自由复制编辑

(2-20)

55y3yyff……将式(2-13)1,(2-18),(2-19),(2-20)代入2-17)按y幂集合得高斯投影坐标反算公式(2-21),f

tf2f

f

y

tfMNf

3f

222ffff

4

tf720MNf

5f

y612t46ffl2NcosNffff

(2-21)

y5120N5ff

5tt4fffff

归纳由

(y)

求

(,l)

的基本思想：由点

(x,)

得到底点

f(

，将底点f作为过渡，也就是说将坐标原点o到点，先求

qQ(x,f1l(,y)

关系式，再将(x,yf

关系式代入

BBQ(q)f3f

关系式得B(x)f

关系式，最后将坐标原点移回到o点,从而求得

(lword档可自由复制编辑

第三章c语言程序码3.1地图影换计程流图开始输入B,L已知，b已知，b由，b得由L得由L得由BB1由得R由e,B得C由L1,L0,B1得A由得M由得N由L得FE由FE，N，A,T,C得X由，得Y输出word档可自由复制编辑结束

开始输入X,Y已知a，b由，b得由X,Y得L0,X0.Y0由X0，Y0得xval,yval由BB1由ycal得M由U由Bf由a,Bf得C,T由a,e2,BfR由xval,ND由公式得L1,B1由B由L1得L输出B,L结束word档可自由复制编辑

3.2C语程源码高斯正算：//斯正算：#include"stdlib.h"#include"math.h"B,L,X,Y;//B纬度L是度；{请输入坐标（经纬度以小数表示纬度B=");//B是纬度scanf("%lf",&B);getchar();经度L=");//L是经度getchar();//京bbeijing=6356863.0188;//北京//西安bxian=6356755.2882;//西安8084京坐标系的转换函数调用printf("\n北京坐标系的转换结果：Y=%lf\n",X,Y);//打印结果安坐标系的转换函数调用西安坐标系的转换结果：X=%lf印结果sort(awgs,bwgs);//WGS坐标系的转换函数调用坐标系的转结果：打印结果}b0)//函，传值进入{doublea,b;//接受参数即长半轴，短半轴偏率e,e2;//a,b为数，e,e2为第一第二偏心率//心率平方e=1-(b/a)*(b/a);e2=(a/b)*(a/b)-1;dblPI;dblPI=3.1415926/180.0;//角度弧度转换参数int带宽度带或度带printf("\n请输入带宽（度带或是6带）ZoneWide=");word档可自由复制编辑

getchar();intZoneNumber;//带号L0;//中央经度if(ZoneWide==3){//3度带ZoneNumber=(int)((L-ZoneWide/2)/ZoneWide+1);}if{//6度带ZoneNumber=(int)(L/ZoneWide+1);}{输入错误!将退出!");exit(0);}//央经度（用字母示）L0=(ZoneNumber-1)*ZoneWide+ZoneWide/2;//度转换成弧度B1=B*dblPI;L0=L0*dblPI;//午圈曲率半径R=a*(1-e)/sqrt((1-e*sin(B1))*(1-e*sin(B1)*sin(B1))*(1-e*sin(B1)*sin(B1)));T=tan(B1)*tan(B1);//午线弧长M=a*((1-e/4-3*e*e/64-5*e*e*e/256)*B1-(3*e/8+3*e*e/32+45*e*e*e/1024)*sin(2*B1)+(15*e*e//酉圈曲率半径N=a/sqrt(1-e*sin(B1)*sin(B1));//纬度偏移k0=1;//Y=FE+k0*N*(A=(1-T+C)*A*A*A/6+(5-18*T*T*T+14*C-58*T*C)*A*A*A*A*A*A/120);*C)*A*A*A*A*A*A/720);//printf("%lf%lf%lf%lf",B,L,a，bX,Y;}word档可自由复制编辑

高斯反算：//斯反算：#include"stdlib.h"#include"math.h"B,L,X,Y;//B纬度L是度；abeijing,bbeijing,axian,bxian,awgs,bwgs;);{请输入直角坐标：");横坐标getchar();纵坐标getchar();//京bbeijing=6356863.0188;//北京//西安bxian=6356755.2882;//西安80//WGS8484京坐标系的转换函数调用printf("\n直角坐标为北京坐标系时：打印结果安坐标系的转换函数调用printf("\n直角坐标为西安坐标系时：打印结果sort(awgs,bwgs);//WGS坐标系的转换函数调用printf("\n直角坐标为WGS坐标系时：B=%lfL=%lf\n",B,L);//印结果}b0)//函，传值进入{doublea,b;//接受参数即长半轴，短半轴偏率e,ee;//a,b参数，为第一第二偏心率//心率平方e=1-(b/a)*(b/a);dblPI=3.1415926/180.0;//角度弧度转换参数dblPI=3.1415926/180.0;//角度弧度转换参数int带宽度带或度带printf("\n请输入带宽（度带或是6带）ZoneWide=");word档可自由复制编辑

getchar();intZoneNumber;//带号if(ZoneWide==3){//3度带ZoneNumber=(int)(Y/500000);}if{//6度带ZoneNumber=(int)(Y/1000000);}{输入错误!将退出!");exit(0);}//央经度（用字母示）L0=(ZoneNumber-1)*ZoneWide+ZoneWide/2;//度转换成弧度L0=L0*dblPI;X0;Y0=ZoneNumber*1000000+500000;X0=0.00;xval=X-X0;子午线弧长//底点纬度Bf=u+(3*e/2-27*e*e*e/32)*sin(2*u)+(21*e*e/16-55*e*e*e*e/32)*sin(4*u)+(151*e*e*e/96)*sin(6*u)+(1097*e*e*e*e/512)*sin(8*u);C=E*cos(Bf)*cos(Bf);T=tan(Bf)*tan(Bf);//卯酉圈曲率半径底点所对曲率半径计算纬度经度s(Bf);转换为度L=L1/dblPI;B=B1/dblPI;//printf("%lf%lf%lf%lf",B,L,a，b）;;}word档可自由复制编辑

3.3序行图3.3.1高斯算程序运行截图一（1）正算程序代码，进过编译运行后，弹出对话框（2）当分别输入纬度和精度L数据后，选择6带，回车，出现下图结果（3）同理，按任意键继续，再继续输入纬度B经度L数据，选择3度带，回车高斯正算程序运行截图二word档可自由复制编辑

高斯反算）反算与正算类似，输入完程序代码后，经过编译运行后，弹出对话框（5）当分别输入横坐标Y和纵坐标数据之后，选择6度带，回车（6）按任意键继续，再次输入纬度和经度L数据，选择3度带，回车图3.3.3高斯算程序运行截图一图3.3.4高斯反算程序运行截二word档可自由复制编辑

第四章算例正算算例如下表，与之输入相近数据，验证，程序运行结果正确。计算次序

测站符号

D1B2934′2L3L01054l125″5l

5114.86636l=l”0.26161869a2,10^-3225411b2,10^-7544891329708a13220.01614a2l17al+a2l

3220.60614b115b2l18b2l

1425526.91807931216’17’

7770.4935834word档可自由复制编辑

计算次序1238971012514

符号xyyA1,10^-6A1A2yA1+A2yBf’(A1+A2y)

点号

A180853.8683271.029-1.76114.44986-0.0057914.4441029355282