GDAL栅格API使用
读取栅格数据
void main()
{
// 注册所有驱动
GDALAllRegister();
// 打开调试信息
CPLSetConfigOption("CPL_DEBUG", "ON");
// 准备数据源
std::string url = u8"C:/Users/admin/Desktop/temp/123.tif";
// 打开栅格数据
GDALDataset *raster = (GDALDataset*)GDALOpenEx(url.c_str(), GDAL_OF_ALL, nullptr, nullptr, nullptr);
if(nullptr == raster)
{
std::cout << CPLGetLastErrorMsg() << std::endl;
return;
}
// 打印栅格驱动信息
printf("Driver: %s/%s\n", raster->GetDriver()->GetDescription(),
raster->GetDriver()->GetMetadataItem(GDAL_DMD_LONGNAME));
// 打印栅格维度(宽、高、波段)
printf("Size is %dx%dx%d\n", raster->GetRasterXSize(),
raster->GetRasterYSize(),
raster->GetRasterCount());
// 打印坐标系信息
if (raster->GetProjectionRef() != NULL)
printf("Projection is `%s'\n", raster->GetProjectionRef());
// 参考变换规则
// Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2)
// Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5)
// In case of north up images,
// the GT(2) and GT(4) coefficients are zero,
// and the GT(1) is pixel width, and GT(5) is pixel height.
// The (GT(0),GT(3)) position is the top left corner of the top left pixel of the raster.
double adfGeoTransform[6];
if (raster->GetGeoTransform(adfGeoTransform) == CE_None)
{
// 打印原点信息
printf("Origin = (%.6f,%.6f)\n",
adfGeoTransform[0], adfGeoTransform[3]);
// 打印像素大小
printf("Pixel Size = (%.6f,%.6f)\n",
adfGeoTransform[1], adfGeoTransform[5]);
}
//遍历波段,下标从1开始
for(int i = 1; i <= raster->GetRasterCount(); i++)
{
// 得到波段,poBand属于dataset，获取后无需释放
GDALRasterBand *poBand = raster->GetRasterBand(i);
// 获取波段的块大小
// 获取波段的数据类型（Byte,UInt16...）
// 获取波段的颜色解释(灰度图、RGB图等)
int nBlockXSize, nBlockYSize;
poBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
printf("Block=%dx%d Type=%s, ColorInterp=%s\n",
nBlockXSize, nBlockYSize,
GDALGetDataTypeName(poBand->GetRasterDataType()),
GDALGetColorInterpretationName(
poBand->GetColorInterpretation()));
// 获取波段最大最小值
int bGotMin, bGotMax;
double adfMinMax[2];
adfMinMax[0] = poBand->GetMinimum(&bGotMin);
adfMinMax[1] = poBand->GetMaximum(&bGotMax);
if(!(bGotMin && bGotMax))
{
//如果获取失败，则计算
GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);
}
printf("Min=%f, Max=%f\n", adfMinMax[0], adfMinMax[1]);
// 金字塔？
if (poBand->GetOverviewCount() > 0)
printf("Band has %d overviews.\n", poBand->GetOverviewCount());
// 获取颜色表
if (poBand->GetColorTable() != NULL)
printf("Band has a color table with %d entries.\n",
poBand->GetColorTable()->GetColorEntryCount());
/*******************************************************
* 读取栅格数据
*******************************************************/
unsigned char *pafScanline;
// 获取波段宽度
int nXSize = poBand->GetXSize();
// 根据宽度开辟缓冲区
pafScanline = (unsigned char *)CPLMalloc(sizeof(unsigned char) * nXSize);
//
/**
*
* RasterIO参数
* CPLErr GDALRasterBand::RasterIO
* (
* GDALRWFlag eRWFlag, //读写标志位:GF_Read or GF_Write
* int nXOff, int nYOff, int nXSize, int nYSize, //选定需要读/写的窗口
* void * pData, // 读取/写入后的数据缓存,其类型必须是参数eBufType的类型,如GDT_Float32,GDT_Byte
* int nBufXSize, int nBufYSize, //缓冲的大小，如果比选定的窗口小，则是低分辨率,此时会利用OverViews做IO操作，以提高效率
* GDALDataType eBufType, //传入的缓冲区类型
* int nPixelSpace,
* int nLineSpace //nPixelSpace,nLineSpace 通常是0， 指明使用的默认值
* )
*/
// 读取第一行数据
poBand->RasterIO(GF_Read,
0, 0, nXSize, 1,
pafScanline,
nXSize, 1,
GDT_Byte, 0, 0);
for(int i = 0; i < nXSize; i++)
{
CPLDebug("TEST", "%d",pafScanline[i]);
}
//读取的数据使用完后释放
CPLFree(pafScanline);
}
//关闭数据源
GDALClose(raster);
}
写入栅格数据
有两种方法可以创建文件:CreateCopy()和Create()，所有支持创建文件的驱动都支持CreateCopy()，但只有一部分驱动支持Create()，通过检查 DCAP_CREATE 和 DCAP_CREATECOPY 元数据来判断该驱动是否支持。
CreateCopy()赋值栅格数据
void main()
{
// 注册所有驱动
GDALAllRegister();
// 打开调试信息
CPLSetConfigOption("CPL_DEBUG", "ON");
// 有两种方法可以创建文件:CreateCopy()和Create()
// 所有支持创建文件的驱动都支持CreateCopy()，但只有一部分驱动支持Create()
// 通过检查 DCAP_CREATE 和 DCAP_CREATECOPY 元数据来判断该驱动是否支持
const char *pszFormat = "GTiff";
GDALDriver *poDriver;
char **papszMetadata;
// 获取GTiff驱动
poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
if (poDriver == NULL)
exit(1);
// 获取驱动的元数据
papszMetadata = poDriver->GetMetadata();
// 从元数据中检查是否支持Create()
if (CSLFetchBoolean(papszMetadata, GDAL_DCAP_CREATE, FALSE))
printf("Driver %s supports Create() method.\n", pszFormat);
// 从元数据中检查是否支持CreateCopy()
if (CSLFetchBoolean(papszMetadata, GDAL_DCAP_CREATECOPY, FALSE))
printf("Driver %s supports CreateCopy() method.\n", pszFormat);
/********************************************
* Use CreateCopy()
*********************************************/
const char *pszSrcFilename = u8"C:/Users/admin/Desktop/temp/123.tif";
const char *pszDstFilename = u8"C:/Users/admin/Desktop/temp/13/copywithoutoptions.tif";
// 打开原始文件
GDALDataset *poSrcDS =
(GDALDataset *)GDALOpen(pszSrcFilename, GA_ReadOnly);
// 使用默认设置赋值到新文件
GDALDataset *poDstDS;
poDstDS = poDriver->CreateCopy(pszDstFilename, poSrcDS,
FALSE, //设为TRUE则表示复制必须严格相等，通常是FALSE以按需适配输出的格式
NULL, //papszOptions: 创建的设置选项
NULL, //pfnProgress:回调函数，展示复制进度
NULL); //pProgressData:传给回调函数的数据
/* Once we're done, close properly the dataset */
if (poDstDS != NULL)
GDALClose((GDALDatasetH)poDstDS);
GDALClose((GDALDatasetH)poSrcDS);
}
Create()创建栅格数据
void main()
{
// 注册所有驱动
GDALAllRegister();
// 打开调试信息
CPLSetConfigOption("CPL_DEBUG", "ON");
const char *pszFormat = "GTiff";
GDALDriver *poDriver;
char **papszMetadata;
// 获取GTiff驱动
poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
if (poDriver == NULL)
exit(1);
// 获取驱动的元数据
papszMetadata = poDriver->GetMetadata();
// 从元数据中检查是否支持Create()
if (CSLFetchBoolean(papszMetadata, GDAL_DCAP_CREATE, FALSE))
printf("Driver %s supports Create() method.\n", pszFormat);
// 从元数据中检查是否支持CreateCopy()
if (CSLFetchBoolean(papszMetadata, GDAL_DCAP_CREATECOPY, FALSE))
printf("Driver %s supports CreateCopy() method.\n", pszFormat);
const char *pszDstFilename = u8"C:/Users/admin/Desktop/temp/13/Create.tif";
GDALDataset *poDstDS;
char **papszOptions = NULL;
poDstDS = poDriver->Create(pszDstFilename,
10, //nXSize:宽度
10, //nYSize:高度
1, //nBands:波段数
GDT_Byte, //栅格类型
papszOptions);
// 参考变换规则
// Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2)
// Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5)
// In case of north up images,
// the GT(2) and GT(4) coefficients are zero,
// and the GT(1) is pixel width, and GT(5) is pixel height.
// The (GT(0),GT(3)) position is the top left corner of the top left pixel of the raster.
double adfGeoTransform[6] = {444720, 30, 0, 3751320, 0, -30};
OGRSpatialReference oSRS;
char *pszSRS_WKT = NULL;
GDALRasterBand *poBand;
GByte abyRaster[10 * 10];
// 设置转换信息
poDstDS->SetGeoTransform(adfGeoTransform);
// 设置空间参考
oSRS.SetUTM(11, TRUE);
oSRS.SetWellKnownGeogCS("NAD27");
oSRS.exportToWkt(&pszSRS_WKT);
poDstDS->SetProjection(pszSRS_WKT);
CPLFree(pszSRS_WKT);
// 获取波段
poBand = poDstDS->GetRasterBand(1);
// 往波段写入信息
for(int i = 0; i < 10*10; i++)
{
abyRaster[i] = rand()%(255-0)+0;
}
poBand->RasterIO(GF_Write,
0, 0, 5, 5, //写入的窗口大小
abyRaster, 10, 10, //待写入的缓冲区
GDT_Byte, //写入的数据类型
0, 0);
/* Once we're done, close properly the dataset */
GDALClose((GDALDatasetH)poDstDS);
}
网友评论