C++调用GDAL库读取并输出tif文件，并计算斑块面积输出景观指数：ＣＳＤ

部分源码选自GDAL库的官方网址：www.gdal.org，其余的代码为笔者自己编写。

// readfile.cpp : 定义控制台应用程序的入口点。
//
/*
part of the codes were cite from:    http://www.gdal.org/gdal_tutorial.html
and remaining of code were created :by www.cnblogs.com/AmatVictorialCuram/‎
and please mark the site if you cite the program fo commercial of publication.
*/
#include "stdafx.h"
#include <iostream>
#include "gdal.h"
#include "gdal_priv.h"
/**/
#include <iomanip>
#include <fstream>
using namespace std;
/**/
int minLabel(int a,int b);

int maxLabel(int a,int b);

double csd(int *Parea,int length);

int main()
{   
    /*
    part of the codes were cite from:    http://www.gdal.org/gdal_tutorial.html
    and remaining of code were created :by www.cnblogs.com/AmatVictorialCuram/‎
    and please mark the site if you cite the program fo commercial of publication.
    */
    const char *pszFilename="E:\tif/fragstats/sample.tif";
    GDALDataset  *poDataset;
    
    GDALAllRegister();

    poDataset = (GDALDataset *) GDALOpen( pszFilename, GA_ReadOnly );
    if( poDataset == NULL )
    {
        cout<<"file read error!"<<endl;;
    }
    double        adfGeoTransform[6];

    printf( "Driver: %s/%s
",
            poDataset->GetDriver()->GetDescription(), 
            poDataset->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME ) );

    printf( "Size is %dx%dx%d
", 
            poDataset->GetRasterXSize(), poDataset->GetRasterYSize(),
            poDataset->GetRasterCount() );

    if( poDataset->GetProjectionRef()  != NULL )
        printf( "Projection is `%s'
", poDataset->GetProjectionRef() );

    if( poDataset->GetGeoTransform( adfGeoTransform ) == CE_None )
    {
        printf( "Origin = (%.6f,%.6f)
",
                adfGeoTransform[0], adfGeoTransform[3] );

        printf( "Pixel Size = (%.6f,%.6f)
",
                adfGeoTransform[1], adfGeoTransform[5] );
    }

    GDALRasterBand  *poBand;
    int             nBlockXSize, nBlockYSize;
    int             bGotMin, bGotMax;
    double          adfMinMax[2];

    poBand = poDataset->GetRasterBand( 1 );
    poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
    printf( "Block=%dx%d Type=%s, ColorInterp=%s
",
        nBlockXSize, nBlockYSize,
        GDALGetDataTypeName(poBand->GetRasterDataType()),
        GDALGetColorInterpretationName(
        poBand->GetColorInterpretation()) );

    adfMinMax[0] = poBand->GetMinimum( &bGotMin );
    adfMinMax[1] = poBand->GetMaximum( &bGotMax );
    if( ! (bGotMin && bGotMax) )
        GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);

    printf( "Min=%.3fd, Max=%.3f
", adfMinMax[0], adfMinMax[1] );

    if( poBand->GetOverviewCount() > 0 )
        printf( "Band has %d overviews.
", poBand->GetOverviewCount() );

    if( poBand->GetColorTable() != NULL )
        printf( "Band has a color table with %d entries.
", 
        poBand->GetColorTable()->GetColorEntryCount() );
    /*****
    Reading Raster Data
    *****/
    float *pafScanline;
    int nXSize = poBand->GetXSize();
    int nYSize=poBand->GetYSize();

    //读取图像的nXSize*nYSize数据
    pafScanline = (float*) CPLMalloc(sizeof(float)*nXSize*nYSize);//创建指针
    poBand->RasterIO(
        GF_Read,//第一个参数表示要读入数据还是写入数据
        0, 0,//nXOff, nYOff表示读取或者写入图像数据的起始坐标图像的左上角坐标为（0,0）
        nXSize, nYSize,/*nXSize, nYSize表示读取或者写入图像数据的窗口大小,nXSize表示宽度,
                       nYSize表示高度,均使用像素为单位,该宽度和高度是从第二个和第三个参数处开始计算。
                       这两个参数和第二第三个参数一起表示就是,读取和写入图像的窗口位置和大小。*/
        pafScanline, //指向存储数据的一个指针
        nXSize, nYSize,//指定缓冲区的大小
        GDT_Float32, 
        0, 0 );
    // poBand->RasterIO( GF_Read, 0, 0, nXSize, nYSize, pafScanline, nXSize, nYSize, GDT_Float32, 0, 0 );
    cout<<"列数：nXSize="<<nXSize<<endl;    
    cout<<"行数：nYSize="<<nYSize<<endl;
    //system("pause");



    //创建nXSize*nYSize的float数组
    float **data=new float *[nXSize];//每列有nXSize列,data数组的大小与dataLabel数组的大小相同,类型不同
    int **dataLabel=new int *[nXSize];//创建标签数组,有nYSize行,nXSize列
    int t=0,i=0,j=0;
    for(int i=0;i<nYSize;i++)//nYSize行
    {
        dataLabel[i]=new int [nXSize];
        data[i]     = new float [nXSize];
    }

    cout<<"输出元数据数组data："<<endl;
    for(int i=0;i<nYSize;i++)
    {
        for(int j=0;j<nXSize;j++)
        {
            data[i][j]=pafScanline[t++];
            cout<<setw(4)<<data[i][j];
        }
        cout<<endl;
    }
    cout<<endl;
    //system("pause");
    cout<<"输出标签数组dataLabel数组："<<endl;
    for(int i=0;i<nYSize;i++)
    {
        cout<<endl;
        for(int j=0;j<nXSize;j++)
        {
            dataLabel[i][j]=nYSize*nXSize;
            cout<<setw(4)<<dataLabel[i][j];
        }
    }
    //system("pause");
    cout<<endl;
    cout<<"联通区域标记算法开始："<<endl;
        dataLabel[0][0]=1;//把左上角的第一个标签赋值为1
    int indexcolor=1;
    for(i=0;i<nYSize;i++)
    {
        for(j=0;j<nXSize;j++)
        {
            if(i==0)
            {
                if(j==0)
                    continue;
                if(data[i][j]==data[i][j-1])
                    dataLabel[i][j]=dataLabel[i][j-1];
                else
                {
                    dataLabel[i][j]=++indexcolor;
                }
                continue;
            }
            if(j==0)
            {
                if(data[i][j]==data[i-1][j])
                {
                    if(dataLabel[i-1][j]<dataLabel[i][j])
                        dataLabel[i][j]=dataLabel[i-1][j];
                }
                if(data[i][j]==data[i-1][j+1])
                {
                    if(dataLabel[i-1][j+1]<dataLabel[i][j])
                        dataLabel[i][j]=dataLabel[i-1][j+1];
                }
                if(dataLabel[i][j]==nYSize*nXSize)
                {
                    dataLabel[i][j]=++indexcolor;
                }
                continue;
            }
            //if(edgeCheckL(i,j) && data[i][j]==data[i][j-1])//左边的,
            if(j>0 && data[i][j]==data[i][j-1])
            {
                if(dataLabel[i][j-1]<dataLabel[i][j])
                    dataLabel[i][j]=dataLabel[i][j-1];


            }
            //if(edgeCheckLU(i,j) && data[i][j]==data[i-1][j-1])//左上角的
            if((i>0 && j>0)&& data[i][j]==data[i-1][j-1])
            {
                if(dataLabel[i-1][j-1]<dataLabel[i][j])
                    dataLabel[i][j]=dataLabel[i-1][j-1];
            }
            //if(edgeCheckU(i,j) && data[i][j]==data[i-1][j])//上面的
            if(i>0 && data[i][j]==data[i-1][j])
            {
                if(dataLabel[i-1][j]<dataLabel[i][j])
                    dataLabel[i][j]=dataLabel[i-1][j];
            }
            //if(edgeCheckUR(i,j) && data[i][j]==data[i-1][j+1])//右上角的
            if((i>0 && j<nYSize-1) && data[i][j]==data[i-1][j+1])
            {
                if(dataLabel[i-1][j+1]<dataLabel[i][j])
                    dataLabel[i][j]=dataLabel[i-1][j+1];
            }
            if(dataLabel[i][j]==nYSize*nXSize)
                dataLabel[i][j]=++indexcolor;
        }
    }
    //system("pause");
    cout<<endl;
    cout<<"输出第一次标记数组"<<endl<<endl;
    for(i=0;i<nYSize;i++)
    {
        for(j=0;j<nXSize;j++)
        {
            cout<<setw(3)<<dataLabel[i][j];
        }
        cout<<endl;
    }
    //system("pause");
    cout<<"合并首次生成的标签数组。。。。"<<endl<<endl;

    //合并：像素值相同，但是标签不同，就把标签值大的变为小的。
    for(i=0;i<nYSize;i++)//行
    {
        for(j=0;j<nXSize;j++)//列
        {
            if(i==0)//第0行,只判左边的！
            {
                if(j==0)//第一个元素
                {
                    j=1;//跳过第一个元素，直接从第二个元素：data[0][1]判断
                    //判断并执行合并：
                    if(data[i][j-1]==data[i][j] && dataLabel[i][j-1]!=dataLabel[i][j])//像素值相同  && 标签不同：合并！
                    {
                        //把所有的大标签替换为当前两个标签中的一个较小的值
                        //执行一次遍历
                        //cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
                        for(int m=0;m<nYSize;m++)
                        {
                            for(int n=0;n<nXSize;n++)
                            {//把大的标签替换为小的标签
                                if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i][j-1]))
                                {
                                    dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i][j-1]);
                                }
                            }
                        }
                    }
                }
            }
            else//非0行
            {
                if(j==0)//第0列，但不是第一行的：只判断上、右上两个方向
                {
                    if(data[i-1][j]==data[i][j] && dataLabel[i-1][j]!=dataLabel[i][j])//像素值相同  && 标签不同：合并！
                    {
                        //把所有的大标签替换为当前两个标签中的一个较小的值
                        //执行一次遍历
                        //cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
                        for(int m=0;m<nYSize;m++)
                        {
                            for(int n=0;n<nXSize;n++)
                            {//把大的标签替换为小的标签
                                if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i-1][j]))
                                {
                                    dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i-1][j]);
                                }
                            }
                        }
                    }
                    if(data[i-1][j+1]==data[i][j] && dataLabel[i-1][j=1]!=dataLabel[i][j])//像素值相同  && 标签不同：合并！
                    {
                        //把所有的大标签替换为当前两个标签中的一个较小的值
                        //执行一次遍历
                        //cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
                        for(int m=0;m<nYSize;m++)
                        {
                            for(int n=0;n<nXSize;n++)
                            {//把大的标签替换为小的标签
                                if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i-1][j+1]))
                                {
                                    dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i-1][j+1]);
                                }
                            }
                        }
                    }
                }
                else if(j==nYSize-1)//非0行且最后一列的：判断左、左上、上三个方向
                {
                    if(data[i][j-1]==data[i][j] && dataLabel[i][j-1]!=dataLabel[i][j])//像素值相同  && 标签不同：合并！
                    {
                        //把所有的大标签替换为当前两个标签中的一个较小的值
                        //执行一次遍历
                        //cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
                        for(int m=0;m<nYSize;m++)
                        {
                            for(int n=0;n<nXSize;n++)
                            {//把大的标签替换为小的标签
                                if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i][j-1]))
                                {
                                    dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i][j-1]);
                                }
                            }
                        }
                    }
                    if(data[i-1][j-1]==data[i][j] && dataLabel[i-1][j-1]!=dataLabel[i][j])//像素值相同  && 标签不同：合并！
                    {
                        //把所有的大标签替换为当前两个标签中的一个较小的值
                        //执行一次遍历
                        //cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
                        for(int m=0;m<nYSize;m++)
                        {
                            for(int n=0;n<nXSize;n++)
                            {//把大的标签替换为小的标签
                                if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i-1][-1]))
                                {
                                    dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i-1][j-1]);
                                }
                            }
                        }
                    }
                    if(data[i-1][j]==data[i][j] && dataLabel[i-1][j]!=dataLabel[i][j])//像素值相同  && 标签不同：合并！
                    {
                        //把所有的大标签替换为当前两个标签中的一个较小的值
                        //执行一次遍历
                        //cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
                        for(int m=0;m<nYSize;m++)
                        {
                            for(int n=0;n<nXSize;n++)
                            {//把大的标签替换为小的标签
                                if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i-1][j]))
                                {
                                    dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i-1][j]);
                                }
                            }
                        }
                    }

                }
                else//非0行且（既不是第一列，也不是最后一列的）：判断左、左上、上、右上四个方向
                {
                    if(data[i][j-1]==data[i][j] && dataLabel[i][j-1]!=dataLabel[i][j])//像素值相同  && 标签不同：合并！
                    {
                        //把所有的大标签替换为当前两个标签中的一个较小的值
                        //执行一次遍历
                        //cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
                        for(int m=0;m<nYSize;m++)
                        {
                            for(int n=0;n<nXSize;n++)
                            {//把大的标签替换为小的标签
                                if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i][j-1]))
                                {
                                    dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i][j-1]);
                                }
                            }
                        }
                    }
                    if(data[i-1][j-1]==data[i][j] && dataLabel[i-1][j-1]!=dataLabel[i][j])//像素值相同  && 标签不同：合并！
                    {
                        //把所有的大标签替换为当前两个标签中的一个较小的值
                        //执行一次遍历
                        //cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
                        for(int m=0;m<nYSize;m++)
                        {
                            for(int n=0;n<nXSize;n++)
                            {//把大的标签替换为小的标签
                                if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i-1][-1]))
                                {
                                    dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i-1][j-1]);
                                }
                            }
                        }
                    }
                    if(data[i-1][j]==data[i][j] && dataLabel[i-1][j]!=dataLabel[i][j])//像素值相同  && 标签不同：合并！
                    {
                        //把所有的大标签替换为当前两个标签中的一个较小的值
                        //执行一次遍历
                        //cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
                        for(int m=0;m<nYSize;m++)
                        {
                            for(int n=0;n<nXSize;n++)
                            {//把大的标签替换为小的标签
                                if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i-1][j]))
                                {
                                    dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i-1][j]);
                                }
                            }
                        }
                    }
                    if(data[i-1][j+1]==data[i][j] && dataLabel[i-1][j+1]!=dataLabel[i][j])//像素值相同  && 标签不同：合并！
                    {
                        //把所有的大标签替换为当前两个标签中的一个较小的值
                        //执行一次遍历
                        //cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
                        for(int m=0;m<nYSize;m++)
                        {
                            for(int n=0;n<nXSize;n++)
                            {//把大的标签替换为小的标签
                                if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i-1][j+1]))
                                {
                                    dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i-1][j+1]);
                                }
                            }
                        }
                    }
                }
            }

        }
    }
    t=1;
    //system("pause");
    cout<<"输出合并后的标签数组dataLabel:"<<endl;
    for(i=0;i<nYSize;i++)
    {
        for(j=0;j<nXSize;j++)
        {
            if(dataLabel[i][j]>t)
            {
                t++;
            }
            cout<<setw(3)<<dataLabel[i][j];

        }
        cout<<endl;
    }
    cout<<endl<<"t="<<t<<endl;
    cout<<"联通区域标记算法结束end！！！"<<endl;
    //system("pause");
    int max=0;
    for(int i=0;i<nYSize;i++)
    {
        for(int j=0;j<nXSize;j++)
        {
            if(dataLabel[i][j]>max)
            {
                max=dataLabel[i][j];
            }
        }
    }
    cout<<"maxLabel="<<max<<endl;
    int *Pid=new int [max+1];
    int *Parea=new int [max+1];
    for(i=0;i<max+1;i++)
    {
        Pid[i]=i;
        Parea[i]=0;
    }
    /////////////////////////////
    for(int i=0;i<nYSize;i++)
    {
        for(int j=0;j<nXSize;j++)
        {
                Parea[dataLabel[i][j]]++;
        }
    }
    t=max+1;
    for(i=1;i<t;i++)
    {
        
        cout<<"dataLabel为"<<i<<"的面积为:"<<Parea[i]<<endl;
    }
    
    double Xi=0,Si=0;
    int NPatch=0;
    for(i=0;i<t;i++)
    {
        if(Parea[i]!=0)
        {
            NPatch++;
            Xi=Xi+Parea[i];//求出板块总面积
        }
    /*cout<<"NPatch="<<NPatch<<endl;
    cout<<"Xi="<<Xi<<endl;*/
    }
    cout<<"NPatch="<<NPatch<<endl;
    //cout<<"Xi="<<Xi<<endl;
    Xi=Xi/NPatch;//面积平均数
    cout<<"Initial Value:Xi="<<Xi<<endl;
        cout<<"Initial Value:Si="<<Si<<endl;
        //计算出所有
    for(i=1;i<t;i++)
    {
        if(Parea[i]!=0)
        {
            Si=Si+((Parea[i]-Xi)*(Parea[i]-Xi))/NPatch;
        }
    }
    cout<<"方差="<<Si<<endl;
    Si=sqrt(Si);
    cout<<"Standard Deviation(标准差)="<<Si<<endl;
    cout<<endl;

    /****输出内容写出到文件当中****/
    ofstream ocout;
    //ocout.open("result.csv");
    //ocout.open("result.txt");
    ocout.open("result.csv");
    /****将计算出的结果输出到屏幕****/
    ocout<<"PatchID"<<","<<"Area"<<","<<"CSD"<<endl;
    for(i=1;i<t;i++)
    {
        if(Parea[i]!=0)
        {
            
            //ocout<<"Patch: Label="<<i<<setw(4)<<"  Area="<<Parea[i]<<setw(10)<<"CSD="<<(Parea[i]-Xi)/Si<<endl;
            ocout<<i<<",";
            ocout<<Parea[i]<<",";
            ocout<<(Parea[i]-Xi)/Si<<","<<endl;
        }
    }
    /*******后面的这一部分没有太大的实质性用处，可有可无**********/
     const char *pszFormat = "GTiff";
    GDALDriver *poDriver;
    char **papszMetadata;

    poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);

    if( poDriver == NULL )
        exit( 1 );

    papszMetadata = poDriver->GetMetadata();
    if( CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATE, FALSE ) )
        printf( "Driver %s supports Create() method.
", pszFormat );
    if( CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATECOPY, FALSE ) )
        printf( "Driver %s supports CreateCopy() method.
", pszFormat );

    /*释放读取文件用的指针*/
    CPLFree(pafScanline);
    //关闭文件
    GDALClose((GDALDatasetH)poDataset);
    free(dataLabel);
    free(data);
    free(Parea);
    free(Pid);
    //
    ocout.close();
}
int minLabel(int a,int b)
{
    return (a>b)?b:a;
}
int maxLabel(int a,int b)
{
    return (a<b)?b:a;
}
double csd(int *Parea,int length)
{
    return 0;
}