注意，对于 conda 环境，会出现找不到 proj 的警告：
处理方法参考 https://stackoverflow.com/questions/72021249/running-a-python-docker-with-miniconda3-and-geopandas
更多说明参考： https://gis.stackexchange.com/questions/364421/how-to-make-proj-work-via-anaconda-in-google-colab/
通过如下方法声明：
import os
# os.environ['PROJ_LIB'] = '/git/usr/miniforge/envs/python312/share/proj'
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 )
/opt/conda/lib/python3.12/site-packages/osgeo/gdal.py:311: FutureWarning: Neither gdal.UseExceptions() nor gdal.DontUseExceptions() has been explicitly called. In GDAL 4.0, exceptions will be enabled by default. warnings.warn(
上面语句创建了一个GeoTIFF格式的栅格影像。 宽度是 $512\times 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
from osgeo import gdal
import numpy
dataset = gdal.Open("/data/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。
from osgeo import osr
from osgeo import gdal
dataset = gdal.Open("/data/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了。
我这里用的空间参考是美国常用的。
如果你使用的空间参考比较麻烦，可以只定义六参数变换。
from osgeo import gdal
gdal.SetConfigOption('HFA_USE_RRD', 'YES')
强制建立pyramids:
tods.BuildOverviews(overviewlist=[2,4, 8,16,32,64,128])
0
使用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 = "/data/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("/data/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("/data/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) 。这样就可以正常显示了。
具体使用哪种存储方式取决于使用的目的。