地理坐标转3度或6度分带的算法以及任意多边形求面积的方法(与Arcmap相差0.81%)

 

一：首先补充下理论知识（网上找的）

对于不同的椭球体之间的转换要用到七参数的方法，这个不还没有研究。本文不是讲的这个。

因为本文章主讲算法，理论知识就不多说了。

分带方法
1．我国采用6度分带和3度分带： 　　
         1∶2.5万及1∶5万的地形图采用6度分带投影，即经差为6度，从零度子午线开始，自西向东每个经差6度为一投影带，全球共分60个带，用1，2， 3，4，5，……表示．即东经0～6度为第一带，其中央经线的经度为东经3度，东经6～12度为第二带，其中央经线的经度为9度。 　　
         1∶1万的地形图采用3度分带，从东经1.5度的经线开始，每隔3度为一带，用1，2，3，……表示，全球共划分120个投影带，即东经1.5～ 4.5度为第1带，其中央经线的经度为东经3度，东经4.5～7.5度为第2带，其中央经线的经度为东经6度．比如福建省位于东经113度－东经120度之间， 跨第38、39、40共计3个带，其中东经115.5度以西为第38带，其中央经线为东经114度；东经115.5～118.5度为39带，其中央经线为 东经117度；东经118.5度以东到山海关为40带，其中央经线为东经120度。 地形图上公里网横坐标前2位就是带号，例如：1∶5万地形图上的横坐标为20345486，其中20即为带号，345486为横坐标值。
2．当地中央经线经度的计算
           六度带中央经线经度的计算：当地中央经线经度＝６°×当地带号－３°，例如：地形图上的横坐标为２０３４５，其所处的六度带的中央经线经度为：６°×２０－３°＝１１７°（适用于１∶２．５万和１∶５万地形图）；６°×49－180－３°＝111°（适用于１∶50万以下地形图）。
三度带中央经线经度的计算：中央经线经度＝３°×当地带号（适用于１∶１万地形图）。

       我国x坐标都是正的，y坐标的最大值（在赤道上）约为330km。为了避免出现负的横坐标，可在横坐标上加上500 OOOm。此外还应在坐标前面再冠以带号。这种坐标称为国家统一坐标。

 

