關於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)); }