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

    }

三、插值结果

原文地址：https://www.cnblogs.com/arxive/p/8192245.html

时间： 2024-08-28 18:39:38

基于GDAL的栅格图像空间插值预处理的相关文章

关于基于GDAL库QT软件平台下C++语言开发使用说明

背景前提 地理空间数据抽象库(GDAL)是一个用于读取和编写栅格和矢量地理空间数据格式的计算机软件库,由开源地理空间基金会在许可的X / MIT风格免费软件许可下发布. 作为一个库,它为调用应用程序提供了一个抽象数据模型,用于所有支持的格式. 它还可以构建有各种有用的命令行接口实用程序,用于数据转换和处理. PROJ.4库支持投影和转换.(摘自维基百科) 相关的OGR库(OGR Simple Features Library [2])是GDAL源代码树的一部分,它为简单的特征矢量图形数据提供了类

Tetrahedron based light probe interpolation(基于四面体的Light Probe插值)

在当前的游戏引擎中,使用Light Probe来计算全局环境光对于动态物体的影响是一种很主流的方法.在预处理阶段生成完场景的Light Probe之后,传统的方法采用查找最近的8个相邻的Probe然后使用三线性的方式(Trilinear Interpolation)进行插值,但是这样的插值代价稍大,不过一个可行的优化就是尽可能地减少插值中使用的Probe的数量,比如由8个减少到4个不等.但是这时就不能再用三线性插值,而需要用其它的插值方法比如Inverse Distance等,不过这样就会带来另

基于Unity3D三维模型的动作插值（空间关键帧动画实现）

1.引言 最近在Unity3D中实现一个基于自定义Mesh网格的骨骼动画.存储关键帧信息,然后通过插值形成中间动画.网格GameObject之间存在父子关系.插值动画对模型骨骼的Position.Sclae.Rotation三个部分分别混合插值. 并且注意,一般选取的是子骨骼相对父骨骼的Transform信息,即上述关键帧信息具体应该为Transform.localPosition.localScale.localRotation. 并且这个系统的关键帧是定义在三维空间中的,用户在一系列关键帧点

基于GDAL的点数据保存

保存数据的第一步是要把数据解析出来,然后根据GDAL的规则进行数据point类型的shapefile数据生成.大概步骤为: 一.定义保存点要素数据的类 这里定义了两个基类: //基类,保存要素类型,点.线.面 class Element { private: char Type; int Code; public: Element(void); ~Element(void); void setType(char Type); char getType(); void setCode(int Co

基于GDAL的线数据保存

保存数据的第一步是要把数据解析出来,然后根据GDAL的规则进行数据point类型的shapefile数据生成.大概步骤为: 一.定义保存点要素数据的类 这里定义了两个基类: //基类,保存要素类型,点.线.面 class Element { private: char Type; int Code; public: Element(void); ~Element(void); void setType(char Type); char getType(); void setCode(int Co

基于GDAL的面数据保存

保存数据的第一步是要把数据解析出来,然后根据GDAL的规则进行数据point类型的shapefile数据生成.大概步骤为: 一.定义保存点要素数据的类 这里定义了两个基类: //基类,保存要素类型,点.线.面 class Element { private: char Type; int Code; public: Element(void); ~Element(void); void setType(char Type); char getType(); void setCode(int Co

基于GDAL库，读取.grd文件（以海洋地形数据为例）C++版

技术背景 海洋地形数据主要是通过美国全球地形起伏数据(GMT)获得,数据格式为grd(GSBG)二进制数据,打开软件通过是Surfer软件,surfer软件可进行数据的编辑处理,以及进一步的可视化表达等功能操作:由于Surfer软件不支持二次开发,没有提供相应的SDK供开发者进行使用,所以这一切只能通过相应类似的技术进行实现,首先,数据的读取,如何通过编程实现数据的读取操作呢?这里就要说一下GIS软件所使用的一个开源库-GDAL,GDAL库的具体解释资料,请查阅官方网站[https://www.

GDAL python教程（1）——用OGR读写矢量数据

本教程的讲义和源码都是取自Utah State University的openGIS课程 相关资料,包括讲义.源码.数据样例,请从此处下载http://www.gis.usu.edu/~chrisg/python/ 本人只是做点翻译,写写学习体会而已,版权属于原作者. 欢迎转载,不过别忘了上面这段话. ================================================== 为什么用open source? 优点 免费,适合个人和小公司 强大的开发工具,找bug更容易

基于Python数据分析与机器学习案例实战教程

课程--基于Python数据分析与机器学习案例实战教程 分享网盘下载--https://pan.baidu.com/s/1jHSaRAY 密码: xk37 课程背景基于数据分析与机器学习领域,使用python作为课程的实战语言,随着大数据与人工智能领域日益火爆,数据分析和机器学习建模成了当下最热门的技术,课程旨在帮助同学们快速掌握python数据分析包以及经典机器学习算法并通过对真实数据集分析进行实战演示. 课程风格通俗易懂,基于真实数据集案例实战. 主体课程分成三个大模块 (1)python数