基于GDAL的栅格图像预处理
前言
栅格数据和矢量数据构成空间数据的主要来源，怎样以开源方式读取并处理这些空间数据？目前有多种开源支持包，这里只介绍GDAL包。GDAL包的优点是支持库简洁、支持栅格和矢量、与多种开发平台结合。OpenGis方式读取空间数据，有利于自己编写程序进行图像预处理和智能识别等等，比如：遥感影像的降噪、锐化；红外图像的林火识别；工厂监控视频识别等等。本文中利用GDAL包读取高程栅格DEM，并添加气象自动站点的数据，进行空间插值研究。
一、程序主要程序功能实现过程
第一步：读取栅格数据，包含坡向和DEM
第二步：读入站点信息数据
第三步：按照行列号读取栅格单元到内存，不考虑高程为0的单元
第三步：每一个坡向情况中，参考固定数目的气象站点。首先确定搜索范围，获取定量数目的监测站点。
第四步:在同一个计算窗口内，通过给各个因子赋权重，依据海拔高程回归关系与加权回归分析得到
温度递减率。
第五步：计算单个站点的温度预测值，然后计算所有站点的距离权重因子，根据因子大小确定综合
影响后的温度值
第四步：开辟新内存存储处理后的栅格数据，然后新建一个tiff格式的文件，把内存
数据导出到该文件中
二、代码示例
int _tmain(int argc, _TCHAR* argv[])
{
/*
DEM、坡向栅格数据的数据框大小
*/
int XsizeDEM;
int YsizeDEM;
int XsizeAspect;
int YsizeAspect;
//double geoTransform[6];
//double Xp,Yp;
/*
DEM、坡向栅格单元对象的VALUE值
*/
short int *pmemDEM;
float *pmemAspect;
float *pmemNew;
//GDAL注册
GDALAllRegister();
/*
栅格单元的底图来源文件名：DEM/ASPECT
*/
const char *pszFileDEM;
pszFileDEM="F:\\beijing_dem\\bj25_CopyRaster11.img";
const char *pszFileAspect;
pszFileAspect="F:\\beijing_dem\\aspect_bj.tif";
/*读取站点信息*/
int n,m;
float site[32][6];
FILE *fp;
if((fp=fopen("F:\\beijing_dem\\site_36\\information.txt","r"))== NULL)
{
printf("cannot open this file\n");
exit(0);
}
for(n=0;n<32;n++) {
for(m=0;m<6;m++) {
fscanf(fp,"%f",&site[n][m]);
}
}
for(n=0;n<32;n++) {
for(m=0;m<6;m++) {
printf("%f",site[n][m]);/*注意这里*/;
printf(" ");
}
printf("\n"); //每输出一行，输出一个换行符
}
fclose(fp);
/*
DEM/ASPECT的数据集读取器
*/
GDALDataset *poDatasetDEM;
GDALRasterBand *poBandDEM;
GDALDataset *poDatasetAspect;
GDALRasterBand *poBandAspect;
/*
判断DEM/ASPECT文件是否存在，不存在错误提示
*/
poDatasetDEM=(GDALDataset*)GDALOpen(pszFileDEM,GA_ReadOnly);
if(poDatasetDEM==NULL)
{
printf("File： %s不能打开",pszFileDEM);
return 0;
}
poDatasetAspect=(GDALDataset*)GDALOpen(pszFileAspect,GA_ReadOnly);
if(poDatasetAspect==NULL)
{
printf("File： %s不能打开",pszFileAspect);
return 0;
}
/*
判断DEM/ASPECT文件如果存在，就把文件输入到波段第一层
*/
poBandDEM=poDatasetDEM->GetRasterBand(1);
poBandAspect=poDatasetAspect->GetRasterBand(1);
/*
被处理栅格单元对象的数据框大小的提取
*/
XsizeDEM=poBandDEM->GetXSize();
YsizeDEM=poBandDEM->GetYSize();
XsizeAspect=poBandAspect->GetXSize();
YsizeAspect=poBandAspect->GetYSize();
/*
被处理栅格单元对象的内存开辟
*/
pmemDEM=(short int *)CPLMalloc(sizeof(short int)*XsizeDEM*YsizeDEM);
poBandDEM->RasterIO(GF_Read,0,0,XsizeDEM,YsizeDEM,pmemDEM,XsizeDEM,YsizeDEM,GDT_Int16,0,0);
pmemAspect=(float*)CPLMalloc(sizeof(float)*XsizeAspect*YsizeAspect);
poBandAspect->RasterIO(GF_Read,0,0,XsizeAspect,YsizeAspect,pmemAspect,XsizeAspect,YsizeAspect,GDT_Float32,0,0);
/*
被处理栅格单元对象的类型提示
*/
printf("Type is: %s\n",GDALGetDataTypeName(poBandDEM->GetRasterDataType()));
//开辟新栅格内存空间
pmemNew=(float *)CPLMalloc(sizeof(float)*XsizeDEM*YsizeDEM);
for(int i=0;i<YsizeDEM;i++)
{
for(int j=0;j<XsizeDEM;j++)
{
int flag=0;
float H_value;
float A_value;
H_value=pmemDEM[i*XsizeDEM+j];
A_value=pmemAspect[i*XsizeDEM+j];
//单个栅格插值处理
//高程没有值的特殊情况
if(H_value==0)
{
pmemNew[i*XsizeDEM+j]=0;
}
else
{
//平地无坡向的情况
if(A_value==-1)
{
//处理站点和插值单元重合的情况，插值单元的值等于站点监测值
for(n=0;n<32;n++)
{
if(((site[n][4]-i)*(site[n][4]-i)+(site[n][5]-j)*(site[n][5]-j))==0)
{
float x1=site[n][0],x2=site[n][3];
pmemNew[i*XsizeDEM+j]=site[n][3];
printf("%f 平地站点数据： ",x1);
printf("%f\n",x2);
flag=1;
}
}
if(flag==0)
{
pmemNew[i*XsizeDEM+j]=plain_calculate(site,H_value,i,j);
//pmemNew[i*XsizeDEM+j]=3;
}
}
else
{
//处理站点和插值单元重合的情况，插值单元的值等于站点监测值
for(n=0;n<32;n++)
{
if(((site[n][4]-i)*(site[n][4]-i)+(site[n][5]-j)*(site[n][5]-j))==0)
{
float x1=site[n][0],x2=site[n][3];
pmemNew[i*XsizeDEM+j]=site[n][3];
printf("%f 非平地站点数据： ",x1);
printf("%f\n",x2);
flag=1;
}
}
if(flag==0)
{
//每一个栅格单元进行插值计算
pmemNew[i*XsizeDEM+j]=calculate(site,H_value,A_value,i,j);
//pmemNew[i*XsizeDEM+j]=1;
}
}
}
//poDatasetDEM->GetGeoTransform( geoTransform);
//Xp = geoTransform[0] +j*geoTransform[1]+i*geoTransform[2];
//Yp = geoTransform[3] + j*geoTransform[4] + i*geoTransform[5];
}
}
/**
创建新的空TIFF栅格文件
*/
GDALAllRegister();
GDALDriver *poDriver;
poDriver=GetGDALDriverManager()->GetDriverByName("GTiff");//AAIGrid(Arc/Info ASCII Grid) HFA (img no lim)
if(poDriver==NULL)
exit(1);
GDALDataset *poDstDS;
poDstDS=poDriver->Create("F:\\beijing_dem\\beijing_mix513.tiff",XsizeDEM,YsizeDEM,1,GDT_Float32,NULL);
double trans[6]={219323.300357,100.00,0.00,4680250.0000,0.0,-100.000};
//如果图像不含地理坐标信息，默认返回值是：(0,1,0,0,0,1)
//In a north up image,
//左上角点坐标(padfGeoTransform[0],padfGeoTransform[3])；
//padfGeoTransform[1]是像元宽度(影像在宽度上的分辨率)；
//padfGeoTransform[5]是像元高度(影像在高度上的分辨率)；
//如果影像是指北的,padfGeoTransform[2]和padfGeoTransform[4]这两个参数的值为0。
poDstDS->SetGeoTransform(trans);
GDALClose((GDALDatasetH)poDstDS);
/*
//处理后的数据层保存
*/
GDALDataset *poDatasetNew;
poDatasetNew=(GDALDataset*)GDALOpen("F:\\beijing_dem\\beijing_mix513.tiff",GA_Update);
GDALRasterBand *poBandNew;
poBandNew=poDatasetNew->GetRasterBand(1);
poBandNew->RasterIO(GF_Write,0,0,XsizeDEM,YsizeDEM,pmemNew,XsizeDEM,YsizeDEM,GDT_Float32,0,0);
GDALClose((GDALDatasetH)poDatasetNew);
//释放数据
CPLFree(pmemDEM);
CPLFree(pmemNew);
printf("处理结束\n");
}三、插值结果
