1. PythonGIS
  2. GDAL
  3. 使用GDAL创建影像

使用GDAL创建影像

使用Create函数创建影像

如果不是使用一个已有的影像文件来创建新的影像，就需要用到 Create 方法了。 在数据处理过程中，这种是主要的方法，它可以把建立在内存中的虚拟数据集输出到实际文件。 也就是栅格数据持久化的概念，将内存中的数据模型（主要是二维数组）转换为存储模型， 对于地理信息，除了数据本身，还有投影、元数据信息等。

help(driver.Create)

可以看出来，这个函数和 CreateCopy 很像，不过它多了几个参数， xsize,ysize是图像大小，bands是图像的波段(通道)数， datatype 就是图像的象元数据类型了。

In [1]:
import gdal
driver = gdal.GetDriverByName( 'GTiff' )
dst_filename = '/tmp/x_tmp.tif'
dst_ds = driver.Create( dst_filename, 512, 512, 1, gdal.GDT_Byte )

上面语句创建了一个GeoTIFF格式的栅格影像。宽度是:math:` 512times 512`,单波段，数据类型是Byte。作为示例，这里没使用源数据。 它创建了一个空的数据集。实际的数据还需要另外的代码。

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

这里的例子设置了空间范围和坐标系(关于这些下面的文章再重点介绍)， 还有一个:math:` 512times 512`的全部都是0的数组数据。 当然这样做，你在打开数据的时候只能看到一片黑色。 如果想要使你的数据有点实质性参考价值的话， 需要把一个丰富多彩的数组塞到图片当中。 当然最好的办法就是从已有的数据中读取了。 除了从原数据中老老实实读出数据后，还可以进行矩阵魔术， 或进行一些计算，如前面提到的NDVI的计算。 使其满足我们的需要。然后再写到磁盘上。

In [7]:
srs.SetWellKnownGeogCS( 'NAD27' )
srs = osr.SpatialReference()
srs.SetUTM( 11, 1 )
srs.SetWellKnownGeogCS( 'NAD27' )
dst_ds.SetProjection( srs.ExportToWkt() )
raster = numpy.zeros( (512, 512) )
dst_ds.GetRasterBand(1).WriteArray( raster )
Out[7]:
0

创建多波段图像

创建多波段影像的方法

大多数的遥感影像都是多波段的。每一个波段记录了不同的波谱信息，

对于影像处理结果，也可以使用多波段的方式来存储。

下面介绍一个建立3波段GTiff的小例子。多波段影像的创建方式与之类似。

In [8]:
from osgeo import gdal
import numpy
dataset = gdal.Open("/gdata/geotiff_file.tif")
width = dataset.RasterXSize
height = dataset.RasterYSize
datas = dataset.ReadAsArray(0,0,width,height)

driver = gdal.GetDriverByName("GTiff")
tods = driver.Create("/tmp/x_K52E015007_3.tif",width,height,3,options=["INTERLEAVE=PIXEL"])
tods.WriteRaster(0,0,width,height,datas.tostring(),width,height,band_list=[1,2,3])
Out[8]:
0

这是一个很简单的另存遥感图像的方法（不包括空间信息）。这里尤其注意 Create 函数中的options= ["INTERLEAVE=PIXEL"] 参数。 没有这个参数，波段像素组织会错，保存出的图像只有横向的1/3。而且彩色完全不对 。

注意向 tods 写入数据时，需要转换数据类型 datas.tostring()。 如果读取数据的时候使用下面的数据：

datas = dataset.ReadData(0,0,width,height)

就不用转换了。

另外要注意 band_list 参数，这个参数是由波段数决定的。 这个需要根据源数据进行判断，防止出现异常。

还有一个问题就是数据在内存中的顺序， 应该与 INTERLEAVE 密切相关的。

当然datas可以另外组织，比如这样也可以：

datas = dataset.ReadAsArray(0,0,width,height)
datas = numpy.concatenate(datas)

分波段处理

当然从波段里读取数据再拼接成完整的图像更是可以的。

In [9]:
datas = []
for i in range(3):
    band = dataset.GetRasterBand(i+1)
    data = band.ReadAsArray(0,0,width,height)
    datas.append(numpy.reshape(data,(1,-1)))
datas = numpy.concatenate(datas)

如果需要各个波段输入的话，可以循环到各个波段中，然后用Band对象的 WriteRaster来操作，而非在Dataset中调用 WriteRaster。

GDAL写操作的其他问题

空间参考与投影信息的写入

使用 GDAL 创建数据时要注意，影像的空间参数信息是单独处理的。 比如我们导入一个NAD27的空间参考，我们可以这样写：

In [10]:
import osr
import gdal
dataset = gdal.Open("/gdata/geotiff_file.tif")
width = dataset.RasterXSize
height = dataset.RasterYSize
datas = dataset.ReadAsArray(0,0,width,height)
driver = gdal.GetDriverByName("GTiff")

