Rasterio快速入门：Python地理空间数据处理指南

Rasterio快速入门：Python地理空间数据处理指南

【免费下载链接】rasterio 项目地址: https://gitcode.com/gh_mirrors/ras/rasterio

前言

Rasterio是一个用于处理地理空间栅格数据的Python库，它基于GDAL构建，提供了简洁高效的API接口。本文将带你快速了解如何使用Rasterio进行栅格数据的读写操作，适合有一定Python基础但初次接触地理空间数据处理的开发者。

环境准备

在开始之前，请确保你已经安装了Rasterio及其依赖项。本文所有示例都基于GeoTIFF格式，但这些概念同样适用于其他栅格数据格式。

读取栅格数据

打开数据集

读取栅格数据的第一步是打开数据集文件。Rasterio的open()函数接受文件路径字符串或类路径对象，返回一个数据集对象。

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

数据集对象具有类似Python文件对象的属性：

print(dataset.name)  # 输出文件名
print(dataset.mode)  # 文件访问模式
print(dataset.closed)  # 文件是否关闭

数据集属性

栅格数据集包含多个重要属性，可以帮助我们理解数据内容：

1. 波段数量：dataset.count返回数据集的波段数
2. 数据维度：dataset.width和dataset.height分别返回栅格的列数和行数
3. 数据类型：dataset.dtypes返回各波段的数据类型

print(f"波段数量: {dataset.count}")
print(f"栅格宽度: {dataset.width}, 高度: {dataset.height}")
print(f"数据类型: {dataset.dtypes}")

地理参考信息

地理空间数据与普通图像的关键区别在于其像素与地球表面的映射关系：

1. 空间范围：dataset.bounds返回数据集的边界框
2. 坐标变换：dataset.transform是仿射变换矩阵，将像素坐标转换为空间坐标
3. 坐标参考系统(CRS)：dataset.crs定义了数据集使用的坐标系统

print(f"空间范围: {dataset.bounds}")
print(f"坐标变换矩阵: {dataset.transform}")
print(f"坐标参考系统: {dataset.crs}")

读取数据

使用read()方法可以获取栅格波段数据，返回的是NumPy数组：

band1 = dataset.read(1)  # 读取第一个波段
print(band1)  # 输出数组内容

空间索引操作

Rasterio提供了便捷的空间索引方法：

1. 坐标转像素：dataset.index(x, y)将地理坐标转换为行列索引
2. 像素转坐标：dataset.xy(row, col)将行列索引转换为地理坐标

# 地理坐标转像素索引
x, y = (dataset.bounds.left + 100000, dataset.bounds.top - 50000)
row, col = dataset.index(x, y)

# 像素索引转地理坐标
center_x, center_y = dataset.xy(dataset.height // 2, dataset.width // 2)

创建和写入栅格数据

准备数据

首先创建一个示例数据数组，这里我们生成一个包含两个高斯分布差异的二维数组：

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)

创建数据集

写入数据时需要指定多个关键参数：

from rasterio.transform import Affine

# 计算分辨率和变换矩阵
res = (x[-1] - x[0]) / 240.0
transform = Affine.translation(x[0] - res / 2, y[0] - res / 2) * Affine.scale(res, res)

# 使用上下文管理器创建数据集
with rasterio.open(
    'output.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)  # 写入数据到第一个波段

最佳实践

1. 使用上下文管理器：确保文件正确关闭，避免资源泄漏
2. 检查数据类型：确保写入的数据类型与声明的dtype匹配
3. 验证变换矩阵：确保地理参考信息准确无误
4. 处理大文件：对于大型数据集，考虑分块读写以提高效率

总结

本文介绍了Rasterio的基本使用方法，包括：

• 如何打开和读取现有栅格数据集
• 如何访问和解释地理参考信息
• 如何进行空间坐标与像素索引的转换
• 如何创建新的栅格数据集并写入数据

掌握了这些基础知识后，你就可以开始使用Python处理各种地理空间栅格数据了。Rasterio还提供了更多高级功能，如波段运算、重采样、数据融合等，值得进一步探索。

【免费下载链接】rasterio 项目地址: https://gitcode.com/gh_mirrors/ras/rasterio