二. 转换算法，用C++来实现的

 

 int zonewide=6;     //3度带比较准，可以用6度带
 int base_type=WGS84;   //投影基准面类型：北京54基准面为54，西安80基准面为80，WGS84基准面
 double L,B;//大地经度，纬度

 //----------初始化已知参数-----------
 L = longitude;//经度
 B = latitude;//纬度

 double l,L0;//投影带中央线起算的经差，投影带中央子午线
 double X,Y;//直角坐标
 double f,e2,a,b,e12;参考椭球体扁率,第一偏心率，长半轴,短半轴，第二偏心率
 double t;//保存tan B
 double m;//保存lcosB
 double q2;//卯酉圈的分量 的平方
 double M,N;//子午，卯酉圈曲半径
 double S;//子午线弧长
 double A0,B0,C0,D0,E0;
 double n;//投影带号
 double a_1_e2;//a*(1-e*e)的值

 if(WGS84 == base_type)
 {
  a = 6378137;
  b = 6356752.3142;
  f = 1/298.257223563;
  e2 = 0.0066943799013;//--第一偏心率
  e12 = 0.00673949674227;//--第二偏心率
 }
 else if(BJ54 == base_type)
 {
  a = 6378245;
  b = 6356863.0187730473;
  f = 1/298.3;
  e2 = 0.006693421622966;//--第一偏心率
  e12 = 0.006738525414683;//--第二偏心率
 }
 else if(XIAN80 == base_type)
 {
  a = 6378140;
  b = 6356755.2881575287;
  f = 1/298.257;
  e2 = 0.006694384999588;//--第一偏心率
  e12 = 0.006739501819473;//--第二偏心率
 }

 //求中央经线
 if (zonewide==6)
 {
  n=(int)(L/6)+1;
  L0=n*zonewide-3;
 }
 else
 {
  n=(int)((L-1.5)/3)+1;
  L0=n*3;
 }

 //---基本计算公式 a 长半轴，b 短半轴--------
 //f = (a-b)/a;//--椭球体扁率
 //e2 = (a*a-b*b)/(a*a);//--第一偏心率
 //e12 = (a*a-b*b)/(b*b);//--第二偏心率
 a_1_e2 = a*(1-e2);

 //---转化为弧度,sin,cos 等运算都是以弧度为基础的
 L0=L0*PI/180;
 L=L*PI/180;
 B=B*PI/180;
 l = L - L0;
 t = tan(B);
 q2 = e12*(cos(B)*cos(B));
 N = a/sqrt(1-e2*sin(B)*sin(B));卯酉圈的曲率半径
 m = l*cos(B);

 //--以下算法是对原始公式的多项式进行了合并处理，保存计算更准
 double m0 = a_1_e2;
 double m2 = 3.0 * e2 * m0 / 2.0;
 double m4 = 5.0 * e2 * m2 / 4.0;
 double m6 = 7.0 * e2 * m4 / 6.0;
 double m8 = 9.0 * e2 * m6 / 8.0;
 double a0 = m0 + m2 / 2.0 + 3.0 * m4 / 8.0 + 5.0 * m6 / 16.0 + 35.0 * m8 / 128;
 double a2 = m2 / 2.0 + m4 / 2.0 + 15.0 * m6 / 32.0 + 7.0 * m8 / 16.0;
 double a4 = m4 / 8.0 + 3.0 * m6 / 16.0 + 7.0 * m8 / 32.0;
 double a6 = m6 / 32.0 + m8 / 16.0;
 double a8 = m8 / 128.0;
 A0 = a0;
 B0 = a2;
 C0 = a4;
 D0 = a6;
 E0 = a8;
 S = (A0 * B - B0* sin(2.0*B)/2.0 + C0 * sin(4.0*B)/4.0 - D0 *sin(6.0*B)/6.0 + E0*sin(8.0*B)/8.0);

 //----高斯投影公式-------
 X = S + N * t *((0.5 + ((5.0 - t * t + 9 * q2 +4 * q2 * q2)/24 + (61.0 - 58.0 * t * t + pow(t,4)) * m * m /720.0) * m * m) * m * m);//pow(t,4)为t的4次方
 Y = N * ((1 + ((1 - t * t +q2)/6.0 +(5.0 - 18.0 * t * t + pow(t,4) + 14 * q2 - 58 * q2 * t * t ) * m * m / 120.0) * m * m ) * m);
 //----根据中国国情坐标纵轴西移500公里当作起始轴
 Y = 1000000 * n + 500000 + Y;

参考网址：http://hi.baidu.com/coldfrost/blog/item/781ad31303b65a29dc5401fc.html

 

三. 面积计算公式

 

注意本公式的首末结点都为（0,0），求之前要把所有的结点都减去第一个结点。

PS：这样求出来的面积是有负值的，因为这是矢量的，和你点的顺序有关系。具体来说顺时针则为负，逆时针则为正。这样又可以多了一种判断地物矢量化方向的方法。

 

四. 结果检测

 

 

 

这是ArcMap中的面积计算结果，我把各个坐标输入进去后得到的结果是654865.33362435829。

误差大约在0.81%。

 

原因分析：

1. 在ArcMap中我选的投影为UTM投影（0.9996），并不是严格的高斯克里格投影。

2.坐标输入的时候，并不是用专门的软件读取出来的，而是手动点击然后copy的。

3.由于简单的结算采用的并不是七参数的转换，所以...

 不过怎么样，证明我的转换算法和面积计算公式都是可行的，哈哈。。。

 

改天把工程和资料打包上传。