高斯投影的邻带坐标换算汇总.

高斯投影6度和3度分带计算公式高斯投影6度和3度分带计算公式什么是高斯投影6度和3度分带？•高斯投影是一种常用于大地测量和地图制图的投影方法。

根据地球的形状和表面特征，我们将地球划分成了若干个分带，每个分带的宽度为6度或3度。

•6度和3度分带指的是每个分带的经度跨度。

例如，6度分带就是每个分带的中央经线与相邻分带的中央经线之间跨越6度。

高斯投影6度和3度分带计算公式6度分带投影计算公式1.计算投影平面与地球经度的差值：L=λ−L02.计算弧长元素：N=a/√1−e2⋅sin2φ3.计算卯酉圈曲率半径：M=N⋅(1−e2)=a⋅(1−e2)/(1−e2⋅sin2φ)4.计算子午线弧长：A=(1+3e2/4+45e4/64+175e6/256+11025e8/16384)⋅N5.计算坐标系原点到点的子午线弧长：S=A−A06.计算纬度差：t=tanφ7.计算坐标Y轴偏移量：y=x⋅cosφ8.计算坐标X、Y（单位：m）：X=S−N⋅tanφ2⋅L2−N⋅tanφ24⋅(5−t2+9C2+4C4)⋅L4−N⋅tanφ720⋅(61−58t2+t4−270C2+330C4)⋅L6Y=N⋅L⋅cosφ1+N⋅L3⋅cosφ6⋅(1−t2+C2)+N⋅L5⋅cosφ120⋅(5−18t2+t4+14C2−58C4)3度分带投影计算公式1.计算投影平面与地球经度的差值：L=λ−L02.计算弧长元素：N=a/√1−e2⋅sin2φ3.计算卯酉圈曲率半径：M=N⋅(1−e2)=a⋅(1−e2)/(1−e2⋅sin2φ)4.计算子午线弧长：A=(1+3e2/4+45e4/64+175e6/256+11025e8/16384)⋅N5.计算坐标系原点到点的子午线弧长：S=A−A06.计算纬度差：t=tanφ7.计算坐标Y轴偏移量：y=x⋅cosφ8.计算坐标X、Y（单位：m）：X=S−N⋅tanφ2⋅L2+N⋅tanφ24⋅(5+t2+9C2+4C4)⋅L4−N⋅tanφ720⋅(61+90t2+45t4+46C2−252C4−90C6)⋅L6Y=N⋅L⋅cosφ1+N⋅L3⋅cosφ6⋅(1+2t2+C2)+N⋅L5⋅cosφ120⋅(5+28t2+24t4+6C2+8C4)示例解释假设我们需要计算某个点在高斯投影6度分带中的投影坐标。

高斯投影邻带换算一、实验目的编写高斯投影邻带换算的程序，掌握高斯投影邻带换算的基本原理和方法。

二、实验内容：已知高斯平面坐标（x0，y0）及其中央子午线经度L0，计算其在相邻投影带（3度或6度或任意带，中央子午线经度为L1）内的高斯平面直角坐标（x1，y1），转换步骤：1）先利用高斯投影反算将高斯平面坐标（x0，y0）转换为经纬度；2）再利用高斯投影正算转换为相邻投影带内的高斯平面直角坐标（x1，y1）。

