使用GDAL创建影像

使用Create函数创建影像

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

help(driver.Create)

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

In [20]:
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 [21]:
import numpy, osr
dst_ds.SetGeoTransform( [ 444720, 30, 0, 3751320, 0, -30 ] )  
Out[21]:
0
In [22]:
srs = osr.SpatialReference()
srs.SetUTM( 11, 1 )
Out[22]:
0
In [23]:
srs.SetWellKnownGeogCS( 'NAD27' )
Out[23]:
0
In [24]:
dst_ds.SetProjection( srs.ExportToWkt() )
Out[24]:
0
In [25]:
raster = numpy.zeros( (512, 512) )
dst_ds.GetRasterBand(1).WriteArray( raster )   
Out[25]:
0

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

In [26]:
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[26]:
0

创建多波段图像

创建多波段影像的方法

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

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

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

In [27]:
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[27]:
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 [28]:
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 [29]:
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[29]:
0

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

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

建立影像金字塔

设定Imagine风格的pyramids:

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

强制建立pyramids:

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

使用CreateCopy函数创建影像

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

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

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

In [32]:
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 [33]:
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 [34]:
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 [35]:
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"]);

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

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