>>> from env_helper import info; info()

页面更新时间： 2024-04-10 14:53:39
运行环境：
    Linux发行版本: Debian GNU/Linux 12 (bookworm)
    操作系统内核: Linux-6.1.0-18-amd64-x86_64-with-glibc2.36
    Python版本: 3.11.2

3.2. 开始使用¶

读写数据文件是空间数据程序员的工作。本文档介绍如何使用栅格读取现有文件和创建新文件。一些高级主题会被掩盖，以便在Rasterio文档的其他地方更详细地介绍。这里只使用 GeoTiff格式，但示例同样适用于其他栅格数据格式。

3.2.1. 以读取模式打开数据集¶

导入栅格开始。

>>> import rasterio

接下来，打开文件。

>>> dataset = rasterio.open('/gdata/rasterio/example.tif')

/usr/lib/python3/dist-packages/rasterio/__init__.py:304: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
  dataset = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)

rasterio open() 函数接受一个路径字符串或类似路径的对象，并返回一个打开的数据集对象。路径可以指向任何支持的栅格格式的文件。 Rasterio将使用适当的gdal格式驱动程序打开它。数据集对象具有与Python文件对象相同的一些属性。

>>> dataset.name

'/gdata/rasterio/example.tif'

>>> dataset.mode

'r'

>>> dataset.closed

False

3.2.2. 数据集属性¶

存储在示例geotiff中的栅格数据的属性可以通过打开的数据集对象的属性访问。数据集对象有条带，此示例的条带计数为1。

>>> dataset.count

数据集带是一个值数组，表示二维空间中单个变量的部分分布。数据集的所有带区数组的行数和列数都相同。示例数据集唯一波段表示的变量是陆地卫星8号操作陆地成像仪（OLI）波段4（波长在640-670纳米之间）的1级数字（DN）。这些值可以缩放为辐射值或反射值。

>>> dataset.width

>>> dataset.height

一些数据集属性通过一个值元组（每个值带一个）公开所有数据集带区的属性。要获取带索引到变量数据类型的映射，请将字典理解应用于 zip() 数据集的乘积 indexes 和 dtypes 属性。

>>> {i: dtype for i, dtype in zip(dataset.indexes, dataset.dtypes)}

{1: 'uint8', 2: 'uint8', 3: 'uint8'}

示例文件的唯一带区包含无符号16位整数值。geotiff格式还支持不同大小的有符号整数和浮点。

3.2.3. 数据集地理参考¶

地理信息系统栅格数据集不同于普通图像；其元素（或“像素”）映射到地球表面的区域。数据集的每个像素都包含在空间边界框中。

>>> dataset.bounds

BoundingBox(left=0.0, bottom=300.0, right=500.0, top=0.0)

我们的示例覆盖了从358485米（在本例中）到590415米（从左到右）以及从下到上4028985米到4265115米的世界。它覆盖了一个区域，宽231.93公里，高236.13公里。

bounds 属性的值来自一个更基本的属性：数据集的地理空间转换。

>>> dataset.transform

Affine(1.0, 0.0, 0.0,
       0.0, 1.0, 0.0)

数据集的 transform 是将（行、列）坐标中的像素位置映射到 (x, y) 空间位置的仿射变换矩阵。

以下的乘积 (0, 0) 数据集左上角的行和列坐标是左上角的空间位置。

>>> dataset.transform * (0, 0)

(0.0, 0.0)

同样获得右下角的位置。

>>> dataset.transform * (dataset.width, dataset.height)

(500.0, 300.0)

但这些数字是什么意思？离哪里4028985米？这些坐标值相对于数据集坐标参考系（CRS）的原点。

>>> dataset.crs

“epsg 32612”表示特定的坐标参考系： UTM 12N区。该系统用于在西108度到114度之间绘制北半球的区域图。示例数据集的左上角， (358485.0, 4265115.0) 位于12区中央子午线以西141.5公里（西111度），赤道以北4265公里。

在 crs 属性与 transform 描述了栅格数据集的地理参考，并将其与其他地理信息系统数据集进行了比较。

3.2.4. 读取栅格数据¶

栅格条带的数据可以通过条带的索引号访问。按照gdal约定，条带从1开始索引。

>>> dataset.indexes

(1, 2, 3)

>>> band1 = dataset.read(1)

这个 read() 方法返回一个numpy n-d数组。

>>> band1

array([[0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       ...,
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0]], dtype=uint8)

数组中的值可以通过行、列索引寻址。

