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

C++调⽤GDAL库读取并输出tif⽂件，并计算斑块⾯积输出景观
指数：ＣＳＤ
部分源码选⾃GDAL库的官⽅⽹址：，其余的代码为笔者⾃⼰编写。

1// readfile.cpp : 定义控制台应⽤程序的⼊⼝点。

2//
3/*
4part of the codes were cite from: /gdal_tutorial.html
5and remaining of code were created :by /AmatVictorialCuram/
6and 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>
15using namespace std;
16/**/
17int minLabel(int a,int b);
18
19int maxLabel(int a,int b);
20
21double csd(int *Parea,int length);
22
23int main()
24 {
25/*
26 part of the codes were cite from: /gdal_tutorial.html
27 and remaining of code were created :by /AmatVictorialCuram/
28 and please mark the site if you cite the program fo commercial of publication.
29*/
30const char *pszFilename="E:\\tif/fragstats/sample.tif";
31 GDALDataset *poDataset;
32
33 GDALAllRegister();
34
35 poDataset = (GDALDataset *) GDALOpen( pszFilename, GA_ReadOnly );
36if( poDataset == NULL )
37 {
38 cout<<"file read error!"<<endl;;
39 }
40double adfGeoTransform[6];
41
42 printf( "Driver: %s/%s\n",
43 poDataset->GetDriver()->GetDescription(),
44 poDataset->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME ) );
45
46 printf( "Size is %dx%dx%d\n",
47 poDataset->GetRasterXSize(), poDataset->GetRasterYSize(),
48 poDataset->GetRasterCount() );
49
50if( poDataset->GetProjectionRef() != NULL )
51 printf( "Projection is `%s'\n", poDataset->GetProjectionRef() );
52
53if( poDataset->GetGeoTransform( adfGeoTransform ) == CE_None )
54 {
55 printf( "Origin = (%.6f,%.6f)\n",
56 adfGeoTransform[0], adfGeoTransform[3] );
57
58 printf( "Pixel Size = (%.6f,%.6f)\n",
59 adfGeoTransform[1], adfGeoTransform[5] );
60 }
61
62 GDALRasterBand *poBand;
63int nBlockXSize, nBlockYSize;
64int bGotMin, bGotMax;
65double adfMinMax[2];
66
67 poBand = poDataset->GetRasterBand( 1 );
68 poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
69 printf( "Block=%dx%d Type=%s, ColorInterp=%s\n",
70 nBlockXSize, nBlockYSize,
71 GDALGetDataTypeName(poBand->GetRasterDataType()),
77if( ! (bGotMin && bGotMax) )
78 GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);
79
80 printf( "Min=%.3fd, Max=%.3f\n", adfMinMax[0], adfMinMax[1] );
81
82if( poBand->GetOverviewCount() > 0 )
83 printf( "Band has %d overviews.\n", poBand->GetOverviewCount() );
84
85if( poBand->GetColorTable() != NULL )
86 printf( "Band has a color table with %d entries.\n",
87 poBand->GetColorTable()->GetColorEntryCount() );
88/*****
89 Reading Raster Data
90 *****/
91float *pafScanline;
92int nXSize = poBand->GetXSize();
93int nYSize=poBand->GetYSize();
94
95//读取图像的nXSize*nYSize数据
96 pafScanline = (float*) CPLMalloc(sizeof(float)*nXSize*nYSize);//创建指针
97 poBand->RasterIO(
98 GF_Read,//第⼀个参数表⽰要读⼊数据还是写⼊数据
990, 0,//nXOff, nYOff表⽰读取或者写⼊图像数据的起始坐标图像的左上⾓坐标为（0,0）
100 nXSize, nYSize,/*nXSize, nYSize表⽰读取或者写⼊图像数据的窗⼝⼤⼩,nXSize表⽰宽度,
101 nYSize表⽰⾼度,均使⽤像素为单位,该宽度和⾼度是从第⼆个和第三个参数处开始计算。

102这两个参数和第⼆第三个参数⼀起表⽰就是,读取和写⼊图像的窗⼝位置和⼤⼩。

*/
103 pafScanline, //指向存储数据的⼀个指针
104 nXSize, nYSize,//指定缓冲区的⼤⼩
105 GDT_Float32,
1060, 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数组
115float **data=new float *[nXSize];//每列有nXSize列,data数组的⼤⼩与dataLabel数组的⼤⼩相同,类型不同116int **dataLabel=new int *[nXSize];//创建标签数组,有nYSize⾏,nXSize列
117int t=0,i=0,j=0;
118for(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;
125for(int i=0;i<nYSize;i++)
126 {
127for(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;
137for(int i=0;i<nYSize;i++)
138 {
139 cout<<endl;
140for(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
150int indexcolor=1;
151for(i=0;i<nYSize;i++)
152 {
153for(j=0;j<nXSize;j++)
154 {
155if(i==0)
156 {
161else
162 {
163 dataLabel[i][j]=++indexcolor;
164 }
165continue;
166 }
167if(j==0)
168 {
169if(data[i][j]==data[i-1][j])
170 {
171if(dataLabel[i-1][j]<dataLabel[i][j])
172 dataLabel[i][j]=dataLabel[i-1][j];
173 }
174if(data[i][j]==data[i-1][j+1])
175 {
176if(dataLabel[i-1][j+1]<dataLabel[i][j])
177 dataLabel[i][j]=dataLabel[i-1][j+1];
178 }
179if(dataLabel[i][j]==nYSize*nXSize)
180 {
181 dataLabel[i][j]=++indexcolor;
182 }
183continue;
184 }
185//if(edgeCheckL(i,j) && data[i][j]==data[i][j-1])//左边的,
186if(j>0 && data[i][j]==data[i][j-1])
187 {
188if(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])//左上⾓的
194if((i>0 && j>0)&& data[i][j]==data[i-1][j-1])
195 {
196if(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])//上⾯的
200if(i>0 && data[i][j]==data[i-1][j])
201 {
202if(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])//右上⾓的
206if((i>0 && j<nYSize-1) && data[i][j]==data[i-1][j+1])
207 {
208if(dataLabel[i-1][j+1]<dataLabel[i][j])
209 dataLabel[i][j]=dataLabel[i-1][j+1];
210 }
211if(dataLabel[i][j]==nYSize*nXSize)
212 dataLabel[i][j]=++indexcolor;
213 }
214 }
215//system("pause");
216 cout<<endl;
217 cout<<"输出第⼀次标记数组"<<endl<<endl;
218for(i=0;i<nYSize;i++)
219 {
220for(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//合并：像素值相同，但是标签不同，就把标签值⼤的变为⼩的。

230for(i=0;i<nYSize;i++)//⾏
231 {
232for(j=0;j<nXSize;j++)//列
233 {
234if(i==0)//第0⾏,只判左边的！
235 {
236if(j==0)//第⼀个元素
237 {
238 j=1;//跳过第⼀个元素，直接从第⼆个元素：data[0][1]判断
239//判断并执⾏合并：
240if(data[i][j-1]==data[i][j] && dataLabel[i][j-1]!=dataLabel[i][j])//像素值相同 && 标签不同：合并！
245for(int m=0;m<nYSize;m++)
246 {
247for(int n=0;n<nXSize;n++)
248 {//把⼤的标签替换为⼩的标签
249if(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 }
258else//⾮0⾏
259 {
260if(j==0)//第0列，但不是第⼀⾏的：只判断上、右上两个⽅向
261 {
262if(data[i-1][j]==data[i][j] && dataLabel[i-1][j]!=dataLabel[i][j])//像素值相同 && 标签不同：合并！263 {
264//把所有的⼤标签替换为当前两个标签中的⼀个较⼩的值
265//执⾏⼀次遍历
266//cout<<"第"<<i<<"⾏"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
267for(int m=0;m<nYSize;m++)
268 {
269for(int n=0;n<nXSize;n++)
270 {//把⼤的标签替换为⼩的标签
271if(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 }
278if(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;
283for(int m=0;m<nYSize;m++)
284 {
285for(int n=0;n<nXSize;n++)
286 {//把⼤的标签替换为⼩的标签
287if(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 }
295else if(j==nYSize-1)//⾮0⾏且最后⼀列的：判断左、左上、上三个⽅向
296 {
297if(data[i][j-1]==data[i][j] && dataLabel[i][j-1]!=dataLabel[i][j])//像素值相同 && 标签不同：合并！298 {
299//把所有的⼤标签替换为当前两个标签中的⼀个较⼩的值
300//执⾏⼀次遍历
301//cout<<"第"<<i<<"⾏"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
302for(int m=0;m<nYSize;m++)
303 {
304for(int n=0;n<nXSize;n++)
305 {//把⼤的标签替换为⼩的标签
306if(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 }
313if(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;
318for(int m=0;m<nYSize;m++)
319 {
320for(int n=0;n<nXSize;n++)
321 {//把⼤的标签替换为⼩的标签
322if(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]);
329if(data[i-1][j]==data[i][j] && dataLabel[i-1][j]!=dataLabel[i][j])//像素值相同 && 标签不同：合并！330 {
331//把所有的⼤标签替换为当前两个标签中的⼀个较⼩的值
332//执⾏⼀次遍历
333//cout<<"第"<<i<<"⾏"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
334for(int m=0;m<nYSize;m++)
335 {
336for(int n=0;n<nXSize;n++)
337 {//把⼤的标签替换为⼩的标签
338if(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 }
347else//⾮0⾏且（既不是第⼀列，也不是最后⼀列的）：判断左、左上、上、右上四个⽅向
348 {
349if(data[i][j-1]==data[i][j] && dataLabel[i][j-1]!=dataLabel[i][j])//像素值相同 && 标签不同：合并！350 {
351//把所有的⼤标签替换为当前两个标签中的⼀个较⼩的值
352//执⾏⼀次遍历
353//cout<<"第"<<i<<"⾏"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
354for(int m=0;m<nYSize;m++)
355 {
356for(int n=0;n<nXSize;n++)
357 {//把⼤的标签替换为⼩的标签
358if(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 }
365if(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;
370for(int m=0;m<nYSize;m++)
371 {
372for(int n=0;n<nXSize;n++)
373 {//把⼤的标签替换为⼩的标签
374if(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 }
381if(data[i-1][j]==data[i][j] && dataLabel[i-1][j]!=dataLabel[i][j])//像素值相同 && 标签不同：合并！382 {
383//把所有的⼤标签替换为当前两个标签中的⼀个较⼩的值
384//执⾏⼀次遍历
385//cout<<"第"<<i<<"⾏"<<"第"<<j<<"列："<<"第"<<++replace<<"次替换"<<endl;
386for(int m=0;m<nYSize;m++)
387 {
388for(int n=0;n<nXSize;n++)
389 {//把⼤的标签替换为⼩的标签
390if(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 }
397if(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;
402for(int m=0;m<nYSize;m++)
403 {
404for(int n=0;n<nXSize;n++)
405 {//把⼤的标签替换为⼩的标签
406if(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]);
414 }
415
416 }
417 }
418 t=1;
419//system("pause");
420 cout<<"输出合并后的标签数组dataLabel:"<<endl;
421for(i=0;i<nYSize;i++)
422 {
423for(j=0;j<nXSize;j++)
424 {
425if(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");
437int max=0;
438for(int i=0;i<nYSize;i++)
439 {
440for(int j=0;j<nXSize;j++)
441 {
442if(dataLabel[i][j]>max)
443 {
444 max=dataLabel[i][j];
445 }
446 }
447 }
448 cout<<"maxLabel="<<max<<endl;
449int *Pid=new int [max+1];
450int *Parea=new int [max+1];
451for(i=0;i<max+1;i++)
452 {
453 Pid[i]=i;
454 Parea[i]=0;
455 }
456/////////////////////////////
457for(int i=0;i<nYSize;i++)
458 {
459for(int j=0;j<nXSize;j++)
460 {
461 Parea[dataLabel[i][j]]++;
462 }
463 }
464 t=max+1;
465for(i=1;i<t;i++)
466 {
467
468 cout<<"dataLabel为"<<i<<"的⾯积为:"<<Parea[i]<<endl; 469 }
470
471double Xi=0,Si=0;
472int NPatch=0;
473for(i=0;i<t;i++)
474 {
475if(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//计算出所有
489for(i=1;i<t;i++)
490 {
491if(Parea[i]!=0)
492 {
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;
508for(i=1;i<t;i++)
509 {
510if(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/*******后⾯的这⼀部分没有太⼤的实质性⽤处，可有可⽆**********/
520const char *pszFormat = "GTiff";
521 GDALDriver *poDriver;
522char **papszMetadata;
523
524 poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
525
526if( poDriver == NULL )
527 exit( 1 );
528
529 papszMetadata = poDriver->GetMetadata();
530if( CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATE, FALSE ) )
531 printf( "Driver %s supports Create() method.\n", pszFormat );
532if( CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATECOPY, FALSE ) )
533 printf( "Driver %s supports CreateCopy() method.\n", 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 }
546int minLabel(int a,int b)
547 {
548return (a>b)?b:a;
549 }
550int maxLabel(int a,int b)
551 {
552return (a<b)?b:a;
553 }
554double csd(int *Parea,int length)
555 {
556return0;
557 }
View Code。