创建以GaussianZoneTtrans 命名的函数，函数输入输出格式为：[x1,y1] = GaussianZoneTtrans (x0,y0, L0, L1, a, f)或[x1,y1] = GaussianZoneTtrans (x0,y0, L0, L1, RefEllipsoid) RefEllipsoid为椭球参数RefEllipsoid = [a, b, c, f, e2, e2_];克拉索夫斯基椭球参数：长半轴 a = 6378245扁率 f = 1/298.3b=a(1-f)c = a*a/b;e2 = (a*a-b*b)/(a*a);e2_ = (a*a-b*b)/(b*b);GaussianZoneTtrans函数如下：function [x1,y1] = GaussianZoneTtrans (x,y, L0, L1, a, f)%已知高斯平面坐标（x0，y0）及其中央子午线经度L0，计算其在相邻投影带（3度或6度或任意带，中央子午线经度为L1）内的高斯平面直角坐标（x1，y1）b=a*(1-f) ; %短半轴e=(sqrt(a^2-b^2))/a;%第一偏心率e_=(sqrt(a^2-b^2))/b;%第二偏心率%根据中央子午线弧长x反算底点纬度BfBf = meridian2latitude(x,a,e);Nf=a./sqrt(1-(e.^2).*(sin(Bf).^2));Mf=a.*(1-e.^2)./sqrt((1-(e.^2).*(sin(Bf).^2)).^3);pf=e_.*cos(Bf);tf=tan(Bf);%已知高斯平面坐标（x，y）及指定中央子午线经度L0，计算大地坐标（B,L）B=Bf-tf.*(y.^2)./(2.*Mf.*Nf)+tf.*(5+3.*(tf.^2)+pf.^2-9.*(tf.^2).*(pf.^2)).*(y.^4)./(2 4.*Mf.*(Nf.^3))...+tf.*(61+90.*(tf.^2)+45.*(tf.^4)).*(y.^6)./(720.*Mf.*(Nf.^5));L_=y./(Nf.*cos(Bf))-(1+2.*(tf.^2)+pf.^2).*y.^3./(6.*(Nf.^3).*cos(Bf))...+(5+28.*tf.^2+24.*tf.^4+6.*pf.^2+8.*(tf.^2).*pf.^2).*y.^5./(120.*Nf.^5.*cos(Bf)); %高斯平面坐标转换后的大地坐标（B,L）L=(L_)+L0.*pi./180;%B,L单位：弧度%计算该点在新投影带内的经度之差,L1为新投影带的中央子午线经度l=L-L1.*pi/180;%单位：弧度N=a./(sqrt(1-(e^2).*(sin(B).^2)));t=tan(B);p=e_.*cos(B);%p表示yita%计算子午线弧长Xm0=a*(1-e^2);m2=(3*(e^2)*m0)/2;m4=(5*(e^2)*m2)/4;m6=(7*(e^2)*m4)/6;m8=(9*(e^2)*m6)/8;a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;a2=m2/2+m4/2+15*m6/32+7*m8/16;a4=m4/8+3*m6/16+7*m8/32;a6=m6/32+m8/16;a8=m8/128;X=a0*B-(a2*sin(2*B))/2+(a4*sin(4*B))/4-(a6*sin(6*B))/6+(a8*sin(8*B))/8;%计算新投影带高斯平面坐标（x1，y1）x1=X+(N.*(sin(B)).*(cos(B)).*(l.^2))./2+(N.*(sin(B)).*((cos(B)).^3).*(5-((t).^2)+9. *(p.^2)+4.*(p.^4)).*(l.^4))./24 ...+(N.*sin(B).*(cos(B)).^5.*(61-58.*(t.^2)+(t.^4)+270.*(p.^2)-330.*(p.^2).*(t.^2)).*(l .^6))./720;y1=N.*cos(B).*l+(N.*(cos(B)).^3.*(1-t.^2+p.^2).*(l.^3))./6+(N.*((cos(B)).^5).*(5-1 8.*(t.^2)+ ...t.^4+14.*(p.^2)-58.*(p.^2).*(t.^2)).*(l.^5))./120;y1=y1+500000;%求横坐标自然值要加上500kmend子午线弧长反算函数meridian2latitude如下：function B = meridian2latitude(x,a,e)m0=a*(1-e^2);m2=(3*(e^2)*m0)/2;m4=(5*(e^2)*m2)/4;m6=(7*(e^2)*m4)/6;m8=(9*(e^2)*m6)/8;a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;a2=m2/2+m4/2+15*m6/32+7*m8/16;a4=m4/8+3*m6/16+7*m8/32;a6=m6/32+m8/16;a8=m8/128;%%纬度B的计算B0=x./a0; %B的初始值while 1F=-(a2.*sin(2.*B0))./2+(a4.*sin(4.*B0))./4-(a6.*sin(6.*B0))./6+(a8.*sin(8.*B0))./8;B=(x-F)./a0;if abs(B0-B)<10.^-6break;endabs(B0-B)B0=B;endend2、实例计算验证某点P在北京54坐标系6度带平面直角坐标为：x1=3589644.286my1=20 679 136.438m求P点在3度带第40带的平面直角坐标（x2,y2）。

高斯投影坐大地坐标系由大地基准面和地图投影确定，由地图投影到特定椭圆柱面后在南北两极剪开展开而成，是对地球表面的逼近，各国或地区有各自的大地基准面，我国目前主要采用的基准面为：1.WGS84基准面，为GPS基准面，17届国际大地测量协会上推荐，椭圆柱长半轴a=6378137m，短半轴b=6356752.3142451m；2.西安80坐标系，1975年国际大地测量协会上推荐，椭圆柱长半轴a=6378140m，短半轴b=6356755.2881575m;3.北京54坐标系，参照前苏联克拉索夫斯基椭球体建立,椭圆柱长半轴6度0-6度，角。

