GDAL API入门

最新推荐文章于 2023-04-01 14:00:04 发布

weixin_30530339

最新推荐文章于 2023-04-01 14:00:04 发布

阅读量114

点赞数

文章标签： python

原文链接：http://www.cnblogs.com/liuyunfeifei/archive/2012/11/29/2794947.html

版权

GDAL API入门

翻译：柴树杉(chaishushan@gmail.com)
原文：http://www.gdal.org/gdal_tutorial.html

打开文件

在打开GDAL所支持的光栅数据之前需要注册驱动。这里的驱动是针对GDAL支持的所有数据格式。通常可以通过调用 GDALAllRegister() 函数来注册所有已知的驱动，同时也包含那些用 GDALDriverManager::AutoLoadDrivers() 从.so文件中自动装载驱动。如果程序需要对某些驱动做限制，可以参考 gdalallregister.cpp 代码。

当驱动被注册之后，我们就可以用 GDALOpen() 函数来打开一个数据集。打开的方式可以是 GA_ReadOnly 或者 GA_Update。

In C++:

 1 #include "gdal_priv.h"
 2 
 3 int main()
 4 {
 5     GDALDataset  *poDataset;
 6 
 7     GDALAllRegister();
 8 
 9     poDataset = (GDALDataset *) GDALOpen( pszFilename, GA_ReadOnly );
10     if( poDataset == NULL )
11     {
12         ...;
13     }

In C:

 1 #include "gdal.h"
 2 
 3 int main()
 4 {
 5     GDALDatasetH  hDataset;
 6 
 7     GDALAllRegister();
 8 
 9     hDataset = GDALOpen( pszFilename, GA_ReadOnly );
10     if( hDataset == NULL )
11     {
12         ...;
13     }

In Python:

1     import gdal
2     from gdalconst import *
3 
4     dataset = gdal.Open( filename, GA_ReadOnly )
5     if dataset is None:
6         ...

如果 GDALOpen() 函数返回NULL则表示打开失败，同时 CPLError() 函数产生相应的错误信息。如果您需要对错误进行处理可以参考 CPLError() 相关文档。通常情况下，所有的 GDAL函数都通过CPLError()报告错误。另外需要注意的是pszFilename并不一定对应一个实际的文件名（当然也可以就是一个文件名）。它的具体解释由相应的驱动程序负责。它可能是一个URL，或者是文件名以后后面带有许多用于控制打开方式的参数。通常建议，不要在打开文件的选择对话框中对文件的类型做太多的限制。

获取Dataset信息

如果GDAL数据模型一节所描述的，一个GDALDataset包含了光栅数据的一系列的波段信息。同时它还包含元数据、一个坐标系统、投影类型、光栅的大小以及其他许多信息。

1     adfGeoTransform[0] /* 左上角 x */
2     adfGeoTransform[1] /* 东西方向一个像素对应的距离 */
3     adfGeoTransform[2] /* 旋转, 0表示上面为北方 */
4     adfGeoTransform[3] /* 左上角 y */
5     adfGeoTransform[4] /* 旋转, 0表示上面为北方 */
6     adfGeoTransform[5] /* 南北方向一个像素对应的距离 */

如果需要输出dataset的基本信息，可以这样：

In C++:

 1  　 double        adfGeoTransform[6];
 2 
 3     printf( "Driver: %s/%s\n",
 4             poDataset->GetDriver()->GetDescription(), 
 5             poDataset->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME ) );
 6 
 7     printf( "Size is %dx%dx%d\n", 
 8             poDataset->GetRasterXSize(), poDataset->GetRasterYSize(),
 9             poDataset->GetRasterCount() );
10 
11     if( poDataset->GetProjectionRef()  != NULL )
12         printf( "Projection is `%s'\n", poDataset->GetProjectionRef() );
13 
14     if( poDataset->GetGeoTransform( adfGeoTransform ) == CE_None )
15     {
16         printf( "Origin = (%.6f,%.6f)\n",
17                 adfGeoTransform[0], adfGeoTransform[3] );
18 
19         printf( "Pixel Size = (%.6f,%.6f)\n",
20                 adfGeoTransform[1], adfGeoTransform[5] );
21     }

In C:

 1 GDALDriverH   hDriver;
 2     double adfGeoTransform[6];
 3 
 4     hDriver = GDALGetDatasetDriver( hDataset );
 5     printf( "Driver: %s/%s\n",
 6             GDALGetDriverShortName( hDriver ),
 7             GDALGetDriverLongName( hDriver ) );
 8 
 9     printf( "Size is %dx%dx%d\n",
10             GDALGetRasterXSize( hDataset ), 
11             GDALGetRasterYSize( hDataset ),
12             GDALGetRasterCount( hDataset ) );
13 
14     if( GDALGetProjectionRef( hDataset ) != NULL )
15         printf( "Projection is `%s'\n", GDALGetProjectionRef( hDataset ) );
16 
17     if( GDALGetGeoTransform( hDataset, adfGeoTransform ) == CE_None )
18     {
19         printf( "Origin = (%.6f,%.6f)\n",
20                 adfGeoTransform[0], adfGeoTransform[3] );
21 
22         printf( "Pixel Size = (%.6f,%.6f)\n",
23                 adfGeoTransform[1], adfGeoTransform[5] );
24     }

In Python:

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

获取一个光栅波段

现在，我们可以通过GDAL获取光栅的一个波段。同样每个波段含有元数据、块大小、颜色表以前其他一些信息。下面的代码从dataset获取一个GDALRasterBand对象，并且显示它的一些信息。

In C++:

 1         GDALRasterBand  *poBand;
 2         int             nBlockXSize, nBlockYSize;
 3         int             bGotMin, bGotMax;
 4         double          adfMinMax[2];
 5         
 6         poBand = poDataset->GetRasterBand( 1 );
 7         poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
 8         printf( "Block=%dx%d Type=%s, ColorInterp=%s\n",
 9                 nBlockXSize, nBlockYSize,
10                 GDALGetDataTypeName(poBand->GetRasterDataType()),
11                 GDALGetColorInterpretationName(
12                     poBand->GetColorInterpretation()) );
13 
14         adfMinMax[0] = poBand->GetMinimum( &bGotMin );
15         adfMinMax[1] = poBand->GetMaximum( &bGotMax );
16         if( ! (bGotMin && bGotMax) )
17             GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);
18 
19         printf( "Min=%.3fd, Max=%.3f\n", adfMinMax[0], adfMinMax[1] );
20         
21         if( poBand->GetOverviewCount() > 0 )
22             printf( "Band has %d overviews.\n", poBand->GetOverviewCount() );
23 
24         if( poBand->GetColorTable() != NULL )
25             printf( "Band has a color table with %d entries.\n", 
26                      poBand->GetColorTable()->GetColorEntryCount() );

