GDAL-影像按照波段进行切片（c#）

       public void SplitImagery(string pFilePath)
        {
            Dataset pDataset = Gdal.Open(pFilePath, Access.GA_ReadOnly);
            double[] pGeoTransform = new double[6];
            pDataset.GetGeoTransform(pGeoTransform);
            string pProjection = pDataset.GetProjectionRef();
            Driver pDriver = pDataset.GetDriver();
            int imgSizeX = pDataset.RasterXSize;
            int imgSizeY = pDataset.RasterYSize;
            int iBandCount = pDataset.RasterCount;
            string pDirectoryPath = System.IO.Path.GetDirectoryName(pFilePath);
            string pFileName = System.IO.Path.GetFileNameWithoutExtension(pFilePath);
            string pFileExtentsion = System.IO.Path.GetExtension(pFilePath);
            for (int iBand = 0; iBand < iBandCount; iBand++)
            {
                Band pBand = pDataset.GetRasterBand(iBand+1);
                const int block_size = 2048;//裁剪块的大小
                int rowChunkIndex = 0;
                for (int yIndex = 0; yIndex < imgSizeY; yIndex += block_size)
                {
                    rowChunkIndex++;
                    int columnChunkIndex = 0;
                    for (int xIndex = 0; xIndex < imgSizeX; xIndex += block_size)
                    {
                        columnChunkIndex++;
                        string pSplitedFilePath = pDirectoryPath + "\\" + pFileName+"_"+ iBand + "_"+ rowChunkIndex+"_"+ columnChunkIndex + pFileExtentsion;
                        //定义两个变量来保存分块大小
                        int nXBK = block_size;//列
                        int nYBK = block_size;//行
                        //如果最下面和最右边的块不够256，剩下多少读取多少
                        if (yIndex + block_size > imgSizeY)     //最下面的剩余块
                            nYBK = imgSizeY - yIndex;
                        if (xIndex + block_size > imgSizeX)     //最右侧的剩余块
                            nXBK = imgSizeX - xIndex;
                        //将图像数据集读到pBuf指针中
                        double[] buffer = new double[nXBK * nYBK];
                        pBand.ReadRaster(xIndex, yIndex, nXBK, nYBK, buffer, nXBK, nYBK, 0, 0);
                        double[] newBuffer = new double[nXBK * nYBK];
                        for (int i = 0; i < nXBK; i++)
                        {
                            for (int j = 0; j < nYBK; j++)
                            {
                                double val =  buffer[j * nXBK + i];

                                newBuffer[j * nXBK + i] = val;//分块读入计算，分块写入
                            }
                        }
                        //定义DataType为64位会产生成倍的数据冗余，但是调试安全
                        Dataset pSpliteDataset = pDriver.Create(pSplitedFilePath, nXBK, nYBK, 1, DataType.GDT_Float64, null);
                        //重新定义裁剪的影像的仿射变换
                        double[] pSpliteGeoTransform = new double[6];
                        pSpliteGeoTransform[0] = pGeoTransform[0] + pGeoTransform[1] * xIndex + pGeoTransform[2] * yIndex;
                        pSpliteGeoTransform[1] = pGeoTransform[1];
                        pSpliteGeoTransform[2] = pGeoTransform[2];
                        pSpliteGeoTransform[3] = pGeoTransform[3] + pGeoTransform[4] * xIndex + pGeoTransform[5] * yIndex;
                        pSpliteGeoTransform[4] = pGeoTransform[4];
                        pSpliteGeoTransform[5] = pGeoTransform[5];
                        pSpliteDataset.SetGeoTransform(pSpliteGeoTransform);
                        pSpliteDataset.SetProjection(pProjection);
                        pSpliteDataset.WriteRaster(0, 0, nXBK, nYBK, newBuffer, nXBK, nYBK, 1, new int[] { 1 }, 0, 0, 0);

                        buffer = null;
                        newBuffer = null;
                        pSpliteDataset.FlushCache();
                        GC.Collect();
                    }
                }
            }
        }