开发者社区> 橘子红了呐> 正文
阿里云
为了无法计算的价值
打开APP
阿里云APP内打开

基于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");  
      
      
      
      
    }  
复制代码

 

 三、插值结果

 

 

 

没有整理与归纳的知识,一文不值!高度概括与梳理的知识,才是自己真正的知识与技能。 永远不要让自己的自由、好奇、充满创造力的想法被现实的框架所束缚,让创造力自由成长吧! 多花时间,关心他(她)人,正如别人所关心你的。理想的腾飞与实现,没有别人的支持与帮助,是万万不能的。





  本文转自wenglabs博客园博客,原文链接:http://www.cnblogs.com/arxive/p/8192245.html,如需转载请自行联系原作者


版权声明:本文内容由阿里云实名注册用户自发贡献,版权归原作者所有,阿里云开发者社区不拥有其著作权,亦不承担相应法律责任。具体规则请查看《阿里云开发者社区用户服务协议》和《阿里云开发者社区知识产权保护指引》。如果您发现本社区中有涉嫌抄袭的内容,填写侵权投诉表单进行举报,一经查实,本社区将立刻删除涉嫌侵权内容。

相关文章
gdb kernel debug的进程断点
gdb调试kernel的时候, 如果设置通用函数断点, 比如vfs_read, 就会遇到一堆撞到断点的地方, 比如tty输入一个字符, 就是vfs_read, 没办法调试具体的某一个进程 一种办法就是条件断点, 其实不是很好用, 比如用pid, 但是有时候这个进程还没启动, 比如task的comm来判定, 但是kernel中是不支持strcmp来判断字符串是否相等, 因为需要跑函数 g
10739 0
gdb kernel debug的进程断点
gdb调试kernel的时候, 如果设置通用函数断点, 比如vfs_read, 就会遇到一堆撞到断点的地方, 比如tty输入一个字符, 就是vfs_read, 没办法调试具体的某一个进程 一种办法就是条件断点, 其实不是很好用, 比如用pid, 但是有时候这个进程还没启动, 比如task的comm来判定, 但是kernel中是不支持strcmp来判断字符串是否相等, 因为需要跑函数 g
1140 0
基于qml创建最简单的图像处理程序(2)-使用c++&qml进行图像处理
《基于qml创建最简单的图像处理程序》系列课程及配套代码基于qml创建最简单的图像处理程序(1)-基于qml创建界面http://www.cnblogs.com/jsxyhelu/p/8343310.
1163 0
海康萤石摄像机远程监控机制分析
版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/voidreturn/article/details/76349814 背景介绍 有一个远程监控的需求,正常情况下采用服务器转发视频流对服务器要求很高,所以p2p看似是一个比较好的方案,但是p2p打洞却不是一件容易的事情,所以打算分析下海康萤石摄像机远程监控的机制,是否采用了p2p的方式。
1856 0
Golang 编写的图片压缩程序,质量、尺寸压缩,批量、单张压缩
目录:   前序   效果图   简介   全部代码   前序:   接触 golang 不久,一直是边学边做,边总结,深深感到这门语言的魅力,等下要跟大家分享是最近项目 服务端 用到的图片压缩程序,我单独分离了出来,做成了 exe 程序,可以在 Window 下运行。
1794 0
10个有用的处理图像的jQuery插件
我们都知道jquery具有灵活性和许多用途,尤其是在一些令人印象深刻的炫彩效果中都必须用jquery的应用图像才能做出来,他真的可以做一些很有用的东西。效果足以让你惊奇,在本文中我推荐10个jquery插件来帮助你处理图像 幻灯片 幻灯片是一个疯狂的简单的jQuery幻灯片插件。
711 0
3402
文章
0
问答
文章排行榜
最热
最新
相关电子书
更多
低代码开发师(初级)实战教程
立即下载
阿里巴巴DevOps 最佳实践手册
立即下载
冬季实战营第三期:MySQL数据库进阶实战
立即下载