GDAL库初学

最新推荐文章于 2024-05-25 01:54:49 发布

cocolalala

最新推荐文章于 2024-05-25 01:54:49 发布

阅读量1.1k

点赞数

分类专栏： Gis 文章标签： gdal Gis

Gis 专栏收录该内容

1 篇文章 0 订阅

订阅专栏

打开文件

在打开一个GDAL支持的栅格资料之前，必需要注册驱动。每个驱动对应各自支持的格式。

通常这个会被GDALAllRegister（）函数完成。

C++

#include"gdal_priv.h"
#include"cpl_conv.h" // forCPLMalloc()
 
int main(){
    GDALDataset  *poDataset;
    GDALAllRegister();
    poDataset = (GDALDataset*) GDALOpen( pszFilename, GA_ReadOnly );
    if( poDataset == NULL )
    {
        ...;
    }

Python

    import gdal
    from gdalconst import *
 
    dataset = gdal.Open( filename, GA_ReadOnly )
    if dataset is None:
        ...

获取数据集信息

// adfGeoTransform[0] /* top left x，图像左上角x坐标值*/

// adfGeoTransform[1] /* w-e pixel resolution，图像横坐标？米/每像素 */

// adfGeoTransform[2] /* rotation, 0 if image is "north up" */

// adfGeoTransform[3] /* top left y 图像左上角y坐标值*/

// adfGeoTransform[4] /* rotation, 0 if image is "north up" */

// adfGeoTransform[5] /* n-s pixel resolution 图像纵坐标？米/每像素*/

//由像素计算地理坐标值：

// Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2)

// Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5)

通过GDAL存取栅格数据，每次只完成一个波段

以波段为基础的元数据，区块，颜色表和其他可见的信息

dataset通过GetRasterCount()方法获取GDALRasterBand对象显示相关信息

C++：

        GDALRasterBand  *poBand;
        int             nBlockXSize, nBlockYSize;
        int             bGotMin, bGotMax;
        double          adfMinMax[2];
       
        poBand = poDataset->GetRasterBand( 1);
        poBand->GetBlockSize(&nBlockXSize, &nBlockYSize );
        printf( "Block=%dx%d Type=%s,ColorInterp=%s\n",
                nBlockXSize, nBlockYSize,
               GDALGetDataTypeName(poBand->GetRasterDataType()),
                GDALGetColorInterpretationName(
                   poBand->GetColorInterpretation()) );
 
        adfMinMax[0] = poBand->GetMinimum(&bGotMin );
        adfMinMax[1] = poBand->GetMaximum(&bGotMax );
        if( ! (bGotMin && bGotMax) )
           GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);
 
        printf( "Min=%.3fd,Max=%.3f\n", adfMinMax[0], adfMinMax[1] );
       
        if( poBand->GetOverviewCount() >0 )
            printf( "Band has %doverviews.\n", poBand->GetOverviewCount() );
 
        if( poBand->GetColorTable() != NULL)
            printf( "Band has a colortable with %d entries.\n",
                    poBand->GetColorTable()->GetColorEntryCount() );

python:

        print 'Driver: ',dataset.GetDriver().ShortName,'/', \
          dataset.GetDriver().LongName
        print'Size is ',dataset.RasterXSize,'x',dataset.RasterYSize, \
          'x',dataset.RasterCount
        print'Projection is ',dataset.GetProjection()
   
        geotransform= dataset.GetGeoTransform()
        if not geotransform is None:
             print 'Origin = (',geotransform[0],',',geotransform[3],')'
             print 'Pixel Size = (',geotransform[1],',',geotransform[5],')'

Reading raster data:读取栅格数据

有很多方法去读取栅格数据，但是最常见的是通过GDALRasterBand::RasterIO()方法。

这个方法将会自动的处理数据类型转换，向上/向下取样和构建窗口。

以下代码将会读取第一行数据到相同大小的缓冲区中，作为操作的一部分将其转换为浮点型。

