由于在一维数据中自然断点法和K-means算法效果一致,因此笔者使用C++语言基于Alglib(用于实现K-mans)和GDAL(处理Tiff格式数据)两个外部库实现Tiff图像的自然断点分类,并按升序分类以GByte格式输出结果。
处理后效果展示:
上代码:
#include "gdal_priv.h"
#include "ogrsf_frmts.h"
#include "stdafx.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "dataanalysis.h"
using namespace alglib;
using namespace std;
//基于Float/Double类型的多/单波段自然间断法分类
//initialTiff:原始数据存放路径 newTiff:自然断点法分类后数据存放路径 nInterval:分类数
void NaturalbreakpointTiff(std::string initialTiff, std::string newTiff, int nInterval)
{
GDALAllRegister();
//设置支持中文路径
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
const char * pszFile = initialTiff.c_str();
GDALDataset *poDataset = (GDALDataset*)GDALOpen(pszFile, GA_ReadOnly);//使用只读方式打开图像
if (!poDataset)
{
printf("File: %s不能打开!\n", pszFile);
}
//仿射参数
double padfTransform0[6];
if (poDataset->GetGeoTransform(padfTransform0) == CE_Failure)
{
printf("获取仿射变换参数失败");
}
int iImgSizeX0 = poDataset->GetRasterXSize();
int iImgSizeY0 = poDataset->GetRasterYSize();
int nCount = poDataset->GetRasterCount(); //影像的波段数
GDALDataType gdal_data_type = poDataset->GetRasterBand(1)->GetRasterDataType();//获取栅格类型
double dnodata = poDataset->GetRasterBand(1)->GetNoDataValue();//获取空值对应大小
int *panBandMap = new int[nCount];
for (int i = 0; i < nCount; i++)
{
panBandMap[i] = i + 1;
}
if(gdal_data_type==GDT_Float32)
{
float *pF_afScanblock1 = new float[iImgSizeX0*iImgSizeY0 * nCount];
poDataset->RasterIO(GF_Read, 0, 0, iImgSizeX0, iImgSizeY0, pF_afScanblock1, iImgSizeX0, iImgSizeY0, gdal_data_type, nCount, panBandMap, 0, 0, 0);
unsigned char *pUC_afScanblock1 = new unsigned char[iImgSizeX0*iImgSizeY0 * nCount];
for (int nBand = 0; nBand < nCount; nBand++)
{
int nEffectNum = 0;
for (int j = 0; j < iImgSizeX0*iImgSizeY0; j++)
{
if (abs((double)pF_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + j] - dnodata) > 0.000001)
{
nEffectNum++;
}
}
clusterizerstate s;
kmeansreport rep;
real_2d_array xy;
xy.setlength(nEffectNum, 1);
int nIndex = 0;
for (int j = 0; j < iImgSizeX0*iImgSizeY0; j++)
{
if (abs((double)pF_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + j] - dnodata) > 0.000001)
{
xy[nIndex][0] = (double)pF_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + j];
nIndex++;
}
}
clusterizercreate(s);
clusterizersetpoints(s, xy, 2);//2 Euclidean distance
clusterizersetkmeanslimits(s, 5, 0);
clusterizerrunkmeans(s, nInterval, rep);
double *pdValue = new double[nInterval];
double *pdValue_sort = new double[nInterval];
nIndex = 0;
for (int i = 0; i < iImgSizeX0*iImgSizeY0; i++)
{
if (abs((double)pF_afScanblock1[i] - dnodata) > 0.000001)
{
int a = rep.cidx[nIndex];
pdValue[rep.cidx[nIndex]] = pF_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + i];
pdValue_sort[rep.cidx[nIndex]] = pF_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + i];
nIndex++;
}
}
sort(pdValue_sort, pdValue_sort + nInterval);
vector<int>vNum;
for (int i = 0; i < nInterval; i++)
{
double _dValue = pdValue[i];
nIndex = 0;
while (abs(pdValue_sort[nIndex] - _dValue) > 0.00001)
{
nIndex++;
}
vNum.push_back(nIndex + 1);
}
nIndex = 0;
for (int i = 0; i < iImgSizeX0*iImgSizeY0; i++)
{
if (abs((double)pF_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + i] - dnodata) > 0.000001)
{
pUC_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + i] = vNum[rep.cidx[nIndex]];
nIndex++;
}
else
{
pUC_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + i] = 100;
}
}
delete[]pdValue_sort;
delete[]pdValue;
}
//写入新的tif
const char* pszDstFilename = newTiff.c_str();
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
GDALDataset* poDstDS = poDriver->Create(pszDstFilename, iImgSizeX0, iImgSizeY0, nCount, gdal_data_type, NULL);
poDstDS->SetProjection(poDataset->GetProjectionRef());//给它设置投影
poDstDS->SetGeoTransform(padfTransform0);//给设置空间转换的六参数
for (int i = 0; i < nCount; i++)
{
poDstDS->GetRasterBand(i + 1)->SetNoDataValue(100);//将空值设置为“无数据值”
}
//保存全波段影像
poDstDS->RasterIO(GF_Write, 0, 0, iImgSizeX0, iImgSizeY0, pUC_afScanblock1, iImgSizeX0, iImgSizeY0, GDT_Byte, nCount, panBandMap, 0, 0, 0);
//释放内存
GDALClose(poDstDS);
GDALClose(poDataset);
delete[]pF_afScanblock1;
delete[]pUC_afScanblock1;
}
else if (gdal_data_type == GDT_Float64)
{
double *pD_afScanblock1 = new double[iImgSizeX0*iImgSizeY0 *nCount];
poDataset->RasterIO(GF_Read, 0, 0, iImgSizeX0, iImgSizeY0, pD_afScanblock1, iImgSizeX0, iImgSizeY0, gdal_data_type, nCount, panBandMap, 0, 0, 0);
unsigned char *pUC_afScanblock1 = new unsigned char[iImgSizeX0*iImgSizeY0 * nCount];
for (int nBand = 0; nBand < nCount; nBand++)
{
int nEffectNum = 0;
for (int j = 0; j < iImgSizeX0*iImgSizeY0; j++)
{
if (abs((double)pD_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 +j] - dnodata) > 0.000001)
{
nEffectNum++;
}
}
clusterizerstate s;
kmeansreport rep;
real_2d_array xy;
xy.setlength(nEffectNum, 1);
int nIndex = 0;
for (int j = 0; j < iImgSizeX0*iImgSizeY0; j++)
{
if (abs((double)pD_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + j] - dnodata) > 0.000001)
{
xy[nIndex][0] = (double)pD_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + j];
nIndex++;
}
}
clusterizercreate(s);
clusterizersetpoints(s, xy, 2);//2 Euclidean distance
clusterizersetkmeanslimits(s, 5, 0);
clusterizerrunkmeans(s, nInterval, rep);
double *pdValue = new double[nInterval];
double *pdValue_sort = new double[nInterval];
nIndex = 0;
for (int i = 0; i < iImgSizeX0*iImgSizeY0; i++)
{
if (abs((double)pD_afScanblock1[i] - dnodata) > 0.000001)
{
int a = rep.cidx[nIndex];
pdValue[rep.cidx[nIndex]] = pD_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + i];
pdValue_sort[rep.cidx[nIndex]] = pD_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + i];
nIndex++;
}
}
sort(pdValue_sort, pdValue_sort + nInterval);
vector<int>vNum;
for (int i = 0; i < nInterval; i++)
{
double _dValue = pdValue[i];
nIndex = 0;
while (abs(pdValue_sort[nIndex] - _dValue) > 0.00001)
{
nIndex++;
}
vNum.push_back(nIndex + 1);
}
nIndex = 0;
for (int i = 0; i < iImgSizeX0*iImgSizeY0; i++)
{
if (abs((double)pD_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + i] - dnodata) > 0.000001)
{
pUC_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + i] = vNum[rep.cidx[nIndex]];
nIndex++;
}
else
{
pUC_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + i] = 100;
}
}
delete[]pdValue_sort;
delete[]pdValue;
}
//写入新的tif
const char* pszDstFilename = newTiff.c_str();
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
GDALDataset* poDstDS = poDriver->Create(pszDstFilename, iImgSizeX0, iImgSizeY0, nCount, gdal_data_type, NULL);
poDstDS->SetProjection(poDataset->GetProjectionRef());//给它设置投影
poDstDS->SetGeoTransform(padfTransform0);//给设置空间转换的六参数
for (int i = 0; i < nCount; i++)
{
poDstDS->GetRasterBand(i + 1)->SetNoDataValue(100);//将空值设置为“无数据值”
}
//保存全波段影像
poDstDS->RasterIO(GF_Write, 0, 0, iImgSizeX0, iImgSizeY0, pUC_afScanblock1, iImgSizeX0, iImgSizeY0, GDT_Byte, nCount, panBandMap, 0, 0, 0);
//释放内存
GDALClose(poDstDS);
GDALClose(poDataset);
delete[]pD_afScanblock1;
delete[]pUC_afScanblock1;
}
delete[]panBandMap;
}