python快速启动¶

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

以读取模式打开数据集¶

考虑一个名为 example.tif 16位陆地卫星8图像覆盖美国科罗拉多高原的一部分 1. 因为图像很大（70 MB），并且动态范围很宽，所以很难在浏览器中显示。重新缩放和动态压缩的版本如下所示。

导入栅格开始。

>>> import rasterio

接下来，打开文件。

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

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

>>> dataset.name
'example.tif'
>>> dataset.mode
'r'
>>> dataset.closed
False

数据集属性¶

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

>>> dataset.count
1

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

>>> dataset.width
7731
>>> dataset.height
7871

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

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

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

数据集地理参考¶

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

>>> dataset.bounds
BoundingBox(left=358485.0, bottom=4028985.0, right=590415.0, top=4265115.0)

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

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

>>> dataset.transform
Affine(30.0, 0.0, 358485.0,
       0.0, -30.0, 4265115.0)

数据集的 transform 是将（行、列）坐标中的像素位置映射到（x、y）空间位置的仿射变换矩阵。这个矩阵的乘积 (0, 0) 数据集左上角的行和列坐标是左上角的空间位置。

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

同样获得右下角的位置。

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

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

>>> dataset.crs
CRS.from_epsg(32612)

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

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

读取栅格数据¶

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

>>> dataset.indexes
(1,)
>>> 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=uint16)

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

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

空间索引¶

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

>>> x, y = (dataset.bounds.left + 100000, dataset.bounds.top - 50000)
>>> row, col = dataset.index(x, y)
>>> row, col
(1666, 3333)
>>> band1[row, col]
7566

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

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

正在创建数据¶

阅读数据只是故事的一半。使用栅格数据集对象，可以将值数组写入栅格数据文件，从而与其他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.033333333333333333, 0.0, -4.0166666666666666,
       0.0, -0.033333333333333333, 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格式可以。

保存栅格数据¶

若要将网格复制到打开的数据集，请调用新数据集的 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 部分。

1: “example.tif”是陆地卫星场景LC80370342016194LGN00第4波段的别名。