C++：

        float *pafScanline;
        int  nXSize = poBand->GetXSize();
 
        pafScanline = (float *)CPLMalloc(sizeof(float)*nXSize);
        poBand->RasterIO( GF_Read, 0, 0,nXSize, 1,
                          pafScanline, nXSize,1, GDT_Float32,
                          0, 0 );

注：padScanline不再使用时用CPLFree()掉

python：

scanline = band.ReadRaster( 0, 0,band.XSize, 1, \
                                    band.XSize, 1, GDT_Float32 )

注：scanline是xsize*4字节的string floating pointer data，以下方法可以转换

        import struct

        tuple_of_floats = struct.unpack('f' * b2.XSize, scanline)

closing the dataset：关闭数据集

调用GDALClose()

Techniques for Creating Files: createCopy() 和 Creat()

CreateCopy 方法需要调用格式驱动上的CreateCopy（）方法，并且通过源dataset来复制

创建的方法需要调用驱动的Create（）方法，然后清楚的写入所有的元数据和栅格波段通过不同的调用

所有支持创建新文件的驱动都支持CreateCopy()方法，但是仅仅一部分支持Create()方法

若要确定一个格式驱动对象是否支持Create()或者CreateCopy()，我们可以考察元数据DCAP_CREATE 和 DCAP_CREATECOPY。确保在调用GetDriverByName()方法之前已经调用了GDALAllRegister()方法。在下面的例子中，循环一个驱动，来确定它是否支持Create()或者CreateCopy()。

c++：

 #include "cpl_string.h"
...
    const char *pszFormat = "GTiff";
    GDALDriver *poDriver;
    char **papszMetadata;
 
    poDriver =GetGDALDriverManager()->GetDriverByName(pszFormat);
 
    if( poDriver == NULL )
        exit( 1 );
 
    papszMetadata = poDriver->GetMetadata();
    if( CSLFetchBoolean( papszMetadata,GDAL_DCAP_CREATE, FALSE ) )
        printf( "Driver %s supports Create()method.\n", pszFormat );
    if( CSLFetchBoolean( papszMetadata,GDAL_DCAP_CREATECOPY, FALSE ) )
        printf( "Driver %s supportsCreateCopy() method.\n", pszFormat );

python：

    format = "GTiff"
    driver = gdal.GetDriverByName( format )
    metadata = driver.GetMetadata()
    if metadata.has_key(gdal.DCAP_CREATE) \
       and metadata[gdal.DCAP_CREATE] == 'YES':
        print 'Driver %s supports Create()method.' % format
    if metadata.has_key(gdal.DCAP_CREATECOPY) \
       and metadata[gdal.DCAP_CREATECOPY] =='YES':
        print 'Driver %s supports CreateCopy()method.' % format

Using CreateCopy()：

一个简单的拷贝，从一个名为pszSrcFilename的文件到一个新文件pszDstFilename,使用驱动是事先初始化好的默认选项格式拷贝方式

c++：

    GDALDataset *poSrcDS =
       (GDALDataset *) GDALOpen( pszSrcFilename,GA_ReadOnly );
    GDALDataset *poDstDS;
 
    poDstDS = poDriver->CreateCopy(pszDstFilename, poSrcDS, FALSE,
                                    NULL, NULL,NULL );
    /* Once we're done, close properly thedataset */
    if( poDstDS != NULL )
        GDALClose( (GDALDatasetH) poDstDS );
    GDALClose( (GDALDatasetH) poSrcDS );

python:

    src_ds = gdal.Open( src_filename )
    dst_ds = driver.CreateCopy( dst_filename,src_ds, 0 )
 
    # Once we're done, close properly thedataset
    dst_ds = None
    src_ds = None

包含自定义选项，使用预先定义的方式来传递选项值

c++:

