• GDAL切割重采样遥感图像


    一个小测试程序开发全过程实录,完全新手入门级的实例,如果你还在为处理大影像而发愁,来试试这个称手的工具吧。

    Imagec 开发日记
    2013-6-25


    需求:


    影像数据切割,重采样
    数据切割的要求是简单的给予矩形的等分切割,并以2的幂次为分割单元,无需使用AOI裁切,
    重采样需要实现多种采样模式,用户可以切换采样模式(下文中所提供的代码只是利用了RasterIO的一个特性使用了默认的最近邻重采样方法)


    基本思路


    考虑是否存在使实用的gdal接口,
    自行设计,利用GDAL的读写接口完成数据输入输出工作
    初步了解,GDAL并不提供现成的切割和重采样的接口。
    /////////////////////////////////////////////////////////////////////
    2013-6-26


    确定使用自己的方法实现:

    裁切


    1、读入数据
    2、设置剪切范围
    3、完成裁切,将数据写入目标文件

    2013-6-27


    重采样

    1.读入数据
    2.借鉴Warp工具代码重写重采样代码
    将数据写入目标文件

    2013-6-28


    重采样2


    1.读入数据
    2.使用RasterIO自带的重采样功能
    3.数据写入目标文件

    ////////////////////////////////////////////////////////////////////////////
    2013-6-27


    实现算法


    裁切:


    读入文件
    GDALDataset *poDataset;
    GDALAllRegister();

    //open a dataset
    poDataset =(GDALDataset *)GDALOpen(fn,GA_ReadOnly);

    if(poDataset == NULL)
    设置相关参数为实现裁切做准备
    char **papsOptions = NULL;
    const int iXSize = poSrcDS->GetRasterXSize()/2;
    const int iYSize = poSrcDS->GetRasterYSize()/2;
    int iSize = GDALGetBandTypeSize(GDT_Byte)/8;
    //生成一个用于存放数据的缓存空间
    poDstDS = poDriver->Create(fnDst,iXSize,iYSize,3,GDT_Byte,papsOptions);
    GByte *abyByte = new GByte[iXSize*iYSize*bandCount];
    int bandMap[3] = {1,2,3};

    const char *pszSRS_WKT = poSrcDS->GetProjectionRef();
    poDstDS->SetProjection(pszSRS_WKT);
    poDstDS->SetGeoTransform(adfGeoTransform);

    //进行裁切工作,读取原始的数据,写入目标文件
    poSrcDS->RasterIO(GF_Read,0,0,iXSize,iYSize,abyRater,
    iXSize,iYSize,GDT_Byte,bandCount,bandMap,iSize*bandCount,iSize*iXSize*3,iSize);
    poDstDS->RasterIO(GF_Write,0,0,iXSize,iYSize,abyRater,
    iXSize,iYSize,GDT_Byte,bandCount,bandMap,iSize*bandCount,iSize*iXSize*3,iSize);
    ////////////////////////////////////////////////////////////////


    //数据重采样实现:


    打开数据
    GDALDataset *SrcDS =  (GDALDataset *)GDALOpen(fn,GA_ReadOnly);

    为RasterIO设置一些必要参数:

    const char *strWkt =pSrcDS->GetProjectionRef();
    int bandCount = pSrcDS->GetRasterCount();
    int iSize = GDALGetDataTypeSize(GDT_Byte)/8;
    int bandMap = {1,2,3};

    double adfGeoTransform[6];
    pSrcDS->GetGeoTransformation(adfGeoTransform);
    double dProj[4] ={0};
    ImageRowCol2Projection(adfGeoTransform,0,0,dProj[0],dProj[2]);//具体参考李明录重采样博文
    ImageRowCol2Projection(adfGeoTransform,pSrcDS->GetRasterXSize(),pSrcDS->GetRasterYSize(),dProj[1],dProj[3]);
    double maxX = dProj[0] >dProj[1]?dProj[0]:dProj[1];
    double minX = dProj[0]<dProj[1]?dProj[1]:dProj[0];
    double maxY = dProj[2]>dProj[3]?dProj[2]:dProj[3];
    double minY = dProj[2]<dProj[3]?dProj[2]:dProj[3];

    double fResX = 2;
    double fResY = 2;
    adfGeoTransform[0] = maxX;
    adfGeoTransform[3] = maxY;
    adfGeoTransform[1] = adfGeoTransform[1]/fResX;
    adfGeoTransform[4] = adfGeoTransform[5]/fResY;

    //得到输出图像长宽
    int iNewWidth = static_cast<int>((maxX-minX)/ABS(adfGeoTransform[1]) + 0.5);
    int iNewHeight = static_cast<int>((maxY-minY)/ABS(adfGeoTransform[5]) + 0.5);
    GByte *abyraster2= new GByte[iNewWidth*iNewHeight*bandCount];
    创建输出文件:
    GDALDataset *pDstDS2 = poDriver->Create("good.img",
    iNewWidth,iNewHeight,3,GDT_Byte,NULL);

    设置输出数据的投影:
    pDstDS2->SetGeoTransform(adfGeoTransform);
    pDstDS2->SetProjection(strWkt);
    使用RasterIO实现数据读写:(重采样全在这读写之间,参考第一篇文献的 第四种方式 一节)
    pSrcDS->RasterIO(GF_Read,0,0,pSrcDS->GetRasterXSize(),pSrcDS->GetRasterYSize(),abyraster2,
    iNewWidth,iNewHeight,GDT_Byte,bandCount,bandMap,iSize*bandCount,
    iSize*iNewWidth*3,iSize);
    pDstDS->RasterIO(GF_Write,0,0,iNewWidth,iNewHeight,GDT_Byte,bandCount,bandMap,iSize*bandCount,
    iSize*iNewWidth*3,iSize);
    完成读写后关闭打开的数据集,释放占用的内存
    if(pDstDS != NULL)
    {
    GDALClose((GDALDatasetH)pDstDS2);
    }
    GDALClose((GDALDatasetH)pSrcDS);

    报错实例


    没有注册驱动导致数据读取失败
    GDALAllRegister();
    ERROR 4'../image/result.img' not recognized as a supported format.


    投影问题
    setGeoTransform 无法为输出图像计算一个投影

    参考文献

    GDAL源码剖析(七)之GDAL RasterIO使用说明

    如何使用GDAL重采样图像

    如何使用GDAL进行AOI裁剪

  • 相关阅读:
    Linux菜鸟级重点
    在与 SQL Server 建立连接时出现与网络相关的或特定于实例的错误
    搭建PHP开发环境
    Struts+Hibernate+Spring实现用户登录功能
    Struts2整合Hibernate3实现用户登录功能
    决战JS(二)
    决战JS
    lightoj-1098
    lightoj-1072
    lightoj-1094 Farthest Nodes in a Tree(求树的直径)
  • 原文地址:https://www.cnblogs.com/dyllove98/p/3181737.html
Copyright © 2020-2023  润新知