In C:

 1         GDALRasterBandH hBand;
 2         int             nBlockXSize, nBlockYSize;
 3         int             bGotMin, bGotMax;
 4         double          adfMinMax[2];
 5         
 6         hBand = GDALGetRasterBand( hDataset, 1 );
 7         GDALGetBlockSize( hBand, &nBlockXSize, &nBlockYSize );
 8         printf( "Block=%dx%d Type=%s, ColorInterp=%s\n",
 9                 nBlockXSize, nBlockYSize,
10                 GDALGetDataTypeName(GDALGetRasterDataType(hBand)),
11                 GDALGetColorInterpretationName(
12                     GDALGetRasterColorInterpretation(hBand)) );
13 
14         adfMinMax[0] = GDALGetRasterMinimum( hBand, &bGotMin );
15         adfMinMax[1] = GDALGetRasterMaximum( hBand, &bGotMax );
16         if( ! (bGotMin && bGotMax) )
17             GDALComputeRasterMinMax( hBand, TRUE, adfMinMax );
18 
19         printf( "Min=%.3fd, Max=%.3f\n", adfMinMax[0], adfMinMax[1] );
20         
21         if( GDALGetOverviewCount(hBand) > 0 )
22             printf( "Band has %d overviews.\n", GDALGetOverviewCount(hBand));
23 
24         if( GDALGetRasterColorTable( hBand ) != NULL )
25             printf( "Band has a color table with %d entries.\n", 
26                      GDALGetColorEntryCount(
27                          GDALGetRasterColorTable( hBand ) ) );

In Python:

 1         band = dataset.GetRasterBand(1)
 2 
 3         print 'Band Type=',gdal.GetDataTypeName(band.DataType)
 4 
 5         min = band.GetMinimum()
 6         max = band.GetMaximum()
 7         if min is not None and max is not None:
 8             (min,max) = ComputeRasterMinMax(1)
 9         print 'Min=%.3f, Max=%.3f' % (min,max)
