// 兲朝火星坐标系偏移公式
// 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));
}