基於GDAL的柵格圖像空間插值預處理

轉自 基於GDAL的柵格圖像空間插值預處理——C語言版

 

基於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");  
      
      
      
      
    }  

 

 三、插值結果