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

部分源码选自GDAL库的官方网址：www.gdal.org，其余的代码为笔者自己编写。
  1 // readfile.cpp : 定义控制台应用程序的入口点。
  2 //
  3 /*
  4 part of the codes were cite from:    http://www.gdal.org/gdal_tutorial.html
  5 and remaining of code were created :by www.cnblogs.com/AmatVictorialCuram/‎
  6 and please mark the site if you cite the program fo commercial of publication.
  7 */
  8 #include "stdafx.h"
  9 #include <iostream>
 10 #include "gdal.h"
 11 #include "gdal_priv.h"
 12 /**/
 13 #include <iomanip>
 14 #include <fstream>
 15 using namespace std;
 16 /**/
 17 int minLabel(int a,int b);
 18 
 19 int maxLabel(int a,int b);
 20 
 21 double csd(int *Parea,int length);
 22 
 23 int main()
 24 {   
 25     /*
 26     part of the codes were cite from:    http://www.gdal.org/gdal_tutorial.html
 27     and remaining of code were created :by www.cnblogs.com/AmatVictorialCuram/‎
 28     and please mark the site if you cite the program fo commercial of publication.
 29     */
 30     const char *pszFilename="E:\tif/fragstats/sample.tif";
 31     GDALDataset  *poDataset;
 32     
 33     GDALAllRegister();
 34 
 35     poDataset = (GDALDataset *) GDALOpen( pszFilename, GA_ReadOnly );
 36     if( poDataset == NULL )
 37     {
 38         cout<<"file read error!"<<endl;;
 39     }
 40     double        adfGeoTransform[6];
 41 
 42     printf( "Driver: %s/%s
",
 43             poDataset->GetDriver()->GetDescription(), 
 44             poDataset->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME ) );
 45 
 46     printf( "Size is %dx%dx%d
", 
 47             poDataset->GetRasterXSize(), poDataset->GetRasterYSize(),
 48             poDataset->GetRasterCount() );
 49 
 50     if( poDataset->GetProjectionRef()  != NULL )
 51         printf( "Projection is `%s'
", poDataset->GetProjectionRef() );
 52 
 53     if( poDataset->GetGeoTransform( adfGeoTransform ) == CE_None )
 54     {
 55         printf( "Origin = (%.6f,%.6f)
",
 56                 adfGeoTransform[0], adfGeoTransform[3] );
 57 
 58         printf( "Pixel Size = (%.6f,%.6f)
",
 59                 adfGeoTransform[1], adfGeoTransform[5] );
 60     }
 61 
 62     GDALRasterBand  *poBand;
 63     int             nBlockXSize, nBlockYSize;
 64     int             bGotMin, bGotMax;
 65     double          adfMinMax[2];
 66 
 67     poBand = poDataset->GetRasterBand( 1 );
 68     poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
 69     printf( "Block=%dx%d Type=%s, ColorInterp=%s
",
 70         nBlockXSize, nBlockYSize,
 71         GDALGetDataTypeName(poBand->GetRasterDataType()),
 72         GDALGetColorInterpretationName(
 73         poBand->GetColorInterpretation()) );
 74 
 75     adfMinMax[0] = poBand->GetMinimum( &bGotMin );
 76     adfMinMax[1] = poBand->GetMaximum( &bGotMax );
 77     if( ! (bGotMin && bGotMax) )
 78         GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);
 79 
 80     printf( "Min=%.3fd, Max=%.3f
", adfMinMax[0], adfMinMax[1] );
 81 
 82     if( poBand->GetOverviewCount() > 0 )
 83         printf( "Band has %d overviews.
", poBand->GetOverviewCount() );
 84 
 85     if( poBand->GetColorTable() != NULL )
 86         printf( "Band has a color table with %d entries.
", 
 87         poBand->GetColorTable()->GetColorEntryCount() );
 88     /*****
 89     Reading Raster Data
 90     *****/
 91     float *pafScanline;
 92     int nXSize = poBand->GetXSize();
 93     int nYSize=poBand->GetYSize();
 94 
 95     //读取图像的nXSize*nYSize数据
 96     pafScanline = (float*) CPLMalloc(sizeof(float)*nXSize*nYSize);//创建指针
 97     poBand->RasterIO(
 98         GF_Read,//第一个参数表示要读入数据还是写入数据
 99         0, 0,//nXOff, nYOff表示读取或者写入图像数据的起始坐标图像的左上角坐标为（0,0）
100         nXSize, nYSize,/*nXSize, nYSize表示读取或者写入图像数据的窗口大小,nXSize表示宽度,
101                        nYSize表示高度,均使用像素为单位,该宽度和高度是从第二个和第三个参数处开始计算。
102                        这两个参数和第二第三个参数一起表示就是,读取和写入图像的窗口位置和大小。*/
103         pafScanline, //指向存储数据的一个指针
104         nXSize, nYSize,//指定缓冲区的大小
105         GDT_Float32, 
106         0, 0 );
107     // poBand->RasterIO( GF_Read, 0, 0, nXSize, nYSize, pafScanline, nXSize, nYSize, GDT_Float32, 0, 0 );
108     cout<<"列数：nXSize="<<nXSize<<endl;    
109     cout<<"行数：nYSize="<<nYSize<<endl;
110     //system("pause");
111 
112 
113 
114     //创建nXSize*nYSize的float数组
115     float **data=new float *[nXSize];//每列有nXSize列,data数组的大小与dataLabel数组的大小相同,类型不同
116     int **dataLabel=new int *[nXSize];//创建标签数组,有nYSize行,nXSize列
117     int t=0,i=0,j=0;
118     for(int i=0;i<nYSize;i++)//nYSize行
119     {
120         dataLabel[i]=new int [nXSize];
121         data[i]     = new float [nXSize];
122     }
123 
124     cout<<"输出元数据数组data："<<endl;
125     for(int i=0;i<nYSize;i++)
126     {
127         for(int j=0;j<nXSize;j++)
128         {
129             data[i][j]=pafScanline[t++];
130             cout<<setw(4)<<data[i][j];
131         }
132         cout<<endl;
133     }
134     cout<<endl;
135     //system("pause");
136     cout<<"输出标签数组dataLabel数组："<<endl;
137     for(int i=0;i<nYSize;i++)
138     {
139         cout<<endl;
140         for(int j=0;j<nXSize;j++)
141         {
142             dataLabel[i][j]=nYSize*nXSize;
143             cout<<setw(4)<<dataLabel[i][j];
144         }
145     }
146     //system("pause");
147     cout<<endl;
148     cout<<"联通区域标记算法开始："<<endl;
149         dataLabel[0][0]=1;//把左上角的第一个标签赋值为1
150     int indexcolor=1;
151     for(i=0;i<nYSize;i++)
152     {
153         for(j=0;j<nXSize;j++)
154         {
155             if(i==0)
156             {
157                 if(j==0)
158                     continue;
159                 if(data[i][j]==data[i][j-1])
160                     dataLabel[i][j]=dataLabel[i][j-1];
161                 else
162                 {
163                     dataLabel[i][j]=++indexcolor;
164                 }
165                 continue;
166             }
167             if(j==0)
168             {
169                 if(data[i][j]==data[i-1][j])
170                 {
171                     if(dataLabel[i-1][j]<dataLabel[i][j])
172                         dataLabel[i][j]=dataLabel[i-1][j];
173                 }
174                 if(data[i][j]==data[i-1][j+1])
175                 {
176                     if(dataLabel[i-1][j+1]<dataLabel[i][j])
177                         dataLabel[i][j]=dataLabel[i-1][j+1];
178                 }
179                 if(dataLabel[i][j]==nYSize*nXSize)
180                 {
181                     dataLabel[i][j]=++indexcolor;
182                 }
183                 continue;
184             }
185             //if(edgeCheckL(i,j) && data[i][j]==data[i][j-1])//左边的,
186             if(j>0 && data[i][j]==data[i][j-1])
187             {
188                 if(dataLabel[i][j-1]<dataLabel[i][j])
189                     dataLabel[i][j]=dataLabel[i][j-1];
190 
191 
192             }
193             //if(edgeCheckLU(i,j) && data[i][j]==data[i-1][j-1])//左上角的
194             if((i>0 && j>0)&& data[i][j]==data[i-1][j-1])
195             {
196                 if(dataLabel[i-1][j-1]<dataLabel[i][j])
197                     dataLabel[i][j]=dataLabel[i-1][j-1];
198             }
199             //if(edgeCheckU(i,j) && data[i][j]==data[i-1][j])//上面的
200             if(i>0 && data[i][j]==data[i-1][j])
201             {
202                 if(dataLabel[i-1][j]<dataLabel[i][j])
203                     dataLabel[i][j]=dataLabel[i-1][j];
204             }
205             //if(edgeCheckUR(i,j) && data[i][j]==data[i-1][j+1])//右上角的
206             if((i>0 && j<nYSize-1) && data[i][j]==data[i-1][j+1])
207             {
208                 if(dataLabel[i-1][j+1]<dataLabel[i][j])
209                     dataLabel[i][j]=dataLabel[i-1][j+1];
210             }
211             if(dataLabel[i][j]==nYSize*nXSize)
212                 dataLabel[i][j]=++indexcolor;
213         }
214     }
215     //system("pause");
216     cout<<endl;
217     cout<<"输出第一次标记数组"<<endl<<endl;
218     for(i=0;i<nYSize;i++)
219     {
220         for(j=0;j<nXSize;j++)
221         {
222             cout<<setw(3)<<dataLabel[i][j];
223         }
224         cout<<endl;
225     }
226     //system("pause");
227     cout<<"合并首次生成的标签数组。。。。"<<endl<<endl;
228 
229     //合并：像素值相同，但是标签不同，就把标签值大的变为小的。
230     for(i=0;i<nYSize;i++)//行
231     {
232         for(j=0;j<nXSize;j++)//列
233         {
234             if(i==0)//第0行,只判左边的！
235             {
236                 if(j==0)//第一个元素
237                 {
238                     j=1;//跳过第一个元素，直接从第二个元素：data[0][1]判断
239                     //判断并执行合并：
240                     if(data[i][j-1]==data[i][j] && dataLabel[i][j-1]!=dataLabel[i][j])//像素值相同  && 标签不同：合并！
241                     {
242                         //把所有的大标签替换为当前两个标签中的一个较小的值
243                         //执行一次遍历
244                         //cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
245                         for(int m=0;m<nYSize;m++)
246                         {
247                             for(int n=0;n<nXSize;n++)
248                             {//把大的标签替换为小的标签
249                                 if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i][j-1]))
250                                 {
251                                     dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i][j-1]);
252                                 }
253                             }
254                         }
255                     }
256                 }
257             }
258             else//非0行
259             {
260                 if(j==0)//第0列，但不是第一行的：只判断上、右上两个方向
261                 {
262                     if(data[i-1][j]==data[i][j] && dataLabel[i-1][j]!=dataLabel[i][j])//像素值相同  && 标签不同：合并！
263                     {
264                         //把所有的大标签替换为当前两个标签中的一个较小的值
265                         //执行一次遍历
266                         //cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
267                         for(int m=0;m<nYSize;m++)
268                         {
269                             for(int n=0;n<nXSize;n++)
270                             {//把大的标签替换为小的标签
271                                 if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i-1][j]))
272                                 {
273                                     dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i-1][j]);
274                                 }
275                             }
276                         }
277                     }
278                     if(data[i-1][j+1]==data[i][j] && dataLabel[i-1][j=1]!=dataLabel[i][j])//像素值相同  && 标签不同：合并！
279                     {
280                         //把所有的大标签替换为当前两个标签中的一个较小的值
281                         //执行一次遍历
282                         //cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
283                         for(int m=0;m<nYSize;m++)
284                         {
285                             for(int n=0;n<nXSize;n++)
286                             {//把大的标签替换为小的标签
287                                 if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i-1][j+1]))
288                                 {
289                                     dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i-1][j+1]);
290                                 }
291                             }
292                         }
293                     }
294                 }
295                 else if(j==nYSize-1)//非0行且最后一列的：判断左、左上、上三个方向
296                 {
297                     if(data[i][j-1]==data[i][j] && dataLabel[i][j-1]!=dataLabel[i][j])//像素值相同  && 标签不同：合并！
298                     {
299                         //把所有的大标签替换为当前两个标签中的一个较小的值
300                         //执行一次遍历
301                         //cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
302                         for(int m=0;m<nYSize;m++)
303                         {
304                             for(int n=0;n<nXSize;n++)
305                             {//把大的标签替换为小的标签
306                                 if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i][j-1]))
307                                 {
308                                     dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i][j-1]);
309                                 }
310                             }
311                         }
312                     }
313                     if(data[i-1][j-1]==data[i][j] && dataLabel[i-1][j-1]!=dataLabel[i][j])//像素值相同  && 标签不同：合并！
314                     {
315                         //把所有的大标签替换为当前两个标签中的一个较小的值
316                         //执行一次遍历
317                         //cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
318                         for(int m=0;m<nYSize;m++)
319                         {
320                             for(int n=0;n<nXSize;n++)
321                             {//把大的标签替换为小的标签
322                                 if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i-1][-1]))
323                                 {
324                                     dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i-1][j-1]);
325                                 }
326                             }
327                         }
328                     }
329                     if(data[i-1][j]==data[i][j] && dataLabel[i-1][j]!=dataLabel[i][j])//像素值相同  && 标签不同：合并！
330                     {
331                         //把所有的大标签替换为当前两个标签中的一个较小的值
332                         //执行一次遍历
333                         //cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
334                         for(int m=0;m<nYSize;m++)
335                         {
336                             for(int n=0;n<nXSize;n++)
337                             {//把大的标签替换为小的标签
338                                 if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i-1][j]))
339                                 {
340                                     dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i-1][j]);
341                                 }
342                             }
343                         }
344                     }
345 
346                 }
347                 else//非0行且（既不是第一列，也不是最后一列的）：判断左、左上、上、右上四个方向
348                 {
349                     if(data[i][j-1]==data[i][j] && dataLabel[i][j-1]!=dataLabel[i][j])//像素值相同  && 标签不同：合并！
350                     {
351                         //把所有的大标签替换为当前两个标签中的一个较小的值
352                         //执行一次遍历
353                         //cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
354                         for(int m=0;m<nYSize;m++)
355                         {
356                             for(int n=0;n<nXSize;n++)
357                             {//把大的标签替换为小的标签
358                                 if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i][j-1]))
359                                 {
360                                     dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i][j-1]);
361                                 }
362                             }
363                         }
364                     }
365                     if(data[i-1][j-1]==data[i][j] && dataLabel[i-1][j-1]!=dataLabel[i][j])//像素值相同  && 标签不同：合并！
366                     {
367                         //把所有的大标签替换为当前两个标签中的一个较小的值
368                         //执行一次遍历
369                         //cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
370                         for(int m=0;m<nYSize;m++)
371                         {
372                             for(int n=0;n<nXSize;n++)
373                             {//把大的标签替换为小的标签
374                                 if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i-1][-1]))
375                                 {
376                                     dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i-1][j-1]);
377                                 }
378                             }
379                         }
380                     }
381                     if(data[i-1][j]==data[i][j] && dataLabel[i-1][j]!=dataLabel[i][j])//像素值相同  && 标签不同：合并！
382                     {
383                         //把所有的大标签替换为当前两个标签中的一个较小的值
384                         //执行一次遍历
385                         //cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
386                         for(int m=0;m<nYSize;m++)
387                         {
388                             for(int n=0;n<nXSize;n++)
389                             {//把大的标签替换为小的标签
390                                 if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i-1][j]))
391                                 {
392                                     dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i-1][j]);
393                                 }
394                             }
395                         }
396                     }
397                     if(data[i-1][j+1]==data[i][j] && dataLabel[i-1][j+1]!=dataLabel[i][j])//像素值相同  && 标签不同：合并！
398                     {
399                         //把所有的大标签替换为当前两个标签中的一个较小的值
400                         //执行一次遍历
401                         //cout<<"第"<<i<<"行"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
402                         for(int m=0;m<nYSize;m++)
403                         {
404                             for(int n=0;n<nXSize;n++)
405                             {//把大的标签替换为小的标签
406                                 if(dataLabel[m][n]==maxLabel(dataLabel[i][j],dataLabel[i-1][j+1]))
407                                 {
408                                     dataLabel[m][n]=minLabel(dataLabel[i][j],dataLabel[i-1][j+1]);
409                                 }
410                             }
411                         }
412                     }
413                 }
414             }
415 
416         }
417     }
418     t=1;
419     //system("pause");
420     cout<<"输出合并后的标签数组dataLabel:"<<endl;
421     for(i=0;i<nYSize;i++)
422     {
423         for(j=0;j<nXSize;j++)
424         {
425             if(dataLabel[i][j]>t)
426             {
427                 t++;
428             }
429             cout<<setw(3)<<dataLabel[i][j];
430 
431         }
432         cout<<endl;
433     }
434     cout<<endl<<"t="<<t<<endl;
435     cout<<"联通区域标记算法结束end！！！"<<endl;
436     //system("pause");
437     int max=0;
438     for(int i=0;i<nYSize;i++)
439     {
440         for(int j=0;j<nXSize;j++)
441         {
442             if(dataLabel[i][j]>max)
443             {
444                 max=dataLabel[i][j];
445             }
446         }
447     }
448     cout<<"maxLabel="<<max<<endl;
449     int *Pid=new int [max+1];
450     int *Parea=new int [max+1];
451     for(i=0;i<max+1;i++)
452     {
453         Pid[i]=i;
454         Parea[i]=0;
455     }
456     /////////////////////////////
457     for(int i=0;i<nYSize;i++)
458     {
459         for(int j=0;j<nXSize;j++)
460         {
461                 Parea[dataLabel[i][j]]++;
462         }
463     }
464     t=max+1;
465     for(i=1;i<t;i++)
466     {
467         
468         cout<<"dataLabel为"<<i<<"的面积为:"<<Parea[i]<<endl;
469     }
470     
471     double Xi=0,Si=0;
472     int NPatch=0;
473     for(i=0;i<t;i++)
474     {
475         if(Parea[i]!=0)
476         {
477             NPatch++;
478             Xi=Xi+Parea[i];//求出板块总面积
479         }
480     /*cout<<"NPatch="<<NPatch<<endl;
481     cout<<"Xi="<<Xi<<endl;*/
482     }
483     cout<<"NPatch="<<NPatch<<endl;
484     //cout<<"Xi="<<Xi<<endl;
485     Xi=Xi/NPatch;//面积平均数
486     cout<<"Initial Value:Xi="<<Xi<<endl;
487         cout<<"Initial Value:Si="<<Si<<endl;
488         //计算出所有
489     for(i=1;i<t;i++)
490     {
491         if(Parea[i]!=0)
492         {
493             Si=Si+((Parea[i]-Xi)*(Parea[i]-Xi))/NPatch;
494         }
495     }
496     cout<<"方差="<<Si<<endl;
497     Si=sqrt(Si);
498     cout<<"Standard Deviation(标准差)="<<Si<<endl;
499     cout<<endl;
500 
501     /****输出内容写出到文件当中****/
502     ofstream ocout;
503     //ocout.open("result.csv");
504     //ocout.open("result.txt");
505     ocout.open("result.csv");
506     /****将计算出的结果输出到屏幕****/
507     ocout<<"PatchID"<<","<<"Area"<<","<<"CSD"<<endl;
508     for(i=1;i<t;i++)
509     {
510         if(Parea[i]!=0)
511         {
512             
513             //ocout<<"Patch: Label="<<i<<setw(4)<<"  Area="<<Parea[i]<<setw(10)<<"CSD="<<(Parea[i]-Xi)/Si<<endl;
514             ocout<<i<<",";
515             ocout<<Parea[i]<<",";
516             ocout<<(Parea[i]-Xi)/Si<<","<<endl;
517         }
518     }
519     /*******后面的这一部分没有太大的实质性用处，可有可无**********/
520      const char *pszFormat = "GTiff";
521     GDALDriver *poDriver;
522     char **papszMetadata;
523 
524     poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
525 
526     if( poDriver == NULL )
527         exit( 1 );
528 
529     papszMetadata = poDriver->GetMetadata();
530     if( CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATE, FALSE ) )
531         printf( "Driver %s supports Create() method.
", pszFormat );
532     if( CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATECOPY, FALSE ) )
533         printf( "Driver %s supports CreateCopy() method.
", pszFormat );
534 
535     /*释放读取文件用的指针*/
536     CPLFree(pafScanline);
537     //关闭文件
538     GDALClose((GDALDatasetH)poDataset);
539     free(dataLabel);
540     free(data);
541     free(Parea);
542     free(Pid);
543     //
544     ocout.close();
545 }
546 int minLabel(int a,int b)
547 {
548     return (a>b)?b:a;
549 }
550 int maxLabel(int a,int b)
551 {
552     return (a<b)?b:a;
553 }
554 double csd(int *Parea,int length)
555 {
556     return 0;
557 }
View Code
【推广】 免费学中医，健康全家人
原文地址：https://www.cnblogs.com/AmatVictorialCuram/p/3584200.html