tods = driver.Create("/tmp/x_K52E015007_3.tif",width,height,3,options=["INTERLEAVE=PIXEL"])
tods.SetGeoTransform( [ 444720, 30, 0, 3751320, 0, -30 ] )
srs = osr.SpatialReference()
srs.SetUTM( 11, 1 )
srs.SetWellKnownGeogCS( 'NAD27' )
tods.SetProjection( srs.ExportToWkt() )
Out[10]:
0

于是TIFF 就变成了 GeoTIFF 。 空间参考系统是NAD27，起点 (444720，3751320) ，像素大小30的TIFF了。 我这里用的空间参考是美国常用的。

如果你使用的空间参考比较麻烦，可以只定义六参数变换。

建立影像金字塔

设定Imagine风格的pyramids:

In [11]:
from osgeo import gdal
gdal.SetConfigOption('HFA_USE_RRD', 'YES')

强制建立pyramids:

In [12]:
tods.BuildOverviews(overviewlist=[2,4, 8,16,32,64,128])
Out[12]:
0

使用CreateCopy函数创建影像

CreateCopy() 的使用比较简单，就是把一个格式的图像直接转化为另一个格式的图像， 相当于一种格式转换。 看看 CreateCopy() 的原型：

CreateCopy(self, filename, source_ds, strict=1, options=[], callback=None, callb

第一个参数是源文件的名称，第二个是目标文件名称。第三及以后的参数都可选。 如果不输入也可以，程序按照默认的方式运行。第三个参数取值是0或者1(True或者False)。 取值为非的时候说明即使不能精确的由原数据转化为目标数据， 程序也照样执行CreateCopy方法，不会产生致命错误。 这种错误有可能是输出格式不支持输入数据格式象元的数据类型， 或者是目标数据不支持写入空间参考等等。下面是一个例子:

In [13]:
import gdal
src_filename = "/gdata/geotiff_file.tif"
dst_filename = "/tmp/x2_geotiff_file.tif"
driver = gdal.GetDriverByName('GTiff')
src_ds = gdal.Open( src_filename )
dst_ds = driver.CreateCopy( dst_filename, src_ds, 0 )

这个函数还有第四个参数， 第四个参数是指在格式转换中所要用到的一些特殊的参数。

In [14]:
dst_filename2 = "/tmp/x3_geotiff_file.tif"
dst_ds = driver.CreateCopy( dst_filename2, src_ds, 0, [ 'TILED=YES', 'COMPRESS=PACKBITS' ] )

dst_filename3 = "/tmp/x4_geotiff_file.tif"
dst_ds3 = driver.CreateCopy( dst_filename3, src_ds, 0, [ 'TILED=YES', 'COMPRESS=PACKBITS' ] )

这个例子就是说明在转换Tiff的过程中用瓦片式存储代替条带式存储。 不同软件支持的存储方式是不一样的。 用的压缩方法是PACKBITS。第四个参数根据各种格式都可能有不同的选项。 所以无法统一列出，在这里可以看到GTiff的全部创建参数， 只有在根据各种不同格式时参考各自的参数和写法。 第五个和第六个参数是登记一个挂钩函数， 使得可以把转换过程反映出来（比如画个进度条）。 这两个函数就不多说了。

像素存储顺序

TIFF格式的文件使用比较多，关于像素存储顺序的问题再多说一点，使用的时候多注意一下。 在转换GTiff的过程中，GTiff是支持使用 TILED 参数的。 如果不加参数，利用缺省方式生成gtiff的过程中有可能出现问题。比如：

In [15]:
dataset = gdal.Open("/gdata/geotiff_file.tif")
width = dataset.RasterXSize
height = dataset.RasterYSize
data = dataset.ReadAsArray(0,0,width,height)
driver = gdal.GetDriverByName("GTiff")
driver.CreateCopy("/tmp/sd.tif",dataset,0);

这样生成的结果，大多数情况下没法直接查看。 条带式存储一般遥感影像使用的比较多， 代码转出的GTiff文件虽然在ESRI的系列软件和其他一些看图程序中可以正常显示。 但是在常用的看图软件中，如Picasa等都不支持，在Windows图像浏览器中也不能正常显示。 究其原因，是GDAL在导出的时候把284号域（PlanarConfiguration域）设成了2 [1]，也就是RRRR……，GGGG……，BBBB……模式显示。 但是在一些软件中只认值1，也就是RGBRGBRGBRGB……，所以上面的代码需要修改成：

In [16]:
dataset = gdal.Open("/gdata/geotiff_file.tif")
width = dataset.RasterXSize
height = dataset.RasterYSize
data = dataset.ReadAsArray(0,0,width,height)
driver = gdal.GetDriverByName("GTiff")
driver.CreateCopy("/tmp/sd2.tif",dataset,0,["INTERLEAVE=PIXEL"]);

默认的 INTERLEAVE 是 BAND(tags[284]=2) ，我们把它改成 PIXEL(tags[284]=1) 。这样就可以正常显示了。

具体使用哪种存储方式取决于使用的目的。

版权归 OSGeo中国中心（http://www.osgeo.cn/ ） 所有
邮箱： bukun@iga.ac.cn