基于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"); }
三、插值结果
