>>> from env_helper import info; info()
页面更新时间： 2024-07-23 22:39:29
运行环境：
Linux发行版本: Debian GNU/Linux 12 (bookworm)
操作系统内核: Linux-6.1.0-23-amd64-x86_64-with-glibc2.36
Python版本: 3.11.2
2.6. 使用GDAL创建影像¶
2.6.1. 使用Create函数创建影像¶
如果不是使用一个已有的影像文件来创建新的影像，就需要用到 Create
方法了。
在数据处理过程中，这种是主要的方法，它可以把建立在内存中的虚拟数据集输出到实际文件。
也就是栅格数据持久化的概念，将内存中的数据模型（主要是二维数组）转换为存储模型，
对于地理信息，除了数据本身，还有投影、元数据信息等。
help(driver.Create)
可以看出来，这个函数和 CreateCopy 很像，不过它多了几个参数，
xsize,ysize是图像大小，bands是图像的波段(通道)数，
datatype 就是图像的象元数据类型了。
>>> from osgeo 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。作为示例，这里没使用源数据。 它创建了一个空的数据集。实际的数据还需要另外的代码。
>>> import numpy
>>> from osgeo import osr
>>> dst_ds.SetGeoTransform( [ 444720, 30, 0, 3751320, 0, -30 ] )
0
>>> srs = osr.SpatialReference()
>>> srs.SetUTM( 11, 1 )
0
>>> srs.SetWellKnownGeogCS( 'NAD27' )
0
>>> dst_ds.SetProjection( srs.ExportToWkt() )
0
>>> raster = numpy.zeros( (512, 512) )
>>> dst_ds.GetRasterBand(1).WriteArray( raster )
0
这里的例子设置了空间范围和坐标系(关于这些下面的文章再重点介绍)，
还有一个 \(512 \times 512\) 的全部都是 0 的数组数据。
当然这样做，你在打开数据的时候只能看到一片黑色。
如果想要使你的数据有点实质性参考价值的话，
需要把一个丰富多彩的数组塞到图片当中。
当然最好的办法就是从已有的数据中读取了。
除了从原数据中老老实实读出数据后，还可以进行矩阵魔术，
或进行一些计算，如前面提到的NDVI的计算。
使其满足我们的需要。然后再写到磁盘上。
>>> 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 )
0
2.6.2. 创建多波段图像¶
创建多波段影像的方法¶
大多数的遥感影像都是多波段的。每一个波段记录了不同的波谱信息，
对于影像处理结果，也可以使用多波段的方式来存储。
下面介绍一个建立3波段GTiff的小例子。多波段影像的创建方式与之类似。
>>> 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.tobytes(),width,height,band_list=[1,2,3])
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)
分波段处理¶
当然从波段里读取数据再拼接成完整的图像更是可以的。
>>> 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。
2.6.3. GDAL写操作的其他问题¶
使用 GDAL 创建数据时要注意，影像的空间参数信息是单独处理的。 比如我们导入一个NAD27的空间参考，我们可以这样写：
>>> from osgeo import osr
>>> from osgeo 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() )
0
于是TIFF 就变成了 GeoTIFF 。 空间参考系统是NAD27，起点
(444720，3751320) ，像素大小30的TIFF了。
我这里用的空间参考是美国常用的。
如果你使用的空间参考比较麻烦，可以只定义六参数变换。
2.6.4. 建立影像金字塔¶
设定Imagine风格的pyramids:
>>> from osgeo import gdal
>>> gdal.SetConfigOption('HFA_USE_RRD', 'YES')
强制建立pyramids:
>>> tods.BuildOverviews(overviewlist=[2,4, 8,16,32,64,128])
0
2.6.5. 使用CreateCopy函数创建影像¶
CreateCopy()
的使用比较简单，就是把一个格式的图像直接转化为另一个格式的图像，
相当于一种格式转换。 看看 CreateCopy() 的原型：
CreateCopy(self, filename, source_ds, strict=1, options=[], callback=None, callb
第一个参数是源文件的名称，第二个是目标文件名称。第三及以后的参数都可选。 如果不输入也可以，程序按照默认的方式运行。第三个参数取值是0或者1(True或者False)。 取值为非的时候说明即使不能精确的由原数据转化为目标数据， 程序也照样执行CreateCopy方法，不会产生致命错误。 这种错误有可能是输出格式不支持输入数据格式象元的数据类型， 或者是目标数据不支持写入空间参考等等。下面是一个例子:
>>> from osgeo 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 )
这个函数还有第四个参数， 第四个参数是指在格式转换中所要用到的一些特殊的参数。
>>> 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的过程中有可能出现问题。比如：
>>> 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 ， 也就是 RRRR……,GGGG……,BBBB…… 模式显示。
但是在一些软件中只认值1，也就是 RGBRGBRGBRGB……
，所以上面的代码需要修改成：
>>> 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) 。这样就可以正常显示了。
具体使用哪种存储方式取决于使用的目的。