• GDAL源码剖析(十)之编写自己的扩展格式


    一、简介

    该节内容参考GDAL的英文原文:http://www.gdal.org/gdal_drivertut.html。

    通常,可以通过从GDALDataset和GDALRasterBand继承来实现GDAL对新的数据格式支持。同时,还需要为这种格式创建一个GDALDriver的实例,让后通过GDALDriverManager将该新的驱动注册给GDAL系统。

    该教程将为JDEM数据格式实现一个简单的只读驱动开始,进而使用RawRasterBand类,实现一个可创建和修改图像数据的格式以及更高级的问题。

    强烈建议在看这篇文章之前,仔细理解GDAL数据模型的相关内容。地址为http://www.gdal.org /gdal_datamodel.html,或者参考上篇博文,GDAL体系架构。

    该文章主要包含下面9部分内容:

    • 从Dataset继承
    • 从RasterBand继承
    • 栅格驱动(Driver)
    • 添加驱动到GDAL中
    • 添加地理参考信息
    • 金字塔(快视图)
    • 创建文件
    • RawDataset和RawRasterBand类
    • 元数据和其他外部扩展

    二、从Dataset继承

    从这里开始,我们将演示如何创建一个日本DEM数据的只读驱动。代码地址为:http:// www.gdal.org/jdemdataset.cpp.html。首先,从GDALDataset继承一个子类JDEMDataset。代码如下:

    class JDEMRasterBand;
    class JDEMDataset : publicGDALPamDataset
    {
        friend class JDEMRasterBand;
        VSILFILE    *fp;
        GByte   abyHeader[1012];
      public:
                        JDEMDataset();
                       ~JDEMDataset();
        static GDALDataset*Open( GDALOpenInfo* );
        static int Identify( GDALOpenInfo* );
        CPLErr GetGeoTransform( double* padfTransform );
        const char *GetProjectionRef();
    };

    我们可以通过重载基类GDALDataset中的一些虚函数来为驱动重新实现这些特殊的功能。但是,Open()相对比较特殊,它不是基类 GDALDataset的虚函数。我们需要一个独立的函数来实现这个功能,因此我们把他声明为静态的(static)。将Open()声明为JDEMDataset的静态函数可以方便的访问JDEMDataset的私有方法去修改内容。

    Open()函数具体实现如下:

    GDALDataset *JDEMDataset::Open(GDALOpenInfo * poOpenInfo)
    {
        if (!Identify(poOpenInfo))
            return NULL;
     
    /*-------------------------------------------------------------------- */
    /*     Confirm the requested access is supported.                      */
    /* --------------------------------------------------------------------*/
        if( poOpenInfo->eAccess == GA_Update)
        {
            CPLError( CE_Failure,CPLE_NotSupported,
                      "The JDEM driver does not support update access toexisting"
                      " datasets.\n" );
            return NULL;
        }
     
    /*-------------------------------------------------------------------- */
    /*     Create a corresponding GDALDataset.                             */
    /* --------------------------------------------------------------------*/
        JDEMDataset     *poDS;
     
        poDS = new JDEMDataset();
     
        poDS->fp = VSIFOpenL( poOpenInfo->pszFilename, "rb");
        if (poDS->fp == NULL)
        {
            delete poDS;
            return NULL;
        }
       
    /*-------------------------------------------------------------------- */
    /*     Read the header.                                               */
    /*-------------------------------------------------------------------- */
        VSIFReadL( poDS->abyHeader, 1, 1012, poDS->fp );
     
        poDS->nRasterXSize= JDEMGetField( (char*) poDS->abyHeader+ 23, 3 );
        poDS->nRasterYSize= JDEMGetField( (char*) poDS->abyHeader+ 26, 3 );
        if  (poDS->nRasterXSize<= 0 || poDS->nRasterYSize<= 0 )
        {
            CPLError( CE_Failure,CPLE_AppDefined,
                      "Invalid dimensions : %d x %d",
                      poDS->nRasterXSize,poDS->nRasterYSize);
            delete poDS;
            return NULL;
        }
     
    /* --------------------------------------------------------------------*/
    /*     Create band information objects.                                */
    /*-------------------------------------------------------------------- */
        poDS->SetBand(1, new JDEMRasterBand(poDS, 1 ));
     
    /*-------------------------------------------------------------------- */
    /*     Initialize any PAM information.                                 */
    /*-------------------------------------------------------------------- */
        poDS->SetDescription(poOpenInfo->pszFilename);
        poDS->TryLoadXML();
     
    /*-------------------------------------------------------------------- */
    /*     Check for overviews.                                            */
    /* --------------------------------------------------------------------*/
        poDS->oOvManager.Initialize( poDS,poOpenInfo->pszFilename);
     
        return( poDS );
    }

    第一步,打开文件,判断驱动是否支持该类型。这里我们应该知道,GDAL在打开文件的时候会调用每个注册驱动的Open函数将文件打开,直到有一个成功为止,可以参考GDAL源代码中的GDALOpen函数实现部分。如果指定的文件不是驱动所支持的类型,则它们必须返回NULL。如果格式是它们所支持的类型,但是数据已经被破坏,那么它们应该产生一个错误。

    在文件打开之后文件的信息会保存在GDALOpenInfo对象中。GDALOpenInfo对象的公有成员如下:

    char        *pszFilename;
    char        **papszSiblingFiles;
    GDALAccess  eAccess;
    int         bStatOK;
    int         bIsDirectory;
    FILE        *fp;
    int         nHeaderBytes;
    GByte       *pabyHeader;
    驱动可以通过检查这些值来判断是否支持这个格式的文件。如果pszFilename引用的是系统中的一个文件,那么bStatOK将被设置为TRUE,同时,如果文件打开成功,程序会读取前一千个字节的数据,并且将数据存储在pabyHeader中,nHeaderBytes保存实际读取的字节数。

    在这个例子中,假设文件已经被打开并且已经读取到文件头部的一些信息。在这里,JDEM并没有魔术数字(magic numbers,这个没搞懂,个人理解是一个文件标识数字),因此我们只能检查不同的数据域是否合理。如果检查文件不是当前驱动所支持的格式,那么将返回NULL。检查函数的代码如下:

    int JDEMDataset::Identify(GDALOpenInfo * poOpenInfo)
    {
    /*-------------------------------------------------------------------- */
    /*     Confirm that the header has what appears to be dates in the     */
    /*     expected locations.  Sadly this isa relatively weak test.      */
    /*-------------------------------------------------------------------- */
        if( poOpenInfo->nHeaderBytes < 50 )
            return FALSE;
     
        /* check if century values seem reasonable */
        if( (!EQUALN((char *)poOpenInfo->pabyHeader+11,"19",2)
             && !EQUALN((char *)poOpenInfo->pabyHeader+11,"20",2))
            || (!EQUALN((char *)poOpenInfo->pabyHeader+15,"19",2)
                 &&!EQUALN((char*)poOpenInfo->pabyHeader+15,"20",2))
            || (!EQUALN((char *)poOpenInfo->pabyHeader+19,"19",2)
                && !EQUALN((char *)poOpenInfo->pabyHeader+19,"20",2)))
        {
            return FALSE;
        }
     
        return TRUE;
    }

    在真正的监测中,检查代码越全面越好。这里的监测相对比较弱。如果一个其他文件的在相应的位置具有相同的内容,那么它们很可能被错误地当作JDEM格式处理,最终导致不可预期的结果。

    对于满足要求的文件格式,还有另外一个必须做的检查,就是检查文件是否可写,对于JDEM格式,需要进行判断,代码如下:

    if( poOpenInfo->eAccess== GA_Update )
    {
        CPLError( CE_Failure,CPLE_NotSupported,
                  "The JDEMdriver does not support update access to existing"
                  "datasets.\n" );
        return NULL;
    }

    如果文件是我们支持的类型,满足上面的监测要求,我们需要创建一个数据集的实例,并在里面记录一些必要的信息。

    JDEMDataset     *poDS;
    poDS = new JDEMDataset();
    poDS->fp = VSIFOpenL(poOpenInfo->pszFilename,"rb" );
    if (poDS->fp == NULL)
    {
        delete poDS;
        return NULL;
    }

    通常在这时我们已经将打开文件了,并且保存文件的指针。在接下来,尽可能的使用GDAL提供的接口中VSI*L的函数来访问文件,这些POSIX风格的API可以访问大文件,内存文件以及Zip压缩的文件。

    接下来是要从读取到的文件头信息中提取图像的宽和高,nRasterXSize和nRasterYSize是在基类GDALDataset中定义的,这两个值必须在Open函数中进行赋值。代码如下:

    VSIFReadL( poDS->abyHeader,1, 1012, poDS->fp);
    poDS->nRasterXSize = JDEMGetField((char *) poDS->abyHeader + 23, 3 );
    poDS->nRasterYSize = JDEMGetField((char *) poDS->abyHeader + 26, 3 );
    if  (poDS->nRasterXSize <= 0 || poDS->nRasterYSize <= 0 )
    {
        CPLError( CE_Failure,CPLE_AppDefined,
                  "Invalid dimensions : %d x %d",
                  poDS->nRasterXSize,poDS->nRasterYSize);
        delete poDS;
        return NULL;
    }
    通过SetBand()函数将所有的波段与当前的GDALDataset对象绑定。在下一节,我们将讨论JDEMRasterBand类的具体细节。

    poDS->SetBand( 1, new JDEMRasterBand( poDS,1 ));

    最后,我们需要给dataset指定一个名称,并且调用GDALPamDataset类的函数TryLoadXML()重.aux.xml文件中读取一些辅助信息,更详细的内容请参阅GDALPamDataset类的相关函数。

    poDS->SetDescription( poOpenInfo->pszFilename );
    poDS->TryLoadXML();
    poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename );
    return( poDS );

    三、从RasterBand继承

    和从GDALDataset类派生子类JDEMDataset一样,我们需要从GDALRasterBand类派生一个JDEMRasterBand类,JDEMRasterBand类的定义如下:

    class JDEMRasterBand : publicGDALPamRasterBand
    {
        friend class JDEMDataset;
     
        int          nRecordSize;
        char*        pszRecord;
       
      public:
     
                JDEMRasterBand(JDEMDataset *, int);
                    ~JDEMRasterBand();
       
        virtual CPLErr IReadBlock( int, int, void * );
    };

    构造函数可以任意定义,这里不写使用默认的就行,但是只能从Open()函数中调用。其他的虚函数必须和gdal_priv.h中定义的一致,例如IReadBlock()。构造函数实现代码如下:

    JDEMRasterBand::JDEMRasterBand( JDEMDataset*poDS, int nBand )
    {
        this->poDS = poDS;
        this->nBand = nBand;
     
        eDataType = GDT_Float32;
     
        nBlockXSize = poDS->GetRasterXSize();
        nBlockYSize = 1;
     
        /* Cannot overflow as nBlockXSize <= 999 */
        nRecordSize = nBlockXSize*5+ 9 + 2;
        pszRecord = NULL;
    }

    下面列举的成员变量在GDALReasterBand中定义的,但是需要在子类的构造函数中进行初始化:

    • poDS: 指向基类GDALDataset的指针。
    • nBand: 对应数据集中的波段序号。
    • eDataType: 当前波段的数据类型。
    • nBlockXSize: 当前波段的块宽度。
    • nBlockYSize: 当前波段的块高度。

    GDALDataType类型的定义在文件gdal.h中,包括GDT_Byte、GDT_UInt16、GDT_Int16和 GDT_Float32。块的大小记录数据的实际或有效的大小。对于分块数据集这个值将是一个分块的大小,而对于其他大多数数据而言这个值将是一行。

    接下来,我们将实现真正读取影像数据的代码,IReadBlock()函数。

    CPLErr JDEMRasterBand::IReadBlock(int nBlockXOff,int nBlockYOff,
                                      void * pImage )
    {
        JDEMDataset *poGDS= (JDEMDataset *) poDS;
        int     i;
       
        if (pszRecord == NULL)
        {
            if (nRecordSize< 0)
                return CE_Failure;
     
            pszRecord = (char*) VSIMalloc(nRecordSize);
            if (pszRecord == NULL)
            {
                CPLError(CE_Failure,CPLE_OutOfMemory,
                        "Cannot allocate scanline buffer");
                nRecordSize = -1;
                return CE_Failure;
            }
        }
     
        VSIFSeekL( poGDS->fp, 1011 + nRecordSize*nBlockYOff, SEEK_SET);
     
        VSIFReadL( pszRecord,1, nRecordSize, poGDS->fp );
     
        if( !EQUALN((char *) poGDS->abyHeader,pszRecord,6))
        {
            CPLError( CE_Failure,CPLE_AppDefined,
                      "JDEM Scanline corrupt.  Perhaps file was not transferred\n"
                      "in binary mode?" );
            return CE_Failure;
        }
       
        if( JDEMGetField( pszRecord + 6, 3 ) != nBlockYOff+ 1 )
        {
            CPLError( CE_Failure,CPLE_AppDefined,
                      "JDEM scanline out of order, JDEM driver doesnot\n"
                      "currently support partial datasets." );
            return CE_Failure;
        }
     
        for( i = 0; i < nBlockXSize;i++ )
            ((float *) pImage)[i] = (float)
                (JDEMGetField( pszRecord+ 9 + 5 * i, 5) * 0.1);
     
        return CE_None;
    }

    在这里需要注意的有:

    • 需要将GDALRasterBand::poDS成员变量传递给子类,如果你的RasterBand类需要访问数据集的私有成员,那么需要吧RasterBand类声明为Dataset的友元类。
    • 如果在处理中出现错误,请使用CPLError()提示,并且返回CE_Failure,否则返回 CE_None即可。
    • pImage缓存应该用一个块的数据来进行赋值,块大的小nBlockXSize 和nBlockYSize在波段中定义。pImage的数据类型应该和栅格波段对象中声明的数据类型eDataType一致。
    • nBlockXOff和nBlockYOff是块的偏移量,因此对于一个128x128大小的块,序号为1和1的块应该对应原始数据是从(128,128)到(255,255)的数据。

    四、栅格驱动(Driver)

    虽然JDEMDataset和JDEMRasterBand已经可以读数据,但是GDAL库仍然不清楚JDEMDataset的任何信息。这可以通过GDALDriverManager来完成。为了将我们自己实现的驱动注册到GDAL库中,我们需要重新实现一个注册函数,函数定义写在gdal_frmts.h文件中:

    CPL_C_START
    void CPL_DLL GDALRegister_JDEM(void);
    CPL_C_END

    函数的实现部分如下:

    void GDALRegister_JDEM()
    {
        GDALDriver  *poDriver;
     
        if( GDALGetDriverByName("JDEM" ) == NULL )
        {
            poDriver = new GDALDriver();
           
            poDriver->SetDescription("JDEM" );
            poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
                                       "Japanese DEM (.mem)" );
            poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
                                       "frmt_various.html#JDEM" );
            poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "mem");
            poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES" );
     
            poDriver->pfnOpen= JDEMDataset::Open;
            poDriver->pfnIdentify= JDEMDataset::Identify;
     
            GetGDALDriverManager()->RegisterDriver( poDriver );
        }
    }

    可以使用GDAL_CHECK_VERSION宏来判断GDAL的版本(从1.5.0版本之后才有),这个宏可以用来使GDAL扩展一些额外的库,比如编写一个GDAL的插件驱动。对于这种情况,尽可能使用这个宏来判断GDAL的版本信息,对于一些老的版本,GDAL可能会不支持。

    当第一次被调用的时候,注册函数将创建一个GDALDriver的对象,并且通过GDALDriverManager来注册。在用GDALDriverManager注册驱动之前,需要设置以下的成员变量:

    驱动对应格式的名称,不能和其他驱动名字冲突。通常在3-5个字符长度,并且匹配驱动类名字的前缀。(手工设置)

    • GDAL_DMD_LONGNAME: 文件格式的详细描述,长度一般在50-60个字符。(手工设置)
    • GDAL_DMD_HELPTOPIC: 如果有的话,是关于这个驱动帮助文档的名称。在这个例子中,JDEM格式包含在frmt_various.html#JDEM位置描述。(可选)
    • GDAL_DMD_EXTENSION: 该类型文件的扩展名。如果扩展名多于一个,则最好选择最通用的那个,或者直接为空。(可选)
    • GDAL_DMD_MIMETYPE: 该格式数据的标准用途类型,例如“image/png”。(可选)
    • GDAL_DMD_CREATIONOPTIONLIST: 用于描述创建时的选项。可以参考geotiff驱动的实现代码。(可选)
    • GDAL_DMD_CREATIONDATATYPES: 支持创建数据集的全部类型列表。如果存在Create()方法这些类型都将被支持。如果存在CreateCopy()方法,该类型列表将是那些支持无损输出的类型。比如,一个格式使用了CreateCopy()方法,如果可以写为Float32类型,那么Byte、Int16和UInt16都应该被支持(因为它们都可以用Float32表示)。
    • pfnOpen: 打开这种格式文件的函数指针。(可选)
    • pfnCreate: 创建这中格式的函数指针。(可选)
    • pfnCreateCopy: 创建一个从其它数据源复制而来的新数据集,但是不需要更新的函数。(可选)
    • pfnDelete: 删除这种格式数据集函数指针。(可选)
    • pfnUnloadDriver: 这个函数只有在驱动被销毁的时候才被调用。在驱动层被用来清除数据。很少用到。(可选)

    五、添加驱动到GDAL库

    GDALRegister_JDEM()函数必须在上层的程序中调用,才能使支持JDEM格式的数据。通常需要以下的工作来完成:

    1. 在gdal/frmts 创建一个目录,目录的名字与驱动的短名称相同。
    2. 在目录中添加文件GNUmakefile和makefile.vc,具体格式可以参考其他目录中,如jdem目录中的文件。
    3. 添加dataset和RasterBand的实现部分(这里值得是源代码)。通常命名规则值<短名称> dataset.cpp,如jdemdataset.cpp。
    4. 在gdal/gcore/gdal_frmts.h文件中添加一个注册函数声明,如GDALRegister_JDEM())。
    5. 在frmts/gdalallregister.cpp文件中的注册函数中调用该注册函数,可以使用ifdef来进行判断,可以参考相关代码。
    6. 在文件GDALmake.opt.in (和GDALmake.opt)中的GDAL_FORMATS宏中添加格式的短名称。
    7. 在文件frmts/makefile.vc中的EXTRAFLAGS宏添加该项,具体参考相关代码。

    完成上面的步骤,重新编译生成GDAL库后,那么所有的工具都会支持这种新的文件格式,可以使用gdalinfo工具来进行测试打开数据和显示数据的信息,可以使用gdal_translate工具测试图像的读取。

    六、添加地理参考信息

    接下来我们继续完善之前的工作,这里将使驱动支持地理参考信息,我们需要重写JDEMDataset类中的两个虚函数,这两个函数的定义在GDALRasterDataset基类中。两个函数定义如下:

    virtual const char *GetProjectionRef(void);
    virtual CPLErr GetGeoTransform(double * );

    在JDEMDataset类中添加函数声明:

    CPLErr GetGeoTransform(double * padfTransform);
    const char *GetProjectionRef();

    重新实现的GetGeoTransform()函数只是复制地理转换矩阵到缓存中。GetGeoTransform()函数可能经常使用,因此一般最好使该函数体小。在很多时候,使用Open()函数已经获取到地理转换的信息,然后用这个函数进行复制里面的信息即可。需要注意的是,转换矩阵中的起始点坐标对应像素左上角的坐标,而不是左上角象元的中心坐标。

    CPLErr JDEMDataset::GetGeoTransform(double * padfTransform)
    {
        double  dfLLLat, dfLLLong,dfURLat, dfURLong;
     
        dfLLLat = JDEMGetAngle((char *) abyHeader+ 29 );
        dfLLLong = JDEMGetAngle((char *) abyHeader+ 36 );
        dfURLat = JDEMGetAngle((char *) abyHeader+ 43 );
        dfURLong = JDEMGetAngle((char *) abyHeader+ 50 );
       
        padfTransform[0] = dfLLLong;
        padfTransform[3] = dfURLat;
        padfTransform[1] = (dfURLong- dfLLLong) / GetRasterXSize();
        padfTransform[2] = 0.0;
           
        padfTransform[4] = 0.0;
        padfTransform[5] = -1 * (dfURLat- dfLLLat) / GetRasterYSize();
     
        return CE_None;
    }

    GetProjectionRef ()方法将返回一个字符串,使用的是OGC WKT格式的坐标系统的定义。在这个例子中,这种格式的所有数据只有一种坐标系统。但是在比较复杂的数据格式中,我们需要使用OGRSpatialReference来针对特定情形来进行导出坐标系统。

    const char *JDEMDataset::GetProjectionRef()
    {
        return( "GEOGCS[\"Tokyo\",DATUM[\"Tokyo\",SPHEROID[\"Bessel1841\",6377397.155,299.1528128,AUTHORITY[\"EPSG\",7004]],TOWGS84[-148,507,685,0,0,0,0],AUTHORITY[\"EPSG\",6301]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",8901]],UNIT[\"DMSH\",0.0174532925199433,AUTHORITY[\"EPSG\",9108]],AUTHORITY[\"EPSG\",4301]]");
    }

    到这里已经完成了JDEM驱动的代码,整个源代码可以在jdemdataset.cpp(地址为http://www.gdal.org/jdemdataset.cpp.html)文件中看到。

    七、金字塔(快视图)

    GDAL可以使用GDALRasterBand::GetOverview()和相关的函数来对文件建立金字塔。然后,这些内容已经超出了本文档的讨论范围,具体可以参考GeoTIFF驱动(gdal/frmts/gtiff/geotiff.cpp)中的相关代码。

    通过GDALRasterBand的函数HasArbitraryOverviews()可以获取文件中是否含有金字塔信息,返回TRUE表示有。在这个例子中,栅格波段对象应该重载实现RasterIO()函数,通过该函数可以有效的访问图像或者重采样等。这个函数比较复杂,需要判断很多参数是否正确等,关于这个的例子可以在OGDI和ECW格式中找到相关的代码。

    然而,很多图像的金字塔都使用GDAL默认的外部金字塔格式来保存,实际上是一个TIFF文件,但是后缀名为.ovr,为了使用读取和创建金字塔,需要对GDALDataset类中的oOvManager进行初始化。通常情况下,在Open()函数的最后调用类似下面的代码(PAM TryLoadXML()函数之前)。

    poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename );

    这样依赖,该格式可以读取和创建金字塔。一般简单的图像格式,都可以使用这种方式,除非一些特殊的格式有自定义的金字塔格式。

    八、创建文件

    有两种方法来创建文件。第一种是调用CreateCopy()函数,可以从源影像中获取影像信息并写入到结果图像中;第二种是调用Create动态创建文件,包含Create函数和用来设置各种信息的函数。

    第一种方法的优点是在创建文件的所有的信息都从源图像中获取并写入到结果图像中。尤其是对于文件那些特别重要的信息例如颜色表、参照系等信息。关于CreateCopy方法的缺点是对于一些信息没有相应的设置方法,例如min/max、scaling、description和GCPs。

    第二种方法的优点在于可以创建一个空的新文件,并且在需要的时候把结果写入其中。对于一个影像来说,动态创建不需要提前获取影像的所有信息。

    对于非常重要的格式而言,两种方法最好都支持。

    1、CreateCopy

    GDALDriver::CreateCopy()方法可以直接调用,因此只需要知道调用的参数既可。不过需要注意以下的细节:

    • 如果参数bStrict为FALSE,那么驱动器将对数据自行做一些合适的处理。特别是对于那些并不完全等同格式。
    • 可以实现CreateCopy进度提示。回调函数的返回值表示是否需要继续进行,并且过程最好是在合理的时间单位中被调用(这个例子没有示范)。
    • 特殊的创建选项必须在帮助中说明。如果存在"NAME=VALUE"的格式,可以使用CPLFetchNameValue()函数来进行设置,比如JPEGCreateCopy()函数就有QUALITY和PROGRESSIVE两个标志来说明相关的创建选项。
    • 返回的GDALDataset指针是只读或者更新模式。在通常情况下返回更新模式,否则返回只读模式即可。

    JPEG格式对应的CreateCopy的代码如下:

    static GDALDataset *
    JPEGCreateCopy( const char * pszFilename, GDALDataset*poSrcDS,
                   int bStrict, char ** papszOptions,
                   GDALProgressFuncpfnProgress, void* pProgressData )
    {
        int  nBands = poSrcDS->GetRasterCount();
        int  nXSize = poSrcDS->GetRasterXSize();
        int  nYSize = poSrcDS->GetRasterYSize();
        int  nQuality = 75;
        int  bProgressive = FALSE;
     
        // --------------------------------------------------
        //      Some somerudimentary checks                  
        // --------------------------------------------------
        if( nBands != 1&& nBands != 3 )
        {
            CPLError( CE_Failure,CPLE_NotSupported,
                "JPEG driver doesn't support %d bands.  Must be 1 (grey) "
                "or 3 (RGB) bands.\n", nBands );
     
            return NULL;
        }
     
        if( poSrcDS->GetRasterBand(1)->GetRasterDataType()!= GDT_Byte && bStrict )
        {
            CPLError( CE_Failure,CPLE_NotSupported,
                "JPEG driver doesn't support data type %s. "
                "Only eight bit byte bands supported.\n",
                GDALGetDataTypeName(
                poSrcDS->GetRasterBand(1)->GetRasterDataType()) );
     
            return NULL;
        }
     
        // --------------------------------------------------
        //      What optionshas the user selected?                            
        // --------------------------------------------------
        if( CSLFetchNameValue(papszOptions,"QUALITY")!= NULL )
        {
            nQuality = atoi(CSLFetchNameValue(papszOptions,"QUALITY"));
            if( nQuality <10 || nQuality > 100 )
            {
                CPLError( CE_Failure,CPLE_IllegalArg,
                    "QUALITY=%s is not a legal value in the range10-100.",
                    CSLFetchNameValue(papszOptions,"QUALITY") );
                return NULL;
            }
        }
     
        if( CSLFetchNameValue(papszOptions,"PROGRESSIVE")!= NULL )
        {
            bProgressive = TRUE;
        }
     
        // --------------------------------------------------
        //      Create thedataset.                                            
        // --------------------------------------------------
        FILE    *fpImage;
     
        fpImage = VSIFOpen(pszFilename, "wb");
        if( fpImage == NULL )
        {
            CPLError( CE_Failure,CPLE_OpenFailed,
                "Unable to create jpeg file %s.\n",
                pszFilename );
            return NULL;
        }
     
        // --------------------------------------------------
        //      InitializeJPG access to the file.                             
        // --------------------------------------------------
        struct jpeg_compress_structsCInfo;
        struct jpeg_error_mgrsJErr;
     
        sCInfo.err = jpeg_std_error( &sJErr);
        jpeg_create_compress( &sCInfo );
     
        jpeg_stdio_dest( &sCInfo,fpImage );
     
        sCInfo.image_width= nXSize;
        sCInfo.image_height= nYSize;
        sCInfo.input_components= nBands;
     
        if( nBands == 1 )
        {
            sCInfo.in_color_space= JCS_GRAYSCALE;
        }
        else
        {
            sCInfo.in_color_space= JCS_RGB;
        }
     
        jpeg_set_defaults( &sCInfo);
     
        jpeg_set_quality( &sCInfo,nQuality, TRUE);
     
        if( bProgressive )
            jpeg_simple_progression( &sCInfo );
     
        jpeg_start_compress( &sCInfo,TRUE );
     
        // --------------------------------------------------
        //      Loop overimage, copying image data.                           
        // --------------------------------------------------
        GByte   *pabyScanline;
        CPLErr      eErr;
     
        pabyScanline = (GByte*) CPLMalloc( nBands* nXSize );
     
        for( int iLine = 0; iLine< nYSize; iLine++)
        {
            JSAMPLE     *ppSamples;
     
            for( int iBand = 0; iBand< nBands; iBand++)
            {
                GDALRasterBand * poBand= poSrcDS->GetRasterBand(iBand+1 );
                eErr = poBand->RasterIO( GF_Read,0, iLine, nXSize,1,
                    pabyScanline + iBand,nXSize, 1, GDT_Byte,
                    nBands, nBands* nXSize );
            }
     
            ppSamples = pabyScanline;
            jpeg_write_scanlines( &sCInfo, &ppSamples, 1 );
        }
     
        CPLFree( pabyScanline);
     
        jpeg_finish_compress( &sCInfo );
        jpeg_destroy_compress( &sCInfo );
     
        VSIFClose( fpImage);
     
        return (GDALDataset*) GDALOpen( pszFilename,GA_ReadOnly );
    }

    2、动态创建

    在动态创建的例子中,没有源数据集。但是提供了文件的大小、波段数和像素的数据类型。对于其他的一些信息(例如地理参照系等),将在以后通过特定的函数进行设置。

    下面使用的是一个简单实现了PCI.aux文件的创建示例。使用自己写的方法创建了一个空白的影像,并且在最后调用GDALOpen,使用更新方式打开数据,返回打开的数据集指针。

    GDALDataset *PAuxDataset::Create(const char * pszFilename,
                                     int nXSize, int nYSize, int nBands,
                                     GDALDataTypeeType,
                                     char ** // papszParmList  )
    {
        char    *pszAuxFilename;
     
        //--------------------------------------------------
        //      Verify inputoptions.                                          
        //--------------------------------------------------
        if( eType != GDT_Byte && eType!= GDT_Float32 && eType != GDT_UInt16
            && eType != GDT_Int16)
        {
            CPLError( CE_Failure,CPLE_AppDefined,
                "Attempt to create PCI .Aux labelled dataset with anillegal\n"
                "data type (%s).\n",
                GDALGetDataTypeName(eType));
     
            return NULL;
        }
     
        //--------------------------------------------------
        //      Try to createthe file.                                        
        //--------------------------------------------------
        FILE    *fp;
     
        fp = VSIFOpen( pszFilename, "w");
     
        if( fp == NULL )
        {
            CPLError( CE_Failure,CPLE_OpenFailed,
                "Attempt to create file `%s' failed.\n",
                pszFilename );
            return NULL;
        }
     
        //--------------------------------------------------
        //      Just writeout a couple of bytes to establish the binary       
        //      file, andthen close it.                                       
        //--------------------------------------------------
        VSIFWrite( (void*) "\0\0", 2, 1, fp );
        VSIFClose( fp);
     
        //--------------------------------------------------
        //      Create theaux filename.                                       
        //--------------------------------------------------
        pszAuxFilename = (char*) CPLMalloc(strlen(pszFilename)+5);
        strcpy( pszAuxFilename,pszFilename );;
     
        for( int i = strlen(pszAuxFilename)-1; i> 0; i-- )
        {
            if( pszAuxFilename[i] == '.' )
            {
                pszAuxFilename[i]= '\0';
                break;
            }
        }
     
        strcat( pszAuxFilename,".aux" );
     
        //--------------------------------------------------
        //      Open thefile.                                                 
        //--------------------------------------------------
        fp = VSIFOpen( pszAuxFilename, "wt");
        if( fp == NULL )
        {
            CPLError( CE_Failure,CPLE_OpenFailed,
                "Attempt to create file `%s' failed.\n",
                pszAuxFilename );
            return NULL;
        }
     
        //--------------------------------------------------
        //      We need towrite out the original filename but without any     
        //      pathcomponents in the AuxilaryTarget line. Do so now.        
        //--------------------------------------------------
        int     iStart;
     
        iStart = strlen(pszFilename)-1;
        while( iStart >0 && pszFilename[iStart-1] != '/'
            && pszFilename[iStart-1]!= '\\' )
            iStart--;
     
        VSIFPrintf( fp,"AuxilaryTarget: %s\n", pszFilename + iStart);
     
        //--------------------------------------------------
        //      Write out theraw definition for the dataset as a whole.       
        //--------------------------------------------------
        VSIFPrintf( fp,"RawDefinition: %d %d %d\n",
            nXSize, nYSize,nBands );
     
        //--------------------------------------------------
        //      Write out adefinition for each band.  We alwayswrite band    
        //      sequentialfiles for now as these are pretty efficiently       
        //      handled byGDAL.                                               
        //--------------------------------------------------
        int     nImgOffset = 0;
     
        for( int iBand = 0; iBand< nBands; iBand++)
        {
            const char * pszTypeName;
            int      nPixelOffset;
            int      nLineOffset;
     
            nPixelOffset = GDALGetDataTypeSize(eType)/8;
            nLineOffset = nXSize* nPixelOffset;
     
            if( eType == GDT_Float32 )
                pszTypeName = "32R";
            else if( eType == GDT_Int16)
                pszTypeName = "16S";
            else if( eType == GDT_UInt16)
                pszTypeName = "16U";
            else
                pszTypeName = "8U";
     
            VSIFPrintf( fp,"ChanDefinition-%d: %s %d %d %d %s\n",
                iBand+1, pszTypeName,
                nImgOffset, nPixelOffset,nLineOffset,
    #ifdef CPL_LSB
                "Swapped"
    #else
                "Unswapped"
    #endif
                );
     
            nImgOffset += nYSize* nLineOffset;
        }
     
        //--------------------------------------------------
        //      Cleanup                                                        
        //--------------------------------------------------
        VSIFClose( fp);
     
        return (GDALDataset*) GDALOpen( pszFilename,GA_Update );
    }

    如果文件格式支持动态创建或支持更新都需要实现GDALRasterBand中的IWriteBlock()方法,它类似于IReadBlock()。并且,在GDALRasterBand的析构函数中调用FlushCache()方法是比较危险的。因此,在析构方法调用之前必须确保任何保存在缓存块的数据都已经清除。

    九、RawDataset和RawRasterBand类

    很多文件格式将图像数据存储在一个使用二进制文件中,使用规则的按照行存储的格式。为了读取这些类似的格式,这里使用了一个RawDataset类和RawRasterBand类来进行处理,关于这两个类的定义可以在gdal/frmts/raw中找到。

    在这个例子中没有用到这个类,如果有用的,可以从RawRasterBand派生,数据集可以从RawDataset派生即可。

    Open()函数调用后通过构造函数实例化栅格波段里面的所有信息。可以查看PNM的驱动就是通过这样的方式来创建他的栅格波段的。

    RawRasterBand可以获取下面的参数:

    • poDS: 该波段的GDALDataset指针,这个数据集必须是从RawRasterDataset派生。
    • nBand: 波段序号,从1开始。
    • fpRaw: 图像数据打开后的文件句柄,类型是FILE *。
    • nImgOffset: 第一个象元的偏移量byte。
    • nPixelOffset: 第一个象元从在一行中的偏移量byte。
    • nLineOffset: 起始扫描线和下一个扫描线的偏移量byte。
    • eDataType: GDALDataType类型的数据类型,用来表示图像在磁盘上的存储类型。
    • bNativeOrder: (这句没看懂啥意思,原文放这里吧)FALSE if the data is not in the same endianness asthe machine GDAL is running on.个人感觉,这个的意思大概是如果GDAL读取的数据的顺序和存储的顺序不同的话,使用TRUE,然后GDAL会对这些数据自动进行交换。

    一个简单的裸数据读取的例子在gdal/frmts/raw目录中,里面有大量的关于裸数据的说明。

    十、元数据和其他外部扩展

    在GDAL数据模型中还有很多其他项目,在GDALDataset和GDALRasterBand中都存在对应的虚函数。它们包括:

    • Metadata: 关于数据集或者波段的Name/value。GDALMajorObject(包括其子类)都能支持元数据,可以在Open()函数中通过调用SetMetadataItem()函数来进行设置。SAR_CEOS(frmts/ceos2/sar_ceosdataset.cpp)和GeoTIFF驱动中都有关于读取元数据的例子。
    • ColorTables: GDT_Byte类型的波段可以包含与它们相关联的颜色表信息。在文件frmts/png/pngdataset.cpp中就是一个关于支持颜色表的例子。
    • ColorInterpretation: PNG驱动中包括一个返回波段被表示为红,绿,蓝,透明度或者灰度的描述信息。
    • GCPs: GDALDataset含有一系列的地面控制点与他们相对应的空间参考坐标进行关联(通过GetGeoTransform来映射转换)。MFF2格式(gdal/frmts.raw.hkvdataset.cpp)就是一个简单的支持GCP的例子。
    • NoDataValue: GetNoDataValue()可以获取波段的NODATA值。具体可以参考frmts/raw/pauxdataset.cpp驱动的相关代码。
    • Category Names: GetCategoryNames()函数可以根据影象的名字将其分类。不过目前还没有哪个驱动用到这个信息。

  • 相关阅读:
    Linux内核基础--事件通知链(notifier chain)good【转】
    10 个迅速提升你 Git 水平的提示【转】
    notifier chain — 内核通知链【转】
    内核通知链 学习笔记 【转】
    Linux内核基础--事件通知链(notifier chain)【转】
    Git 使用规范流程【转】
    Linux中断(interrupt)子系统之二:arch相关的硬件封装层【转】
    学习 Linux,101: 自定义或编写简单脚本【转】
    MySQL数据处理函数
    Effective JavaScript Item 36 实例状态仅仅保存在实例对象上
  • 原文地址:https://www.cnblogs.com/xiaowangba/p/6314017.html
Copyright © 2020-2023  润新知