值为正，Y增加500公里；反算则是由高斯平面坐标（X，Y）求解大地坐标（L，B）。

二、计算模/***************************************本文直接依据空间立体三角函数关系得出结果。

*****/（一）正算由图表1，由方程式(1),令，可得在图表2中，，则由椭圆方程，令可知：（三、程序代doubleL=(m_L-6.0*L0//换算成弧度doublexita=atan(b*b*tan(B)/a/a/cos(L));doubledxita=0.000001;doublexi=dxita;x=0.0;doublec=a*a/b/b;while(xi<xita){x+=dxita/sqrt(c*sin(xi)*sin(xi)+cos(xi)*cos(xi));xi+=dxita;坐标*&B,do ubledoubledxi=0.000001;doublexi=dxi;doubleX=0.0;doublec=a*a/b/b;while(X<x/a){X+=dxi/sqrt(c*sin(xi)*sin(xi)+cos(xi)*cos(xi));xi+=dxi;}doubler=a/sqrt(c*sin(xi)*sin(xi)+cos(xi)*cos(xi));。

工程测量中高斯-克吕格投影换带计算的应用在测量工程中，有些测区刚好处于投影带边缘，甚至有些工程横跨两个或两个以上投影带，如交通、水利、电力等较长的线路，为了坐标统一的需要，可以进行坐标换带，将相邻带的坐标换成同一系统的数据。

坐标换带有直接换带计算法和间接换带计算法两种，间接换带计算法就是根据第一带的平面坐标某1,y1和中央子午线的经度L1，按高斯-克吕格投影坐标反算公式求得大地坐标B、L，然后根据B,L和第二带的中央子午线经度L2，按高斯-克吕格投影坐标正算公式求得在第二带中的平面坐标某2、y2。

由于在换带计算中，把椭球面上的大地坐标作为过渡坐标，因而称为间接换带法。

这种方法理论上严密，精度高，而且通用性强，虽然计算量较大，但可用电子计算机计算来克服，已成为坐标换带中最根本的方法。

2、换带计算公式用a表示椭球长半轴，b表示椭球短半轴，f=为扁率，e=为第一偏心率，eˊ=为第二偏心率，N=为卯酉圈曲率半径，R=为子午圈曲率半径，B表示经度，L表示纬度。

2.1高斯-克吕格投影反算公式：B=Bf-[-(5+3Tf+Cf-9TfCf)+(61+90Tf+45Tf2)]L=L0+[D-(1+2Tf+Cf)+(5+28Tf+6Cf+8TfCf+24Tf2)]式中：Nf==,Rf=Bf=+(-)in2+(-)in4+in6,e1==,Tf=tg2Bf,Cf=eˊ2co2Bf,D=2.2高斯-克吕格投影正算公式：某N=k0{M+NtgB[+(5-T+9C+4C2)]+(61-58T+T2+270C-330TC)YE=FE+k0N[A+(1-T+C)+(5-18T+T2+14C-58TC)]式中：k0=1,T=tg2B,C=eˊ2co2B,A=(L-L0)coB,N==M=a[(1---)B-(--)in2B+(-)in4B-in6B东西偏移量FE=500000米+带号某10000003.1进行换带计算的步骤分析通常建立独立坐标系的方法是以一个国家控制点和方位角作为起算数据，观测边投影到平均高程面。

§8.3高斯投影坐标正反算公式任何一种投影①坐标对应关系是最主要的；②如果是正形投影，除了满足正形投影的条件外（C-R 偏微分方程），还有它本身的特殊条件。

8.3.1高斯投影坐标正算公式： B,l ⇒ x,y高斯投影必须满足以下三个条件：①中央子午线投影后为直线；②中央子午线投影后长度不变；③投影具有正形性质，即正形投影条件。

由第一条件知中央子午线东西两侧的投影必然对称于中央子午线，即(8-10)式中，x 为l 的偶函数，y 为l 的奇函数；0330'≤l ，即20/1/≈''''ρl ，如展开为l 的级数，收敛。

ΛΛ+++=++++=553316644220l m l m l m y l m l m l m m x （8-33）式中Λ,,10m m 是待定系数，它们都是纬度B 的函数。

由第三个条件知：qyl x l y q x ∂∂-=∂∂∂∂=∂∂, (8-33)式分别对l 和q 求偏导数并代入上式ΛΛΛΛ----=++++++=+++5533156342442204523164253l dqdm l dq dm l dq dm l m l m l m l dqdm l dq dm dq dm l m l m m (8-34) 上两式两边相等，其必要充分条件是同次幂l 前的系数应相等，即ΛΛΛΛΛΛdq dm m dqdm m dqdm m 2312013121⋅=⋅-==(8-35)(8-35)是一种递推公式，只要确定了0m 就可依次确定其余各系数。

由第二条件知:位于中央子午线上的点，投影后的纵坐标x 应等于投影前从赤道量至该点的子午线弧长X ，即(8-33)式第一式中，当0=l时有：0m X x == (8-36) 顾及(对于中央子午线)B V Mr M B N dq dB M dBdXcos cos 2==== 得：B V cB N r dq dB dB dX dq dX dq dm m cos cos 01===⋅===(8-37,38)B B Ndq dB dB dm dq dm m cos sin 22121112=⋅-=⋅-= (8-39)依次求得6543,,,m m m m 并代入(8-33)式，得到高斯投影正算公式6425644223422)5861(cos sin 720)495(cos 24cos sin 2lt t B B N lt B simB N l B B N X x ''+-''+''++-''+''⋅''+=ρηηρρ5222425532233)5814185(cos 120)1(cos 6cos l t t t B N lt B N l B N y ''-++-''+''+-''+''⋅''=ηηρηρρ (8-42) 8.3.2高斯投影坐标反算公式x,y ⇒B,l投影方程：),(),(21y x l y x B ϕϕ== (8-43)满足以下三个条件：①x 坐标轴投影后为中央子午线是投影的对称轴；② x 坐标轴投影后长度不变；③投影具有正形性质，即正形投影条件。

一、用EXCEL进行高斯投影换算从经纬度BL换算到高斯平面直角坐标XY（高斯投影正算），或从XY换算成BL（高斯投影反算），一般需要专用计算机软件完成，在目前流行的换算软件中，存在一个共同的不足之处，就是灵活性较差，大都需要一个点一个点地进行，不能成批量地完成，给实际工作带来许多不便。

笔者发现，用EXCEL可以很直观、方便地完成坐标换算工作，不需要编制任何软件，只需要在EXCEL的相应单元格中输入相应的公式即可。

下面以54系为例，介绍具体的计算方法。

完成经纬度BL到平面直角坐标XY的换算，在EXCEL中大约需要占用21列，当然读者可以通过简化计算公式或考虑直观性，适当增加或减少所占列数。

在EXCEL中，输入公式的起始单元格不同，则反映出来的公式不同，以公式从第2行第1列（A2格）为起始单元格为例，各单元格的公式如下：单元格单元格内容说明A2输入中央子午线，以度.分秒形式输入，如115度30分则输入115.30起算数据L0=INT(A2)+(INT(A2*100)-INT(A2)*100)/60+(A2*10000-INT(A2*100)*100)/3600把L0化成度C2以度小数形式输入纬度值，如38°14’20〃则输入38.1420起算数据BD2以度小数形式输入经度值起算数据LE2=INT(C2)+(INT(C2*100)-INT(C2)*100)/60+(C2*10000-INT(C2*100)*100)/3600把B化成度F2=INT(D2)+(INT(D2*100)-INT(D2)*100)/60+(D2*10000-INT(D2*100)*100)/3600把L化成度G2=F2-B2L-L0H2=G2/57.2957795130823化作弧度I2=TAN(RADIANS(E2))Tan(B)=COS(RADIANS(E2))COS(B)K2=0.006738525415*J2*J2L2=I2*I2M2=1+K2N2=6399698.9018/SQRT(M2)O2=H2*H2*J2*J2P2=I2*J2Q2=P2*P2R2=(32005.78006+Q2*(133.92133+Q2*0.7031))S2=6367558.49686*E2/57.29577951308-P2*J2*R2+((((L2-58)*L2+61)* O2/30+(4*K2+5)*M2-L2)*O2/12+1)*N2*I2*O2/2计算结果XT2=((((L2-18)*L2-(58*L2-14)*K2+5)*O2/20+M2-L2)*O2/6+1)*N2*(H2*J2)计算结果Y表中公式的来源及EXCEL软件的操作方法，请参阅有关资料，这里不再赘述。