GDAL庫學習筆記(1)：無縫拼接Google衛星圖

    開工之前要先了解一下瓦片地圖，瓦片地圖金字塔模型是一種多分辨率層次模型，從瓦片金字塔的底層到頂層，分辨率越來越低，但表示的地理范圍不變。實現原理就是，首先確定地圖服務平台所要提供的縮放級別的數量N，把縮放級別最低、地圖比例尺最大的地圖圖片作為金字塔的底層，即第0層，並對其進行分塊，從地圖圖片的左上角開始，從左至右、從上到下進行切割，分割成相同大小(比如256x256像素)的正方形地圖瓦片，形成第0層瓦片矩陣;在第O層地圖圖片的基礎上，按每2x2像素合成為一個像素的方法生成第1層地圖圖片，並對其進行分塊，分割成與下一層相同大小的正方形地圖瓦片，形成第1層瓦片矩陣;采用同樣的方法生成第2層瓦片矩陣;…;如此下去，直到第N一1層，構成整個瓦片金字塔。

  在Google地圖中，地圖數據由大量的正方形圖片組成。共有2級縮放比例，每個地圖圖片都有坐標值，由X和Y值構成。比例因子zoom取值范圍是(0-22)。操作地圖滑竿顯示更大比例尺地圖時，圖片的數量發生裂變。兩種模式在請求及響應的速度方面有明顯的差異，基於地圖瓦片服務框架的響應速度要快於傳統的WebGIS，同時對地圖服務器的負載也相應小一些。

     現在可以開工，首先准備好瓦片地圖，如果是用來學習研究的話，我們可以通過一些第三方工具（當然自己寫一個也可以）從谷歌地圖抓取瓦片地圖存放到本地，例如以C:\googlemas\satelite\x\y.jpg的格式存放，其中x，y，z分別代瓦片所在的列、行、縮放等級。

 

有了這些圖我們開始拼接了，關於gdal庫的介紹的文章很多了，現在我主要是介紹如何用C#來拼接谷歌地圖，廢話不多說了，直接上代碼：

void SaveBitmapBuffered(Dataset src, Dataset dst, int x, int y)
{
    if (src.RasterCount < 3) {
        System.Environment.Exit(-1);
    }

    // Get the GDAL Band objects from the Dataset
    Band redBand = src.GetRasterBand(1);
    Band greenBand = src.GetRasterBand(2);
    Band blueBand = src.GetRasterBand(3);

    // Get the width and height of the raster
    int width = redBand.XSize;
    int height = redBand.YSize;

    byte[] r = new byte[width * height];
    byte[] g = new byte[width * height];
    byte[] b = new byte[width * height];

    redBand.ReadRaster(0, 0, width, height, r, width, height, 0, 0);
    greenBand.ReadRaster(0, 0, width, height, g, width, height, 0, 0);
    blueBand.ReadRaster(0, 0, width, height, b, width, height, 0, 0);

    Band wrb = dst.GetRasterBand(1);
    wrb.WriteRaster(x * width, y * height, width, height, r, width, height, 0, 0);
    Band wgb = dst.GetRasterBand(2);
    wgb.WriteRaster(x * width, y * height, width, height, g, width, height, 0, 0);
    Band wbb = dst.GetRasterBand(3);
    wbb.WriteRaster(x * width, y * height, width, height, b, width, height, 0, 0);
}

/// <summary>
/// 拼接瓦片
/// </summary>
/// <param name="tilesBounds"></param>
/// <param name="tilePath"></param>
/// <param name="outPutFileName"></param>
public void CombineTiles(TilesBounds tilesBounds, string tilePath, string outPutFileName)
{
    if (File.Exists(outPutFileName))
    {
        File.Delete(outPutFileName);
    }
    int imageWidth = 256 * (tilesBounds.maxCol - tilesBounds.minCol + 1);
    int imageHeight = 256 * (tilesBounds.maxRow - tilesBounds.minRow + 1);

    //Register driver(s). 
    Gdal.AllRegister();
    Driver driver = Gdal.GetDriverByName("GTiff");
    Dataset destDataset = driver.Create(outPutFileName, imageWidth, imageHeight, 3, DataType.GDT_Byte, null);


    for (int col = tilesBounds.minCol; col <= tilesBounds.maxCol; col++)
    {
        for (int row = tilesBounds.minRow; row <= tilesBounds.maxRow; row++)
        {
            try
            {
                string sourceFileName = tilePath + tilesBounds.zoomLevel.ToString() + "\\" + col.ToString() + "\\" + row.ToString() + ".jpg";
                if (File.Exists(sourceFileName))
                {
                    Dataset sourceDataset = Gdal.Open(sourceFileName, Access.GA_ReadOnly);
                    if (sourceDataset != null)
                    {
                        SaveBitmapBuffered(sourceDataset, destDataset, col - tilesBounds.minCol, row - tilesBounds.minRow);
                    }
                }
            }
            catch (Exception ex)
            {
                MessageBox.Show(ex.ToString());
            }
        }
    }
    destDataset.Dispose();

}

//調用
TilesBounds tilesBounds = new TilesBounds();
  tilesBounds.minCol = 107901;
  tilesBounds.maxCol = 107926;
  tilesBounds.minRow = 49652;
  tilesBounds.maxRow = 49668;
  tilesBounds.zoomLevel = 17;
string outPutFileName = "c:\\北京衛星圖" + tilesBounds.zoomLevel.ToString() + ".tif";
string tilePath = "C:\\googlemaps\\\satelite\\";
CombineTiles(tilesBounds, tilePath, outPutFileName);

ok，大功告成了，看看效果圖吧。