>>> band1[dataset.height // 2, dataset.width // 2]

3.2.5. 空间索引¶

数据集具有 index() 获取地理参考空间中点对应的数组索引的方法。要获取数据集左上角以东100公里和以南50公里的像素值，请执行以下操作。

>>> x, y = (dataset.bounds.left + 100, dataset.bounds.top - 50)
>>> row, col = dataset.index(x, y)
>>> row, col

(-50, 100)

>>> band1[row, col]

要获取像素的空间坐标，请使用数据集的 xy() 方法。图像中心的坐标可以这样计算。

>>> dataset.xy(dataset.height // 2, dataset.width // 2)

(250.5, 150.5)

3.2.6. 创建数据¶

读取数据只是故事的一半。使用栅格数据集对象，可以将值数组写入栅格数据文件，从而与其他GIS应用程序（如QGIS）共享。

作为一个例子，考虑一组浮点值，代表温度或压力异常场，该场在规则网格上测量或建模，240列，180行。横轴上的第一个和最后一个网格点位于西4.0度和东经4.0度，纵轴上的第一个和最后一个网格点位于南纬3度和北纬3度。

>>> import numpy as np
>>> x = np.linspace(-4.0, 4.0, 240)
>>> y = np.linspace(-3.0, 3.0, 180)
>>> X, Y = np.meshgrid(x, y)
>>> Z1 = np.exp(-2 * np.log(2) * ((X - 0.5) ** 2 + (Y - 0.5) ** 2) / 1 ** 2)
>>> Z2 = np.exp(-3 * np.log(2) * ((X + 0.5) ** 2 + (Y + 0.5) ** 2) / 2.5 ** 2)
>>> Z = 10.0 * (Z2 - Z1)

本例中的虚拟字段由两个高斯分布的差组成，并由数组表示。 Z . 其轮廓如下所示。

以写入模式打开数据集

要将此阵列以及地理参考信息保存到新的栅格数据文件，请调用 rasterio.open() 使用要创建的新文件的路径， ‘w’ 指定写入模式和几个关键字参数。

驱动：所需格式驱动程序的名称
宽度：数据集的列数
高度：数据集的行数
计数：数据集带区计数
D型：数据集的数据类型
crs ：坐标参考系标识符或描述
转型：仿射变换矩阵，以及
NODATA ：一个“nodata”值

这些关键字参数的前5个参数对数据文件的特定格式属性进行参数化，并且在打开要写入的文件时是必需的。最后3个是可选的。

在这个例子中，坐标参考系是 ‘+proj=latlong’ 它描述了一个以十进制度数为单位的等矩形坐标参考系。右仿射变换矩阵可以用 from_origin() .

>>> from rasterio.transform import from_origin
>>> res = (x[-1] - x[0]) / 240.0
>>> transform = from_origin(x[0] - res / 2, y[-1] + res / 2, res, res)
>>> transform

Affine(0.03333333333333333, 0.0, -4.016666666666667,
       0.0, -0.03333333333333333, 3.0166666666666666)

示例网格中的左上角位于西3度和北2度。以该网格点为中心的栅格像素将延伸 res / 2 或1/60度，在每个方向，因此在上述表达式中的移动。

用于存储示例网格的数据集是这样打开的

>>> new_dataset = rasterio.open(
>>>      '/tmp/new.tif',
>>>      'w',
>>>      driver='GTiff',
>>>      height=Z.shape[0],
>>>      width=Z.shape[1],
>>>      count=1,
>>>      dtype=Z.dtype,
>>>      crs='+proj=latlong',
>>>      transform=transform,
>>>  )

的值高度，宽度和 D型关键字参数直接取自二维数组的属性， Z . 并非所有栅格格式都支持中的64位浮点值。 Z 但geotiff格式可以。

3.2.7. 保存栅格数据¶

若要将网格复制到打开的数据集，请调用新数据集的 write() 方法，将网格和目标带号作为参数。

>>> new_dataset.write(Z, 1)

然后调用 close() 方法将数据同步到磁盘并完成。

>>> new_dataset.close()

因为Rasterio的数据集对象模仿了Python的文件对象并实现了Python的上下文管理器协议，所以可以改为执行以下操作。

>>> with rasterio.open(
>>>     '/tmp/new.tif',
>>>     'w',
>>>     driver='GTiff',
>>>     height=Z.shape[0],
>>>     width=Z.shape[1],
>>>     count=1,
>>>     dtype=Z.dtype,
>>>     crs='+proj=latlong',
>>>     transform=transform,
>>> ) as dst:
>>>     dst.write(Z, 1)

这些是读取和写入栅格数据文件的基础。更多功能和示例包含在 advanced topics 部分。

3.1. RaterIO介绍

3.3. NumPy 掩膜

Python与开源GIS