光栅API教程 python

原地址：https://gdal.org/tutorials/raster_api_tut.html

打开文件：

在打开GDAL支持的栅格数据存储之前，需要注册驱动程序。每种支持的格式都有一个驱动程序。通常这是通过GDALAllRegister()函数完成的，该函数尝试注册所有已知的驱动程序，包括那些使用GDALDriverManager::AutoLoadDrivers()从.so文件自动加载的驱动程序。如果某些应用程序需要限制驱动程序的设置，那么检查来自gdalallregister.cpp的代码可能会有帮助。当导入gdal模块时，Python自动调用GDALAllRegister()。一旦注册了驱动程序，应用程序应该调用独立的GDALOpen()函数来打开数据集，传递数据集的名称和所需的访问权限(GA_ReadOnly或GA_Update)。

from osgeo import gdal
dataset = gdal.Open(filename, gdal.GA_ReadOnly)

注意，如果GDALOpen()返回NULL，则意味着打开失败，并且已经通过CPLError()发出了错误消息。如果您想控制如何向用户报告错误，请查看CPLError()文档。一般来说，所有的GDAL都使用CPLError()来进行错误报告。另外，请注意，pszFilename实际上不需要是物理文件的名称(尽管它通常是)。它的解释是依赖于驱动程序的，它可能是一个URL，一个在最后添加了控制打开或几乎任何其他参数的文件名。请尽量不要将GDAL文件选择对话框限制为只选择物理文件。

获取数据集信息

如栅格数据模型中所述，GDALDataset包含栅格波段列表，所有栅格波段都属于同一区域，具有相同的分辨率。它还拥有元数据、坐标系统、地理坐标变换、光栅大小和其他各种信息。

没有任何旋转或剪切的“北上”图像的特殊但常见的情况下，地理坐标变换Geotransform教程采用以下形式:

adfGeoTransform[0] /* top left x */
adfGeoTransform[1] /* w-e pixel resolution */
adfGeoTransform[2] /* 0 */
adfGeoTransform[3] /* top left y */
adfGeoTransform[4] /* 0 */
adfGeoTransform[5] /* n-s pixel resolution (negative value) */

在一般情况下，这是一个仿射变换。
如果我们想打印一些关于数据集的一般信息，我们可以做以下事情:

print("Driver: {}/{}".format(dataset.GetDriver().ShortName,
                            dataset.GetDriver().LongName))
print("Size is {} x {} x {}".format(dataset.RasterXSize,
                                    dataset.RasterYSize,
                                    dataset.RasterCount))
print("Projection is {}".format(dataset.GetProjection()))
geotransform = dataset.GetGeoTransform()
if geotransform:
    print("Origin = ({}, {})".format(geotransform[0], geotransform[3]))
    print("Pixel Size = ({}, {})".format(geotransform[1], geotransform[5]))

获取光栅波段

在这个时候，通过GDAL访问光栅数据是在一个波段完成的。此外，还有元数据、块大小、颜色表和不同波段的其他信息。以下代码从数据集(通过GDALRasterBand::GetRasterCount()编号为1)获取GDALRasterBand对象，并显示有关该对象的一些信息。

band = dataset.GetRasterBand(1)
print("Band Type={}".format(gdal.GetDataTypeName(band.DataType)))

min = band.GetMinimum()
max = band.GetMaximum()
if not min or not max:
    (min,max) = band.ComputeRasterMinMax(True)
print("Min={:.3f}, Max={:.3f}".format(min,max))

if band.GetOverviewCount() > 0:
    print("Band has {} overviews".format(band.GetOverviewCount()))

if band.GetRasterColorTable():
    print("Band has a color table with {} entries".format(band.GetRasterColorTable().GetCount()))

阅读栅格数据

有几种读取光栅数据的方法，但最常见的是通过GDALRasterBand::RasterIO()方法。此方法将自动处理数据类型转换、上下采样和窗口设置。下面的代码将把数据的第一个扫描行读入大小类似的缓冲区，并作为操作的一部分将其转换为浮点数。

scanline = band.ReadRaster(xoff=0, yoff=0,
                        xsize=band.XSize, ysize=1,
                        buf_xsize=band.XSize, buf_ysize=1,
                        buf_type=gdal.GDT_Float32)

注意，返回的扫描线类型为string，并包含xsize*4字节的原始二进制浮点数据。这可以转换为Python值使用struct模块从标准库:

import struct
tuple_of_floats = struct.unpack('f' * b2.XSize, scanline)

RasterIO调用接受以下参数:

CPLErr GDALRasterBand::RasterIO( GDALRWFlag eRWFlag,
                                int nXOff, int nYOff, int nXSize, int nYSize,
                                void * pData, int nBufXSize, int nBufYSize,
                                GDALDataType eBufType,
                                int nPixelSpace,
                                int nLineSpace )

注意，相同的RasterIO()调用用于基于eRWFlag的设置进行读取或写入(要么是GF_Read，要么是GF_Write)。nXOff、nYOff、nXSize、nYSize参数描述了磁盘上要读(或写)的光栅数据窗口。它不必落在平铺边界，但访问可能更有效，如果它。

