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 )