关于GCJ02和WGS84坐标系对比

大家都知道，在兲朝的电子地图的坐标都是经过了一个坐标偏移，叫GCJ_02的东西。在网上发现了将WGS84经纬度转成GCJ02的一个代码，写了个小程序测试了下看看全国各地的偏移量有多大。

关于WGS84转GCJ02的资料网上很多，我参考的是https://on4wp7.codeplex.com/SourceControl/changeset/view/21483#EvilTransform.cs。

先放一个处理的结果图，大概说明一下，绿色为偏移量最小的地方，红色为偏移量最大地方。右下角为图例，最小值为0.000213487度大概为20多米，最大值为0.0104393大概为1公里多（不过这个地点已经超过我国范围了，下图右上角的区域）。关于GCJ02这个坐标系这里不讨论，下面主要说一下怎么用GDAL之类的库生成下面这个彩色的误差图像。

首先，找个中国的四至范围（陆地区域） 最西为东经 73°，最东为东经 135.5°。最男为北纬 18°，最北为北纬 54°，然后指定一个输出图像的格网大小，也就是分辨率，上面这个图大致为10000米也就是10公里一个像素。这样就可以得到这个图像的大小和仿射变换的参数了。

接下来，创建图像，然后遍历图像的每一个像素值，并且计算得到该像素值行列号对应的真实的WGS84经纬度坐标。

然后将WGS84经纬度通过上面的网址里面的转换关系计算转换后的GCJ02坐标系下的经纬度，然后计算这两个经纬度之间的距离，这里简单起见，直接用经纬度的欧拉距离，实际上应该用椭球上的两点大圆距离。

最后将每个点的距离计算出来，写出到图像即可。下面是全部代码。

最后牢骚几句。其实从上面这个图以及下面的公式可以发现，兲超某单位搞出来的这个坐标转换还是很厉害的，误差每个地方都不一样，而且还是连续的，从公式中可以看出，坐标转换的公式由一个关于经纬度的线性多项式（次数从0.5到2）加上经纬度的正弦函数组成。如果都是线性多项式的话，可以很容易推到反函数，但是后面加了一个非线性的函数（正弦函数应该是为了周期性的增加误差用的），这样反函数就非常不容易推导出来。所以关于从GCJ02坐标系转到WGS84只能用迭代法来进行求解了。

[cpp] view plaincopyprint?