10 
11         if band.GetOverviewCount() > 0:
12             print 'Band has ', band.GetOverviewCount(), ' overviews.'
13 
14         if not band.GetRasterColorTable() is None:
15             print 'Band has a color table with ', \
16             band.GetRasterColorTable().GetCount(), ' entries.'

读光栅数据

GDAL有几种读光栅数据的方法，但是GDALRasterBand::RasterIO()是最常用的一种。该函数可以自动转换数据类型、采样以及裁剪。下面的代码读光栅的第1行数据，同时转换为float保存到缓冲。

In C++:

1         float *pafScanline;
2         int   nXSize = poBand->GetXSize();
3 
4         pafScanline = (float *) CPLMalloc(sizeof(float)*nXSize);
5         poBand->RasterIO( GF_Read, 0, 0, nXSize, 1, 
6                           pafScanline, nXSize, 1, GDT_Float32, 
7                           0, 0 );

In C:

1         float *pafScanline;
2         int   nXSize = GDALGetRasterBandXSize( hBand );
3 
4         pafScanline = (float *) CPLMalloc(sizeof(float)*nXSize);
5         GDALRasterIO( hBand, GF_Read, 0, 0, nXSize, 1, 
6                       pafScanline, nXSize, 1, GDT_Float32, 
7                       0, 0 );

In Python:

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

返回的是一个string，包含了xsize*4大小的二进制数据，是float类型指针。可以使用python的struct模块转换为python数据类型：

1         import struct
2 
3         tuple_of_floats = struct.unpack('f' * b2.XSize, scanline)

RasterIO函数的完整说明如下：

1         CPLErr GDALRasterBand::RasterIO( GDALRWFlag eRWFlag,
2                                  int nXOff, int nYOff, int nXSize, int nYSize,
3                                  void * pData, int nBufXSize, int nBufYSize,
4                                  GDALDataType eBufType,
5                                  int nPixelSpace,
6                                  int nLineSpace )

RasterIO()可以通过指定eRWFlag参数来确定是读/写数据（GF_Read或GF_Write）。参数nXOff/nYOff/nXSize/nYSize描述了要读的影象范围（或者是写）。同时它也可以自动处理边界等特殊情况。

参数pData指定读/写对应的缓冲。缓冲的类型必须是eBufType中定义的，例如GDT_Float32、GDT_Byte等。RasterIO ()会自动转换缓冲和波段的类型，使它们一致。当数据向下转换时，或者是数据超出转换后的数据类型可以表示的范围时，将会用最接近的数据来代替。例如一个 16位的整数被转换为GDT_Byte时，所有大于255的值都会用255代替（数据并不会被缩放）。

参数nBufXSize和nBufYSize描述了缓冲的大小。当时读写是是全部数据时，该值和影象的大小相同。当需要对影象抽样的时候，缓冲也可以比真实的影象小。因此，利用RasterIO()实现预览功能是很方便的。

参数nPixelSpace和nLineSpace通常被设置为0。当然，也可以使用他们来控制内存中的数据。关闭Dataset

需要强调的一点是：GDALRasterBand对象属于相应的dataset，用户不能私自delete 任何GDALRasterBand对象。GDALDataset可以用GDALClose()关闭数据，或者是直接 delete GDALDataset对象。关闭GDALDataset的时候会进行相关的清除操作和刷新一些写操作。

创建文件的技巧

如果相应格式的驱动支持写操作的话，则可以创建文件。GDAL有两函数可以创建文件： CreateCopy()和Create()。 CreateCopy()函数直接从参数给定的数据集复制数据。 Create()函数则需要用户明确地写入各种数据（元数据、光栅数据等）。所有支持创建的格式驱动都支持CreateCopy()函数，但是并不一定支持Create()函数。

为了确定数据格式是否支持Create或CreateCopy，可以检查驱动对象中的DCAP_CREATE 和DCAP_CREATECOPY元数据。在使用GetDriverByName()函数之前确保GDALAllRegister() 已经被调用过。

In C++:

 1 #include "cpl_string.h"
 2 ...
 3     const char *pszFormat = "GTiff";
 4     GDALDriver *poDriver;
 5     char **papszMetadata;
 6 
 7     poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
 8 
 9     if( poDriver == NULL )
