基於GDAL的柵格圖像空間插值預處理
阿新 • • 發佈：2018-01-04
存儲處理 format value ast pri double 錯誤 依據 emd
轉自 基於GDAL的柵格圖像空間插值預處理——C語言版
基於GDAL的柵格圖像預處理
前言
柵格數據和矢量數據構成空間數據的主要來源，怎樣以開源方式讀取並處理這些空間數據？目前有多種開源支持包，這裏只介紹GDAL包。GDAL包的優點是支持庫簡潔、支持柵格和矢量、與多種開發平臺結合。OpenGis方式讀取空間數據，有利於自己編寫程序進行圖像預處理和智能識別等等，比如：遙感影像的降噪、銳化；紅外圖像的林火識別；工廠監控視頻識別等等。本文中利用GDAL包讀取高程柵格DEM，並添加氣象自動站點的數據，進行空間插值研究。
一、程序主要程序功能實現過程
第一步：讀取柵格數據，包含坡向和DEM
第二步：讀入站點信息數據
第三步：按照行列號讀取柵格單元到內存，不考慮高程為0的單元
第三步：每一個坡向情況中，參考固定數目的氣象站點。首先確定搜索範圍，獲取定量數目的監測站點。
第四步:在同一個計算窗口內，通過給各個因子賦權重，依據海拔高程回歸關系與加權回歸分析得到
溫度遞減率。
第五步：計算單個站點的溫度預測值，然後計算所有站點的距離權重因子，根據因子大小確定綜合
影響後的溫度值
第四步：開辟新內存存儲處理後的柵格數據，然後新建一個tiff格式的文件，把內存
數據導出到該文件中
二、代碼示例
int _tmain(intargc, _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"); }
三、插值結果
基於GDAL的柵格圖像空間插值預處理