十五、landsat水体提取部分c++代码

GDALAllRegister();
    CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
    GDALDataset* poDataset = (GDALDataset*)GDALOpen(opt_water_filename1.toLocal8Bit().data(), GA_ReadOnly);
    GDALRasterBand *raseterBandGreen, *raseterBandNir;
    unsigned short* bufferBlockGreen = (unsigned short*)CPLMalloc(sizeof(unsigned short)* block_size*block_size);
    unsigned short* bufferBlockNIR = (unsigned short*)CPLMalloc(sizeof(unsigned short)* block_size*block_size);
    QStringList sections = opt_water_filename1.split(".");
    QString processfilename = sections.at(0).trimmed().append("temp.tif");
    opt_water_filename2 = processfilename;

    GDALDataset *poDataset2 = GetGDALDriverManager()->GetDriverByName(pszFormat)->Create(processfilename.toLocal8Bit().data(), poDataset->GetRasterXSize(), poDataset->GetRasterYSize(), 1, GDT_Float64, NULL);
    GDALRasterBand * processRasterBand = poDataset2->GetRasterBand(1);
    double* processBufferBlock = (double*)CPLMalloc(sizeof(double)*block_size*block_size);

    poDataset->GetGeoTransform(padfTransform);
    
    if (water_exaction == 1 || water_exaction == 2){
        if (water_exaction == 1){
            raseterBandGreen = poDataset->GetRasterBand(3);
            raseterBandNir = poDataset->GetRasterBand(5);
        }
        else{
            raseterBandGreen = poDataset->GetRasterBand(3);
            raseterBandNir = poDataset->GetRasterBand(6);
        }
        for (int i = 0; i < poDataset->GetRasterYSize(); i += block_size)
        {
            for (int j = 0; j < poDataset->GetRasterXSize(); j += block_size)

            {
                int nXBK = block_size;
                int nYBK = block_size;
                if (i + block_size > poDataset->GetRasterYSize())
                    nYBK = poDataset->GetRasterYSize() - i;
                if (j + block_size > poDataset->GetRasterXSize())
                    nXBK = poDataset->GetRasterXSize() - j;
                raseterBandGreen->RasterIO(GF_Read, j, i, nXBK, nYBK, bufferBlockGreen, nXBK, nYBK, GDT_UInt16, 0, 0);
                raseterBandNir->RasterIO(GF_Read, j, i, nXBK, nYBK, bufferBlockNIR, nXBK, nYBK, GDT_UInt16, 0, 0);
                for (int n = 0; n < nYBK*nXBK; n++){
                    processBufferBlock[n] = 1.0*(bufferBlockGreen[n] - bufferBlockNIR[n]) / (bufferBlockGreen[n] + bufferBlockNIR[n]);
                }
                processRasterBand->RasterIO(GF_Write, j, i, nXBK, nYBK, processBufferBlock, nXBK, nYBK, GDT_Float64, 0, 0);
            }
        }
    }
    else if (water_exaction == 3){
        raseterBandGreen = poDataset->GetRasterBand(3);
        raseterBandNir = poDataset->GetRasterBand(5);
        GDALRasterBand *raseterBandSWI = poDataset->GetRasterBand(6);
        unsigned short* bufferBlockSWI = (unsigned short*)CPLMalloc(sizeof(unsigned short)* block_size*block_size);
        for (int i = 0; i < poDataset->GetRasterYSize(); i += block_size)
        {
            for (int j = 0; j < poDataset->GetRasterXSize(); j += block_size)
            {
                int nXBK = block_size;
                int nYBK = block_size;
                if (i + block_size > poDataset->GetRasterYSize())
                    nYBK = poDataset->GetRasterYSize() - i;
                if (j + block_size > poDataset->GetRasterXSize())
                    nXBK = poDataset->GetRasterXSize() - j;
                raseterBandGreen->RasterIO(GF_Read, j, i, nXBK, nYBK, bufferBlockGreen, nXBK, nYBK, GDT_UInt16, 0, 0);
                raseterBandNir->RasterIO(GF_Read, j, i, nXBK, nYBK, bufferBlockNIR, nXBK, nYBK, GDT_UInt16, 0, 0);
                for (int n = 0; n < nYBK*nXBK; n++){
                    processBufferBlock[n] = 1.0*(2 *bufferBlockGreen[n] - bufferBlockNIR[n] - bufferBlockSWI[n]) / ( 2*bufferBlockGreen[n] + bufferBlockNIR[n] + bufferBlockSWI[n]);
                }
                processRasterBand->RasterIO(GF_Write, j, i, nXBK, nYBK, processBufferBlock, nXBK, nYBK, GDT_Float64, 0, 0);
            }
        }
    }
    Pro = poDataset->GetProjectionRef();
    nodata = poDataset->GetRasterBand(1)->GetNoDataValue();
    poDataset2->GetRasterBand(1)->SetNoDataValue(nodata);
    poDataset2->SetGeoTransform(padfTransform);
    poDataset2->SetProjection(Pro);

    double dMinMax[2] = { -1.0, 1.0 };
    processRasterBand->ComputeRasterMinMax(FALSE, dMinMax);
    GUIntBig panHistogram[1500] = { 0 };
    processRasterBand->GetHistogram(dMinMax[0] - 0.5, dMinMax[1] + 0.5, 1500, panHistogram, TRUE, FALSE, NULL, NULL);

    QImage img = Paint(imageHist, panHistogram, dMinMax[0], dMinMax[1]);//开始绘图
    CPLFree(bufferBlockGreen);
    CPLFree(bufferBlockNIR);
    GDALClose(poDataset);
    CPLFree(processBufferBlock);
    GDALClose(poDataset2);