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

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

>>> dataset.width
500
>>> dataset.height
300

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

>>> {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]
0

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]
0

要获取像素的空间坐标,请使用数据集的 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 . 其轮廓如下所示。

_images/field.png

以写入模式打开数据集

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

  1. 驱动 :所需格式驱动程序的名称

  2. 宽度 :数据集的列数

  3. 高度 :数据集的行数

  4. 计数 :数据集带区计数

  5. D型 :数据集的数据类型

  6. crs :坐标参考系标识符或描述

  7. 转型 :仿射变换矩阵,以及

  8. 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 部分。