#include "cpl_string.h"
    ...
    char **papszOptions = NULL;
   
    papszOptions = CSLSetNameValue(papszOptions, "TILED", "YES" );
    papszOptions = CSLSetNameValue(papszOptions, "COMPRESS", "PACKBITS" );
    poDstDS = poDriver->CreateCopy(pszDstFilename, poSrcDS, FALSE,
                                   papszOptions, GDALTermProgress, NULL );
 
    /* Once we're done, close properly thedataset */
    if( poDstDS != NULL )
        GDALClose( (GDALDatasetH) poDstDS );
    CSLDestroy( papszOptions );

python:

    src_ds = gdal.Open( src_filename )
    dst_ds = driver.CreateCopy( dst_filename,src_ds, 0,
                                [ 'TILED=YES','COMPRESS=PACKBITS' ] )
    # Once we're done, close properly thedataset
    dst_ds = None
    src_ds = None

Using Create()

有时候，你不仅输出一个现有文件到另一个新文件，一般来讲需要使用GDALDriver::Create()方法。该方法有跟CreateCopy()方法类似的属性选择列

表，但是图像大小，波段数和波段类型必须准确给出。

C++:

    GDALDataset *poDstDS;      
    char **papszOptions = NULL;
 
    poDstDS = poDriver->Create(pszDstFilename, 512, 512, 1, GDT_Byte,
                                papszOptions );

一旦数据集创建成功，所有的元数据和光栅数据必须被写入文件。这些数据的多样化依赖于个人设置习惯。这里是一个简单的例子。

    double adfGeoTransform[6] = { 444720, 30,0, 3751320, 0, -30 };
    OGRSpatialReference oSRS;
    char *pszSRS_WKT = NULL;
    GDALRasterBand *poBand;
    GByte abyRaster[512*512];
 
    poDstDS->SetGeoTransform( adfGeoTransform);
   
    oSRS.SetUTM( 11, TRUE );
    oSRS.SetWellKnownGeogCS( "NAD27");
    oSRS.exportToWkt( &pszSRS_WKT );
    poDstDS->SetProjection( pszSRS_WKT );
    CPLFree( pszSRS_WKT );
 
    poBand = poDstDS->GetRasterBand(1);
    poBand->RasterIO( GF_Write, 0, 0, 512,512,
                      abyRaster, 512, 512,GDT_Byte, 0, 0 );   
 
    /* Once we're done, close properly thedataset */
    GDALClose( (GDALDatasetH) poDstDS );

python:

    dst_ds = driver.Create( dst_filename, 512,512, 1, gdal.GDT_Byte )
 
    import osr
    import numpy
 
    dst_ds.SetGeoTransform( [ 444720, 30, 0,3751320, 0, -30 ] )
   
    srs = osr.SpatialReference()
    srs.SetUTM( 11, 1 )
    srs.SetWellKnownGeogCS( 'NAD27' )
    dst_ds.SetProjection( srs.ExportToWkt() )
 
    raster = numpy.zeros( (512, 512),dtype=numpy.uint8 )   
    dst_ds.GetRasterBand(1).WriteArray( raster)
 
    # Once we're done, close properly thedataset
    dst_ds = None

GDAL库的细节

关于ReadRaster

ReadAsArray，返回的是Numeric模块定义的Array

例：

dataset =gdal.Open("f:/pyproj/gtif/aspect.tif")

band =dataset.GetRasterBand(1)

band.ReadAsArray(100,100,5,5,10,10)

参数说明：

头两个100是取值窗口左上角在实际数据中所处象元的x，y位置

两个5是取值窗口覆盖面积的大小

两个10是取值窗口取出数组进行缩放后数组的大小（可以不指定，不指定同前两个参数）

注意：不要访问越界，越界会报错

cocolalala

关注

0
点赞
踩
2

收藏

觉得还不错? 一键收藏
0
评论
GDAL库初学

打开文件在打开一个GDAL支持的栅格资料之前，必需要注册驱动。每个驱动对应各自支持的格式。通常这个会被GDALAllRegister（）函数完成。 C++Python import gdal from gdalconst import * dataset = gdal.Open( filename, GA_ReadOnly )
复制链接

扫一扫