10         exit( 1 );
11 
12     papszMetadata = poDriver->GetMetadata();
13     if( CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATE, FALSE ) )
14         printf( "Driver %s supports Create() method.\n", pszFormat );
15     if( CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATECOPY, FALSE ) )
16         printf( "Driver %s supports CreateCopy() method.\n", pszFormat );

In C:

 1 #include "cpl_string.h"
 2 ...
 3     const char *pszFormat = "GTiff";
 4     GDALDriver hDriver = GDALGetDriverByName( pszFormat );
 5     char **papszMetadata;
 6 
 7     if( hDriver == NULL )
 8         exit( 1 );
 9 
10     papszMetadata = GDALGetMetadata( hDriver, NULL );
11     if( CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATE, FALSE ) )
12         printf( "Driver %s supports Create() method.\n", pszFormat );
13     if( CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATECOPY, FALSE ) )
14         printf( "Driver %s supports CreateCopy() method.\n", pszFormat );

In Python:

1     format = "GTiff"
2     driver = gdal.GetDriverByName( format )
3     metadata = driver.GetMetadata()
4     if metadata.has_key(gdal.DCAP_CREATE) \
5        and metadata[gdal.DCAP_CREATE] == 'YES':
6         print 'Driver %s supports Create() method.' % format
7     if metadata.has_key(gdal.DCAP_CREATECOPY) \
8        and metadata[gdal.DCAP_CREATECOPY] == 'YES':
9         print 'Driver %s supports CreateCopy() method.' % format

我们可以看出有些格式不支持Create()或CreateCopy()调用。

使用CreateCopy()

GDALDriver::CreateCopy()函数使用比较简单，并且原先数据中的所有信息都被正确的设置。函数还可以指定某些可的选择参数，也通过一个回调函数来获得数据复制的进展情况。下面的程序用默认的方式copy一个pszSrcFilename文件，保存为 pszDstFilename 文件。

In C++:

1     GDALDataset *poSrcDS = (GDALDataset *) GDALOpen( pszSrcFilename, GA_ReadOnly );
2     GDALDataset *poDstDS = poDriver->CreateCopy( pszDstFilename, poSrcDS, FALSE, NULL, NULL, NULL );
3     if( poDstDS != NULL )
4         delete poDstDS;

In C:

1     GDALDatasetH hSrcDS = GDALOpen( pszSrcFilename, GA_ReadOnly );
2     GDALDatasetH hDstDS = GDALCreateCopy( hDriver, pszDstFilename, hSrcDS, FALSE, NULL, NULL, NULL );
3     if( hDstDS != NULL )
4         GDALClose( hDstDS );

In Python:

1     src_ds = gdal.Open( src_filename )
2     dst_ds = driver.CreateCopy( dst_filename, src_ds, 0 )

CreateCopy()返回一个可写入的dataset，并且返回的dataset最终需要用户自己关闭（和delete）以保证数据被真正地写入磁盘（dataset本身可能有缓冲）。参数FALSE表示当转换到输出格式时遇到不匹配或者丢失数据时，CreateCopy()宽大处理。这主要是因为输出格式可能不支持输入的数据类型，或者是不支持写操作。

一个更复杂的处理方式是指定某些选项，并且用预定义的回调函数获得进度。

In C++:

 1     #include "cpl_string.h"
 2     ...
 3     char **papszOptions = NULL;
 4     
 5     papszOptions = CSLSetNameValue( papszOptions, "TILED", "YES" );
 6     papszOptions = CSLSetNameValue( papszOptions, "COMPRESS", "PACKBITS" );
 7     poDstDS = poDriver->CreateCopy( pszDstFilename, poSrcDS, FALSE, 
 8                                     papszOptions, GDALTermProgress, NULL );
 9     if( poDstDS != NULL )
10         delete poDstDS;

In C:

 1     #include "cpl_string.h"
 2     ...
 3     char **papszOptions = NULL;
 4     
 5     papszOptions = CSLSetNameValue( papszOptions, "TILED", "YES" );
 6     papszOptions = CSLSetNameValue( papszOptions, "COMPRESS", "PACKBITS" );
 7     hDstDS = GDALCreateCopy( hDriver, pszDstFilename, hSrcDS, FALSE, 
 8                              papszOptions, GDALTermProgres, NULL );
 9     if( hDstDS != NULL )
