GDAL 文档 » 教程 » 栅格 API 教程

目录

打开文件

获取数据集信息

获取栅格波段

读取栅格数据

关闭数据集

创建文件的技术

使用 CreateCopy（）

使用 Create（）

---

打开文件

---

---

在打开 GDAL 支持的栅格数据存储之前，必须注册驱动程序。每种支持的格式都有一个驱动程序。通常，这是通过GDALAllRegister（）函数完成的，该函数尝试注册所有已知的驱动程序，包括使用GDALDriverManager：：AutoLoadDrivers（）从.so文件自动加载的驱动程序。如果对于某些应用程序，有必要限制驱动程序集，则查看gdalallregister.cpp的代码可能会有所帮助。Python 在导入 gdal 模块时自动调用 GDALAllRegister（）。

注册驱动程序后，应用程序应调用独立的 GDALOpen（） 函数来打开数据集，传递数据集的名称和所需的访问权限（GA_ReadOnly 或 GA_Update）。

在 Python 中：

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

这段代码使用了 Python 中的 gdal 库来读取文件 (filename) 中的栅格数据。
它首先从 osgeo 库中导入 gdal，然后使用 gdal.Open() 函数打开文件，并
将打开方式设置为只读 (gdal.GA_ReadOnly)。如果文件打开失败，dataset 
将为 None。所以在 if not dataset: 语句中，如果 dataset 是 None，
则代码执行相应错误处理的代码。

请注意，如果 GDALOpen（） 返回 NULL，则表示打开失败，并且已经通过 CPLError（） 发出错误消息。如果要控制向用户报告错误的方式，请查看 CPLError（） 文档。一般来说，所有的GDAL都使用CPLError（）进行错误报告。另外，请注意，psz文件名实际上不一定是物理文件的名称（尽管通常是）。它的解释取决于驱动程序，它可能是一个URL，一个文件名，末尾添加了其他参数，控制打开或几乎任何东西。请尽量不要将GDAL文件选择对话框限制为仅选择物理文件。

获取数据集信息

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

在没有任何旋转或切变的“北上”影像的特殊但常见情况下，地理配准变换地理变换教程采用以下形式：

地理转换参数（Geotransform），是一个数组，包含了六个元素。
adfGeoTransform[0] /* top left x */                    
表示左上角象素点的 x 坐标。
adfGeoTransform[1] /* w-e pixel resolution */        
表示从左到右（也就是从西到东）每个象素点的宽度（水平分辨率）。   
adfGeoTransform[2] /* 0 */
表示轴的旋转角度，通常为 0。
adfGeoTransform[3] /* top left y */
表示左上角象素点的 y 坐标。
adfGeoTransform[4] /* 0 */
表示轴的旋转角度，通常为 0。
adfGeoTransform[5] /* n-s pixel resolution (negative value) */
表示从上到下（也就是从北到南）每个象素点的高度（垂直分辨率），是一个负值，因为地理坐标系是从上到下逐渐减小的。

在一般情况下，这是一个仿射变换。

如果我们想打印一些关于数据集的一般信息，我们可以执行以下操作：

在 Python 中：

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]))
这段代码主要是用来读取一个栅格数据集的基本信息，包括驱动、大小、投影以及地理转换等。

第一行代码中，通过 `dataset.GetDriver()` 方法获取数据集的驱动信息，
然后使用字符串格式化 `format()` 函数将短名称和长名称打印出来。

第二行代码中，使用 `dataset.RasterXSize`、`dataset.RasterYSize` 和 
`dataset.RasterCount` 分别获取栅格数据集的 X、Y 方向的大小和波段数，
然后使用字符串格式化打印出来。

第三行代码中，使用 `dataset.GetProjection()` 方法获取数据集的投影信息，
然后使用字符串格式化打印出来。

第四行代码中，使用 `dataset.GetGeoTransform()` 方法获取数据集的地理转换信息，
如果存在地理转换信息，就按照需要的格式打印出来。其中，`geotransform[0]` 和 
`geotransform[3]` 分别表示原点的左上角 X、Y 坐标，`geotransform[1]` 和 
`geotransform[5]` 分别表示像素在 X、Y 方向上的大小。

获取栅格波段

目前通过GDAL访问栅格数据是一次一个波段进行的。此外，元数据、块大小、颜色表和其他各种信息都可以在每个波段的基础上获得。以下代码从数据集中获取GDALRasterBand对象（编号从1到）并显示有关其的一些信息。GDALRasterBand :: GetRasterCount()。

在 Python 中：

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()))
这段代码用于获取和输出栅格数据集的属性信息。首先通过 `GetRasterBand()` 
方法获取数据集的第一个波段，然后通过 `GetDataTypeName()` 方法获取该波段
的数据类型并打印输出。接着使用 `GetMinimum()` 和 `GetMaximum()` 方法
获取该波段的最小值和最大值，如果这两个值任意一个为 None，则调用 
`ComputeRasterMinMax()` 方法计算栅格的最小值和最大值。最后分别输出
最小值和最大值。如果该波段有概览，输出概览数量；如果该波段有颜色表，输出颜色表的条目数。