pData是读入或写入数据的内存缓冲区。它的真实类型必须是作为eBufType传递的任何类型，如GDT_Float32或GDT_Byte。RasterIO()调用将负责缓冲区的数据类型和频带的数据类型之间的转换。注意，在将浮点数据转换为整数时，RasterIO()将四舍五入，并且在转换输出的合法范围之外的源值时，将使用最接近的合法值。这意味着，例如，读入GDT_Byte缓冲区的16位数据将把所有大于255的值映射为255，该数据没有缩放!

nBufXSize和nBufYSize值描述了缓冲区的大小。当加载完全分辨率的数据时，这将是相同的窗口大小。但是，要加载一个简化的分辨率概览，可以将其设置为小于磁盘上的窗口。在这种情况下，如果概述合适，RasterIO()将利用概述来更有效地执行IO。

nPixelSpace和nLineSpace通常为零，表示应该使用默认值。但是，它们可以用来控制对内存数据缓冲区的访问，例如允许读取包含其他像素交叉数据的缓冲区。

关闭数据集

请记住，GDALRasterBand对象属于它们的数据集，它们不应该被c++ delete操作符销毁。GDALDataset s可以通过调用GDALClose()来关闭(由于在跨模块边界分配和释放内存时存在已知问题，不建议对Windows用户使用GDALDataset上的delete操作符。参见FAQ中的相关主题)。调用GDALClose将导致正确的清理，并刷新所有挂起的写操作。对于在更新模式下以GTiff等流行格式打开的数据集，忘记调用GDALClose可能会导致以后无法打开它。

创建文件的技术

如果格式驱动程序支持创建，可能会创建GDAL支持的格式的新文件。创建文件有两种常用技术，使用CreateCopy()和Create()。CreateCopy方法涉及到调用格式驱动程序上的CreateCopy()方法，并传递一个应该复制的源数据集。Create方法包括调用驱动程序上的Create()方法，然后显式地编写所有元数据，以及使用单独的调用编写栅格数据。所有支持创建新文件的驱动程序都支持CreateCopy()方法，但只有少数驱动程序支持Create()方法。

要确定特定格式是否支持创建或CreateCopy，可以检查格式驱动程序对象上的DCAP_CREATE和DCAP_CREATECOPY元数据。确保在调用GDALDriverManager::GetDriverByName()之前已经调用了GDALAllRegister()。在这个示例中，我们获取一个驱动程序，并确定它是否支持Create()和/或CreateCopy()。

fileformat = "GTiff"
driver = gdal.GetDriverByName(fileformat)
metadata = driver.GetMetadata()
if metadata.get(gdal.DCAP_CREATE) == "YES":
    print("Driver {} supports Create() method.".format(fileformat))

if metadata.get(gdal.DCAP_CREATECOPY) == "YES":
    print("Driver {} supports CreateCopy() method.".format(fileformat))

注意，许多驱动程序是只读的，并且不支持Create()或CreateCopy()。

使用CreateCopy ()

GDALDriver::CreateCopy()方法使用起来相当简单，因为大多数信息是从源数据集收集的。但是，它包含一些选项，用于传递特定于格式的创建选项，以及在进行长数据集复制时向用户报告进度。简单的复制一个名为pszSrcFilename的文件，到一个名为pszDstFilename的新文件，使用默认选项的格式，其驱动程序之前可能是这样的:

src_ds = gdal.Open(src_filename)
dst_ds = driver.CreateCopy(dst_filename, src_ds, strict=0)
# Once we're done, close properly the dataset
dst_ds = None
src_ds = None

请注意，CreateCopy()方法返回一个可写数据集，必须正确关闭该数据集，才能完成对数据集的写入并将其刷新到磁盘。在Python中，当dst_ds超出范围时，就会自动发生这种情况。在CreateCopy()调用中的目标文件名之后，bStrict选项使用的FALSE(或0)值表明，即使无法创建目标数据集来精确匹配输入数据集，也应该继续执行CreateCopy()调用，而不会出现致命错误。这可能是因为输出格式不支持输入数据集的像素数据类型，或者是因为目标不支持写入地理坐标。

更复杂的情况可能包括传递创建选项，并使用预定义的进度监视器，如下所示:

src_ds = gdal.Open(src_filename)
dst_ds = driver.CreateCopy(dst_filename, src_ds, strict=0,
                        options=["TILED=YES", "COMPRESS=PACKBITS"])
# Once we're done, close properly the dataset
dst_ds = None
src_ds = None

使用Create ()

对于不只是将现有文件导出到新文件的情况，通常需要使用GDALDriver::Create()方法(尽管可以通过使用虚拟文件或内存文件来实现一些有趣的选项)。Create()方法接受一个类似于CreateCopy()的选项列表，但是必须显式地提供图像大小、频带数量和频带类型。

dst_ds = driver.Create(dst_filename, xsize=512, ysize=512,
                    bands=1, eType=gdal.GDT_Byte)

成功创建数据集后，必须将所有适当的元数据和栅格数据写入文件。这是什么将根据使用情况而有所不同，但是这里介绍一个简单的案例，包括投影、geotransform和光栅数据.

from osgeo import osr
import numpy
dst_ds.SetGeoTransform([444720, 30, 0, 3751320, 0, -30])
srs = osr.SpatialReference()
srs.SetUTM(11, 1)
srs.SetWellKnownGeogCS("NAD27")
dst_ds.SetProjection(srs.ExportToWkt())
raster = numpy.zeros((512, 512), dtype=numpy.uint8)
dst_ds.GetRasterBand(1).WriteArray(raster)
# Once we're done, close properly the dataset
dst_ds = None