10         GDALClose( hDstDS );

In Python:

1     src_ds = gdal.Open( src_filename )
2     dst_ds = driver.CreateCopy( dst_filename, src_ds, 0, 
3                                 [ 'TILED=YES', 'COMPRESS=PACKBITS' ] )

使用Create()

如果你不是简单地复制一个文件的话，就可能需要使用GDALDriver::Create() 来创建文件。Create()的参数列表和CreateCopy()相似，但是需要明确指定影象的大小、波段数以及波段数据类型。

In C++:

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

In C:

1     GDALDatasetH hDstDS;        
2     char **papszOptions = NULL;
3 
4     hDstDS = GDALCreate( hDriver, pszDstFilename, 512, 512, 1, GDT_Byte,  papszOptions );

In Python:

1     dst_ds = driver.Create( dst_filename, 512, 512, 1, gdal.GDT_Byte )

当dataset被正确地创建之后，特定的元数据和光栅数据都要被写到文件中。这些操作一般需要依赖用户的具体选择，下边的代码是一个简单示例。

In C++:

 1     double adfGeoTransform[6] = { 444720, 30, 0, 3751320, 0, -30 };
 2     OGRSpatialReference oSRS;
 3     char *pszSRS_WKT = NULL;
 4     GDALRasterBand *poBand;
 5     GByte abyRaster[512*512];
 6 
 7     poDstDS->SetGeoTransform( adfGeoTransform );
 8     
 9     oSRS.SetUTM( 11, TRUE );
10     oSRS.SetWellKnownGeogCS( "NAD27" );
11     oSRS.exportToWkt( &pszSRS_WKT );
12     poDstDS->SetProjection( pszSRS_WKT );
13     CPLFree( pszSRS_WKT );
14 
15     poBand = poDstDS->GetRasterBand(1);
16     poBand->RasterIO( GF_Write, 0, 0, 512, 512, 
17                       abyRaster, 512, 512, GDT_Byte, 0, 0 );    
18 
19     delete poDstDS;

In C:

 1     double adfGeoTransform[6] = { 444720, 30, 0, 3751320, 0, -30 };
 2     OGRSpatialReferenceH hSRS;
 3     char *pszSRS_WKT = NULL;
 4     GDALRasterBandH hBand;
 5     GByte abyRaster[512*512];
 6 
 7     GDALSetGeoTransform( hDstDS, adfGeoTransform );
 8 
 9     hSRS = OSRNewSpatialReference( NULL );
10     OSRSetUTM( hSRS, 11, TRUE );
11     OSRSetWellKnownGeogCS( hSRS, "NAD27" );                     
12     OSRExportToWkt( hSRS, &pszSRS_WKT );
13     OSRDestroySpatialReference( hSRS );
14 
15     GDALSetProjection( hDstDS, pszSRS_WKT );
16     CPLFree( pszSRS_WKT );
17 
18     hBand = GDALGetRasterBand( hDstDS, 1 );
19     GDALRasterIO( hBand, GF_Write, 0, 0, 512, 512, 
20                   abyRaster, 512, 512, GDT_Byte, 0, 0 );    
21 
22     GDALClose( hDstDS );

In Python:

 1     import Numeric, osr
 2 
 3     dst_ds.SetGeoTransform( [ 444720, 30, 0, 3751320, 0, -30 ] )
 4     
 5     srs = osr.SpatialReference()
 6     srs.SetUTM( 11, 1 )
 7     srs.SetWellKnownGeogCS( 'NAD27' )
 8     dst_ds.SetProjection( srs.ExportToWkt() )
 9 
10     raster = Numeric.zeros( (512, 512) )    
11     dst_ds.GetRasterBand(1).WriteArray( raster )

转载于:https://www.cnblogs.com/liuyunfeifei/archive/2012/11/29/2794947.html

weixin_30530339

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
GDAL API入门

GDAL API入门翻译：柴树杉(chaishushan@gmail.com)原文：http://www.gdal.org/gdal_tutorial.html打开文件在打开GDAL所支持的光栅数据之前需要注册驱动。这里的驱动是针对GDAL支持的所有数据格式。通常可以通过调用GDALAllRegister()函数来注册所有已知的驱动，同时也包含那些用GDALDr...
复制链接

扫一扫