1. // 兲朝火星坐标系偏移公式  
2. // https://on4wp7.codeplex.com/SourceControl/changeset/view/21483#EvilTransform.cs  
3. static double transformLat(double x, double y)  
4. {  
5.     double ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * sqrt(abs(x));  
6.     ret += (20.0 * sin(6.0 * x * M_PI) + 20.0 * sin(2.0 * x * M_PI)) * 2.0 / 3.0;  
7.     ret += (20.0 * sin(y * M_PI) + 40.0 * sin(y / 3.0 * M_PI)) * 2.0 / 3.0;  
8.     ret += (160.0 * sin(y / 12.0 * M_PI) + 320 * sin(y * M_PI / 30.0)) * 2.0 / 3.0;  
9.     return ret;  
10. }  
11.   
12. static double transformLon(double x, double y)  
13. {  
14.     double ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * sqrt(abs(x));  
15.     ret += (20.0 * sin(6.0 * x * M_PI) + 20.0 * sin(2.0 * x * M_PI)) * 2.0 / 3.0;  
16.     ret += (20.0 * sin(x * M_PI) + 40.0 * sin(x / 3.0 * M_PI)) * 2.0 / 3.0;  
17.     ret += (150.0 * sin(x / 12.0 * M_PI) + 300.0 * sin(x / 30.0 * M_PI)) * 2.0 / 3.0;  
18.     return ret;  
19. }  
20.   
21. // World Geodetic System ==> Mars Geodetic System  
22. void WGS2GCJTransform(double wgLon, double wgLat, double &mgLon, double &mgLat)  
23. {  
24.     const double a = 6378245.0;  
25.     const double ee = 0.00669342162296594323;  
26.   
27.     double dLat = transformLat(wgLon - 105.0, wgLat - 35.0);  
28.     double dLon = transformLon(wgLon - 105.0, wgLat - 35.0);  
29.   
30.     double radLat = wgLat / 180.0 * M_PI;  
31.     double magic = sin(radLat);  
32.     magic = 1 - ee * magic * magic;  
33.   
34.     double sqrtMagic = sqrt(magic);  
35.     dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * M_PI);  
36.     dLon = (dLon * 180.0) / (a / sqrtMagic * cos(radLat) * M_PI);  
37.   
38.     mgLat = wgLat + dLat;  
39.     mgLon = wgLon + dLon;  
40. }  
41.   
42. void CalcGCJ02ToWGS84Error(const char *pszFile, double dRes)  
43. {  
44.     GDALAllRegister();  
45.     GDALDriver *pDirver = (GDALDriver *)GDALGetDriverByName("GTiff");  
46.   
47.     //中国国土面积的四至（陆地区域） 最西为东经 73°，最东为东经 135.5°。最男为北纬 18°，最北为北纬 54°  
48.     double dMinX = 73.0;  
49.     double dMaxX = 135.5;  
50.     double dMinY = 18.0;  
51.     double dMaxY = 54.0;  
52.   
53.     // 根据指定的分辨率计算输出图像大小  
54.     int nWidth = static_cast<int>((dMaxX - dMinX) / dRes + .5);  
55.     int nHeight = static_cast<int>((dMaxY - dMinY) / dRes + .5);  
56.   
57.     printf("%dx%d\n", nWidth, nHeight);  
58.   
59.     // 构造输出图像的仿射变换参数  
60.     double dGeoTransform[6] = {dMinX, dRes, 0, dMaxY, 0, -dRes};  
61.   
62.     // 创建输出图像  
63.     GDALDataset* poDS = pDirver->Create(pszFile, nWidth, nHeight, 1, GDT_Float32, NULL);  
64.     poDS->SetGeoTransform(dGeoTransform);  
65.     poDS->SetProjection(SRS_WKT_WGS84);  
66.   
67.     //GDALRasterBand *pBandLon = poDS->GetRasterBand(1);  
68.     //GDALRasterBand *pBandLat = poDS->GetRasterBand(2);  
69.     GDALRasterBand *pBandDis = poDS->GetRasterBand(1);  
70.     double *pBuffer = new double[nWidth];  
71.     double *pDstLon = new double[nWidth];  
72.     double *pDstLat = new double[nWidth];  
73.   
74.     for (int i=0; i<nHeight; i++)  
75.     {  
76. #pragma omp parallel for  
77.         for (int j=0; j<nWidth; j++)  
78.         {  
79.             double dSrcLon, dSrcLat;  
80.             GDALApplyGeoTransform(dGeoTransform, j, i, &dSrcLon, &dSrcLat);  
81.             WGS2GCJTransform(dSrcLon, dSrcLat, pDstLon[j], pDstLat[j]);  
82.   
83.             pDstLon[j] = dSrcLon - pDstLon[j];  
84.             pDstLat[j] = dSrcLat - pDstLat[j];  
85.             double dDis  = sqrt(pDstLon[j]*pDstLon[j] + pDstLat[j]*pDstLat[j]);  
86.             pBuffer[j] = dDis;  
87.         }  
88.   
89.         //pBandLon->RasterIO(GF_Write, 0, i, nWidth, 1, pDstLon, nWidth, 1, GDT_Float64, 0, 0);  
90.         //pBandLat->RasterIO(GF_Write, 0, i, nWidth, 1, pDstLat, nWidth, 1, GDT_Float64, 0, 0);  
91.         pBandDis->RasterIO(GF_Write, 0, i, nWidth, 1, pBuffer, nWidth, 1, GDT_Float64, 0, 0);  
92.     }  
93.   
94.     RELEASE(pDstLon);  
95.     RELEASE(pDstLat);  
96.     RELEASE(pBuffer);  
97.   
98.     GDALClose(GDALDatasetH(poDS));  
99. }  

转载自 http://blog.csdn.net/liminlu0314/article/details/42564039