读取栅格数据

有几种方法可以读取栅格数据，但最常见的是通过 GDALRasterBand：：RasterIO（） 方法。此方法将自动处理数据类型转换，向上/向下采样和窗口化。以下代码会将数据的第一扫描行读入大小相似的缓冲区，并在操作过程中将其转换为浮点数。

在 Python 中：

scanline = band.ReadRaster(xoff=0, yoff=0,
                        xsize=band.XSize, ysize=1,
                        buf_xsize=band.XSize, buf_ysize=1,
                        buf_type=gdal.GDT_Float32)
该代码使用 GDAL 库读取一个栅格数据集（例如遥感图像）中的一个波段，具体过程如下：

- 通过 `band` 变量获取一个波段对象，该对象包含有关该波段的元数据和像素值数据。
- 使用 `ReadRaster` 方法读取指定区域（即从横坐标为 0、纵坐标为 0 的像素开始，
横向跨越整个图像，纵向只有一个像素）内的像素值，并将这些值存储在一维数组 `scanline` 中。
- 具体参数含义如下：
  - `xoff`：指定读取像素起始位置的横坐标偏移量；
  - `yoff`：指定读取像素起始位置的纵坐标偏移量；
  - `xsize`：指定要读取的像素区域的宽度（即沿横向方向读取的像素数）；
  - `ysize`：指定要读取的像素区域的高度（即沿纵向方向读取的像素数）；
  - `buf_xsize`：指定存储读取结果的缓冲区（如数组或字符串）的宽度（即一维数组的长度）；
  - `buf_ysize`：指定存储读取结果的缓冲区的高度（即读取结果的维数，本例中为 1）；
  - `buf_type`：指定存储读取结果的缓冲区元素的数据类型（本例中为浮点型，即 `gdal.GDT_Float32`）。

请注意，返回的扫描线是字符串类型，包含 xsize*4 字节的原始二进制浮点数据。这可以使用标准库中的结构模块转换为 Python 值：

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

这段代码使用Python的标准库struct来解包二进制数据。具体来说，
它将一个字节数组scanline中的数据按照大小为4字节的浮点数格式
（'f'）进行解包，并将结果存储在一个元组tuple_of_floats中。
解包时使用了*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 )
这是GDALRasterBand类的一个方法，用于在栅格数据集上进行IO（input/output）操作。
该方法将指定区域的栅格数据从数据集中读取或写入到缓冲区。

参数说明：
- eRWFlag：读写标志，指定操作为读、写或读写。取值为GF_Read（读）、GF_Write（写）或GF_Update（读写）。
- nXOff/nYOff：指定读取或写入区域的左上角像素在数据集中的位置。
- nXSize/nYSize：指定需要读取或写入的像素数据集大小（宽度和高度）。
- pData：用于存储读取或写入的像素数据的缓冲区指针。
- nBufXSize/nBufYSize：缓冲区的大小，单位为像素。
- eBufType：缓冲区数据类型，与数据集中的数据类型进行匹配。
- nPixelSpace：缓冲区中相邻像素之间的偏移量（以字节为单位）。
- nLineSpace：缓冲区中相邻行之间的偏移量（以字节为单位）。

返回值为CPLErr类型，表示操作是否成功。如果成功，返回CE_None。如果发生错误，则返回表示错误类型的值。

请注意，相同的 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++ 删除运算符销毁它们。GDALDataset可以通过调用GDALClose（）来关闭（不建议在适用于Windows用户的GDALDataset上使用删除运算符，因为在跨模块边界分配和释放内存时存在已知问题）。请参阅常见问题解答中的相关主题）。调用 GDALClose 将导致正确清理和刷新任何挂起的写入。忘记在以更新模式以流行格式（如 GTiff）打开的数据集上调用 GDALClose 可能会导致之后无法打开它。

创建文件的技术

如果格式驱动程序支持创建，则可以创建 GDAL 支持格式的新文件。创建文件有两种通用技术，使用 CreateCopy（） 和 Create（）。CreateCopy 方法涉及在格式驱动程序上调用 CreateCopy（） 方法，并传入应复制的源数据集。Create 方法涉及在驱动程序上调用 Create（） 方法，然后使用单独的调用显式写入所有元数据和栅格数据。所有支持创建新文件的驱动程序都支持 CreateCopy（） 方法，但只有少数支持 Create（） 方法。

若要确定特定格式是否支持“创建”或“创建复制”，可以检查格式驱动程序对象上的DCAP_CREATE和DCAP_CREATECOPY元数据。在调用 GDALDriverManager：：GetDriverByName（） 之前，请确保已调用 GDALAllRegister（）。在这个例子中，我们获取一个驱动程序，并确定它是否支持 Create（） 和/或 CreateCopy（）。

在 Python 中：

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))

该代码使用Python中的GDAL库，主要功能是判断是否支持使用GDAL库中的Create()和
CreateCopy()方法来创建GDAL数据集。

首先，代码定义了一个变量fileformat表示文件格式，此处为GTiff格式。然后调用
gdal.GetDriverByName()方法将fileformat格式的驱动程序存储在driver变量中。
接着，调用driver.GetMetadata()方法获取该驱动程序的元数据信息，存储在metadata变量中。

