關於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只能用迭代法來進行求解了。

// 兲朝火星坐標系偏移公式
// https://on4wp7.codeplex.com/SourceControl/changeset/view/21483#EvilTransform.cs
static double transformLat(double x, double y)
{
	double ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * sqrt(abs(x));
	ret += (20.0 * sin(6.0 * x * M_PI) + 20.0 * sin(2.0 * x * M_PI)) * 2.0 / 3.0;
	ret += (20.0 * sin(y * M_PI) + 40.0 * sin(y / 3.0 * M_PI)) * 2.0 / 3.0;
	ret += (160.0 * sin(y / 12.0 * M_PI) + 320 * sin(y * M_PI / 30.0)) * 2.0 / 3.0;
	return ret;
}

static double transformLon(double x, double y)
{
	double ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * sqrt(abs(x));
	ret += (20.0 * sin(6.0 * x * M_PI) + 20.0 * sin(2.0 * x * M_PI)) * 2.0 / 3.0;
	ret += (20.0 * sin(x * M_PI) + 40.0 * sin(x / 3.0 * M_PI)) * 2.0 / 3.0;
	ret += (150.0 * sin(x / 12.0 * M_PI) + 300.0 * sin(x / 30.0 * M_PI)) * 2.0 / 3.0;
	return ret;
}

// World Geodetic System ==> Mars Geodetic System
void WGS2GCJTransform(double wgLon, double wgLat, double &mgLon, double &mgLat)
{
	const double a = 6378245.0;
	const double ee = 0.00669342162296594323;

	double dLat = transformLat(wgLon - 105.0, wgLat - 35.0);
	double dLon = transformLon(wgLon - 105.0, wgLat - 35.0);

	double radLat = wgLat / 180.0 * M_PI;
	double magic = sin(radLat);
	magic = 1 - ee * magic * magic;

	double sqrtMagic = sqrt(magic);
	dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * M_PI);
	dLon = (dLon * 180.0) / (a / sqrtMagic * cos(radLat) * M_PI);

	mgLat = wgLat + dLat;
	mgLon = wgLon + dLon;
}

void CalcGCJ02ToWGS84Error(const char *pszFile, double dRes)
{
	GDALAllRegister();
	GDALDriver *pDirver = (GDALDriver *)GDALGetDriverByName("GTiff");

	//中國國土面積的四至（陸地區域） 最西為東經 73°，最東為東經 135.5°。最男為北緯 18°，最北為北緯 54°
	double dMinX = 73.0;
	double dMaxX = 135.5;
	double dMinY = 18.0;
	double dMaxY = 54.0;

	// 根據指定的分辨率計算輸出圖像大小
	int nWidth = static_cast<int>((dMaxX - dMinX) / dRes + .5);
	int nHeight = static_cast<int>((dMaxY - dMinY) / dRes + .5);

	printf("%dx%d\n", nWidth, nHeight);

	// 構造輸出圖像的仿射變換參數
	double dGeoTransform[6] = {dMinX, dRes, 0, dMaxY, 0, -dRes};

	// 創建輸出圖像
	GDALDataset* poDS = pDirver->Create(pszFile, nWidth, nHeight, 1, GDT_Float32, NULL);
	poDS->SetGeoTransform(dGeoTransform);
	poDS->SetProjection(SRS_WKT_WGS84);

	//GDALRasterBand *pBandLon = poDS->GetRasterBand(1);
	//GDALRasterBand *pBandLat = poDS->GetRasterBand(2);
	GDALRasterBand *pBandDis = poDS->GetRasterBand(1);
	double *pBuffer = new double[nWidth];
	double *pDstLon = new double[nWidth];
	double *pDstLat = new double[nWidth];

	for (int i=0; i<nHeight; i++)
	{
#pragma omp parallel for
		for (int j=0; j<nWidth; j++)
		{
			double dSrcLon, dSrcLat;
			GDALApplyGeoTransform(dGeoTransform, j, i, &dSrcLon, &dSrcLat);
			WGS2GCJTransform(dSrcLon, dSrcLat, pDstLon[j], pDstLat[j]);

			pDstLon[j] = dSrcLon - pDstLon[j];
			pDstLat[j] = dSrcLat - pDstLat[j];
			double dDis  = sqrt(pDstLon[j]*pDstLon[j] + pDstLat[j]*pDstLat[j]);
			pBuffer[j] = dDis;
		}

		//pBandLon->RasterIO(GF_Write, 0, i, nWidth, 1, pDstLon, nWidth, 1, GDT_Float64, 0, 0);
		//pBandLat->RasterIO(GF_Write, 0, i, nWidth, 1, pDstLat, nWidth, 1, GDT_Float64, 0, 0);
		pBandDis->RasterIO(GF_Write, 0, i, nWidth, 1, pBuffer, nWidth, 1, GDT_Float64, 0, 0);
	}

	RELEASE(pDstLon);
	RELEASE(pDstLat);
	RELEASE(pBuffer);

	GDALClose(GDALDatasetH(poDS));
}