墨卡托投影瓦片转换成WGS84投影瓦片

墨卡托,投影,瓦片,转换成,wgs84 · 浏览次数 : 19

小编点评

**图中两级附近是很色的原因是：** 由于 Mercator 投影无数据完整，所以图中两级附近的地图区域的经度范围在 > 85° 或者 < -85°。这些区域的经度值不在 Mercator 投影范围内，导致图中这些区域无法得到正常渲染。

正文

如题,相信任何一个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;
}