接下来，代码针对获取到的元数据信息，使用metadata.get()方法查询驱动程序是否支持
Create()和CreateCopy()方法。若支持，则在控制台输出相应的信息，并将
fileformat格式作为参数传递。输出的信息为："Driver GTiff supports Create() method."
或"Driver GTiff supports CreateCopy() method."。


总结：该代码使用Python中的GDAL库来判断是否支持使用Create()和CreateCopy()方法创建GDAL数据集，
输出支持情况的信息提示。

请注意，许多驱动程序是只读的，不支持 Create（） 或 CreateCopy（）。

使用 CreateCopy（）

GDALDriver：：CreateCopy（） 方法的使用非常简单，因为大多数信息都是从源数据集中收集的。但是，它包括用于传递特定于格式的创建选项的选项，以及在进行长数据集复制时向用户报告进度的选项。从名为 pszSrcFilename 的文件到名为 pszDstFilename 的新文件的简单副本，使用先前获取其驱动程序的格式的默认选项可能如下所示：

在 Python 中：

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

这段代码使用了 Python 中的 GDAL 库来处理地理数据，具体解释如下：

1. `src_ds = gdal.Open(src_filename)`：打开名为 `src_filename` 
的地理数据文件，并将其存储为 GDAL 数据集对象 `src_ds`。
2. `dst_ds = driver.CreateCopy(dst_filename, src_ds, strict=0)`：
利用 GDAL 中提供的 `CreateCopy` 方法，基于 `src_ds` 创建一个全新的数据集 `dst_ds`，
并将其存储为名为 `dst_filename` 的新文件。这里的 `driver` 可以是 GDAL 支持的
不同地理文件格式的驱动器。`strict=0` 表示在创建新数据集时，忽略数据类型等方面的严格控制，
以提高灵活性。
3. `dst_ds = None` 和 `src_ds = None`：释放创建好的数据集对象，在代码结束前将其设为空值。
这样做可以释放内存空间，防止数据泄漏和造成程序错误。

请注意，CreateCopy（） 方法返回一个可写数据集，并且必须正确关闭它才能完成将数据集写入和刷新到磁盘。在 Python 的情况下，当“dst_ds”超出范围时，这会自动发生。在 CreateCopy（） 调用中，紧跟在目标文件名之后的 bStrict 选项使用的 FALSE（或 0）值表示，即使无法创建目标数据集以完全匹配输入数据集，CreateCopy（） 调用也应继续进行而不会发生致命错误。这可能是因为输出格式不支持输入数据集的像素数据类型，或者因为目标不支持例如写入地理配准。

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

在 Python 中：

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

这段代码首先使用Open函数打开一个名为src_filename的GDAL支持的栅格数据集，
并将其赋值给src_ds变量。

接下来，使用CreateCopy函数创建一个名为dst_filename
的栅格数据集，并将src_ds作为输入源，将新的数据集赋值给dst_ds变量。

此函数还可以接受可选参数，比如strict和options，这里设置了TILED和COMPRESS选项。
TILED选项表示将数据集存储为分块的方式，以提高IO效率，而COMPRESS选项表示
使用PACKBITS压缩算法对栅格数据进行压缩。

最后，我们将dst_ds和src_ds设置为None，以关闭和释放栅格数据集相关的资源。

使用 Create（）

对于不仅仅是将现有文件导出到新文件的情况，通常需要使用 GDALDriver：：Create（） 方法（尽管通过使用虚拟文件或内存中文件可以实现一些有趣的选项）。Create（） 方法采用的选项列表与 CreateCopy（） 非常相似，但必须明确提供图像大小、波段数和波段类型。

在 Python 中：

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

这段代码用于创建一个新的栅格数据集。其中，变量 `dst_filename` 是输出文件的路径和文件名；
`xsize` 和 `ysize` 分别是数据集的宽度和高度，单位为像素；`bands` 表示数据集中的波段数，
这里设置为1；`eType` 表示数据类型，这里设置为 `gdal.GDT_Byte` 表示数据为字节类型。
`driver` 是 GDAL 驱动程序实例，用于创建新的数据集。这段代码使用 `Create` 方法创建一个
具有指定属性的新栅格数据集，并返回一个指向该数据集的句柄。

成功创建数据集后，必须将所有适当的元数据和栅格数据写入文件。这将根据使用情况而有所不同，但此处介绍了投影、地理变换和栅格数据的简单情况。

在 Python 中：

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

这段代码主要是用来创建一个栅格数据集，并将其投影设置为 NAD27 UTM zone 11N。具体解释如下：

1. 首先导入 osr 模块并创建一个 SpatialReference 实例 srs。
2. 设置投影信息，包括地理转换（设置坐标系、偏移量和像素大小等）和投影方式（设置投影信息）。
3. 通过将 numpy 数组传递给数据集的 RasterBand 写入栅格数据。
4. 最后关闭数据集。

可以看出该段代码主要是在进行地图投影和栅格数据的写入，用于创建一个 NAD27 UTM zone 11N 坐标系下的栅格数据集。