如题,相信任何一个GIS引擎开发者都会遇到的问题，要解决这个问题，首先要了解两者的区别(mercator投影).WGS84投影
如图所示，是一张Mercator投影的地图瓦片，这种投影将地球投影成 W 与 H 相等的图片，是一个正方形(宽度与高度相等)
经纬度方式：-180°，-85°，+180°，+85°, 两极数据超过85°的地方无数据。
下图所示WGS84投影瓦片，图中宽度 = 2倍高度，这种投影比较容易理解
经纬度方式：-180°，-90°，+180°，+90°, 数据无丢失。
目前主流GIS球体都内部都采用WGS84投影（例如 Google Earth,Arg Earth,FastEarth OsgEarth,Cesuim)，这种投影更加容易的被球体适配，所以作为主流使用方式。
二维地图则更多的使用Mercaror投影（例如 Google Map,Baidu,Bing,Arc Map),对于两种投影的数据都在大规模的使用，如何能在三维球体上也能使用Mercator投影的
地图数据，就是本文所要解决的问题。
常见的解决方式有如下几种：
1. 实时CPU计算，用CPU方式实时的将mercator瓦片转换成 wgs84瓦片（代表作: FastEarth)，也是本文的实现方式
优点： 通用性强
缺点: 速度相比GPU要慢一点，因为墨卡托投影 范围 - 180°，+180°， -85°，+ 85°，对于两极数据 没有数据输出，转换后数据丢失。
处理流程如下图所示：
2. 实时GPU计算，利用GPU机型实时计算，该方式速度快（代表作：cesium）。
3. 非实时计算，与第一种类似，提前计算转换便存储计算结果。
4. 参考代码如下：
计算经纬度范围函数,输入的是瓦片的Id，输出瓦片的经纬度范围
/// <summary> /// 获取瓦片编号 /// </summary> aabb2dr getRange(const int3& tileKey) override { unsigned xTiles = 2 << (tileKey.z); unsigned yTiles = 1 << (tileKey.z); double xTileWidth = TWO_PI / xTiles; double yTileHeight = PI / yTiles; double dLon0 = tileKey.x * xTileWidth -PI; double dLat0 = HALF_PI - tileKey.y * yTileHeight; aabb2dr abbb; abbb.setExtents(dLon0,dLat0 - yTileHeight,dLon0 + xTileWidth,dLat0); return abbb; }
根据经纬度范围计算Mercator瓦片Id
aabb2di calcMercatorRange(aabb2dr wgsAabb)
{ wgsAabb._minimum *= RAD2DEG; wgsAabb._maximum *= RAD2DEG; double resX = wgsAabb.getSize().x / tileWidth; double resY = wgsAabb.getSize().y / tileHeight; real2 pxOffset(resX,resY); /// 2. 根据经纬度范围查询mercator id list; int2 tileStart = merctor.getTileId(wgsAabb._minimum + pxOffset,lev); int2 tileEnd = merctor.getTileId(wgsAabb._maximum - pxOffset,lev); int minX = std::min(tileStart.x,tileEnd.x); int minY = std::min(tileStart.y,tileEnd.y); int maxX = std::max(tileStart.x,tileEnd.x); int maxY = std::max(tileStart.y,tileEnd.y);
return aabb2di(minX,minY,maxX,maxY);
}
合成大的瓦片
FEImage toImage(aabb2di box,int tileWidth,int tileHeight) { int width = (box._maximum.x - box._minimum.x + 1) * tileWidth; int height = (box._maximum.y - box._minimum.y + 1) * tileHeight; FEImage srcImage; FEImage tmpImage; srcImage.create(width,height,FEImage::FORMAT_RGB); srcImage.fillColor(Rgba(0,0,0,0)); tmpImage.create(tileWidth,tileHeight,FEImage::FORMAT_RGB); aabb2dr srcBox; /// 3. 加载所有mercator id list /// 4. 合并image 生成一个大的 image; for (int c = minX ; c <= maxX;++ c) { for (int r = minY ; r <= maxY;++ r) { FEObject* ptr = grp->readObject(app,int3(c,r,lev),user,userArg,&tmpImage); if(ptr == nullptr) continue; byte* pBuffer = (byte*)srcImage.dataPtr((c - minX) * tileWidth, (r - minY) * tileHeight); byte* pSrc = (byte*)tmpImage.data(); for (int h = 0; h < tileHeight; ++h) { memcpy(pBuffer,pSrc,tmpImage.pitchSize()); pBuffer += srcImage.pitchSize(); pSrc += tmpImage.pitchSize(); } aabb2dr box = merctor.getRange(int3(c,r,lev)); srcBox.merge(box); } } }
萃取函数：
bool convert( const FEImage& srcImage , FEImage& dstImage , const aabb2dr& srcBox , aabb2dr& dstBox , int dstW =256 , int dstH = 256) { OGRSpatialReference srcSRS; srcSRS.importFromEPSG(3857); OGRSpatialReference dstSRS; dstSRS.importFromEPSG(4326); char* srcPrj = nullptr; char* dstPrj = nullptr; srcSRS.exportToWkt(&srcPrj); dstSRS.exportToWkt(&dstPrj); /// image 数据转换成 gdal数据源 uint w = srcImage.width(); uint h = srcImage.height(); GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("MEM"); GDALDataset* poSrcDS = poDriver->Create( "MEM", w, h,3, GDT_Byte,nullptr); byte* pData = (byte*)srcImage.data(); auto rBand = (GDALRasterBand*)GDALGetRasterBand(poSrcDS, 1 ); auto gBand = (GDALRasterBand*)GDALGetRasterBand(poSrcDS, 2 ); auto bBand = (GDALRasterBand*)GDALGetRasterBand(poSrcDS, 3 ); GDALRasterIOExtraArg exterArg; INIT_RASTERIO_EXTRA_ARG(exterArg); exterArg.eResampleAlg = GDALRIOResampleAlg::GRIORA_Cubic; rBand->RasterIO(GF_Write, 0, 0, w, h, pData + 0, w,h, GDT_Byte, 3, w * 3,&exterArg); gBand->RasterIO(GF_Write, 0, 0, w, h, pData + 1, w,h, GDT_Byte, 3, w * 3,&exterArg); bBand->RasterIO(GF_Write, 0, 0, w, h, pData + 2, w,h, GDT_Byte, 3, w * 3,&exterArg); double srcTrans[6] = {0}; double resX = srcBox.getSize().x / w; double resY = srcBox.getSize().y / h; /// min longitude srcTrans[0] = srcBox._minimum.x; srcTrans[1] = resX; srcTrans[2] = 0; /// max latitude srcTrans[3] = srcBox._maximum.y; srcTrans[4] = 0; srcTrans[5] = -resY; poSrcDS->SetGeoTransform(srcTrans); poSrcDS->SetProjection(srcPrj); double dstTrans[6] = {0}; double dstResX = dstBox.getSize().x / dstW; double dstResY = dstBox.getSize().y / dstH; /// min longitude dstTrans[0] = dstBox._minimum.x; dstTrans[1] = dstResX; dstTrans[2] = 0; /// max latitude dstTrans[3] = dstBox._maximum.y; dstTrans[4] = 0; dstTrans[5] = -dstResY; /// 做转换 GDALDataset* dstDS = extractDataset(poSrcDS,dstTrans,dstW,dstH,dstPrj); if (dstDS == nullptr) { GDALClose(poSrcDS); return false; } byte* pDstBuf = (byte*)dstImage.data(); auto rDst = (GDALRasterBand*)GDALGetRasterBand(dstDS, 1 ); auto gDst = (GDALRasterBand*)GDALGetRasterBand(dstDS, 2 ); auto bDst = (GDALRasterBand*)GDALGetRasterBand(dstDS, 3 ); rDst->RasterIO(GF_Read, 0, 0, dstW, dstH, pDstBuf + 0, dstW, dstH, GDT_Byte, 3, dstW * 3); gDst->RasterIO(GF_Read, 0, 0, dstW, dstH, pDstBuf + 1, dstW, dstH, GDT_Byte, 3, dstW * 3); bDst->RasterIO(GF_Read, 0, 0, dstW, dstH, pDstBuf + 2, dstW, dstH, GDT_Byte, 3, dstW * 3); GDALClose(dstDS); return true; }
结果展示：
图中两级附近是很色的，原因即纬度在 > 85° 或者 < -85°，Mercator投影无数据
完整代码因依赖工程，有需要的小伙伴可以微信联系我,或者留下邮箱
微信号：JingWeiKuaiTu
