导入栅格开始。
import rasterio
接下来，打开文件。
dataset = rasterio.open('/data/gdata/geotiff_file.tif')
rasterio open() 函数接受一个路径字符串或类似路径的对象，并返回一个打开的数据集对象。
路径可以指向任何支持的栅格格式的文件。
Rasterio将使用适当的gdal格式驱动程序打开它。数据集对象具有与Python文件对象相同的一些属性。
dataset.name
'/data/gdata/geotiff_file.tif'
dataset.mode
'r'
dataset.closed
False
dataset.count
3
数据集带是一个值数组，表示二维空间中单个变量的部分分布。 数据集的所有带区数组的行数和列数都相同。 示例数据集唯一波段表示的变量是陆地卫星8号操作陆地成像仪（OLI）波段4（波长在640-670纳米之间）的1级数字（DN）。 这些值可以缩放为辐射值或反射值。
dataset.width
1500
dataset.height
900
一些数据集属性通过一个值元组（每个值带一个）公开所有数据集带区的属性。
要获取带索引到变量数据类型的映射，请将字典理解应用于 zip() 数据集的乘积 indexes 和 dtypes 属性。
{i: dtype for i, dtype in zip(dataset.indexes, dataset.dtypes)}
{1: 'uint8', 2: 'uint8', 3: 'uint8'}
示例文件的唯一带区包含无符号16位整数值。geotiff格式还支持不同大小的有符号整数和浮点。
dataset.bounds
BoundingBox(left=1868454.913, bottom=5326126.266, right=1913454.913, top=5353126.266)
我们的示例覆盖了从358485米（在本例中）到590415米（从左到右）以及从下到上4028985米到4265115米的世界。 它覆盖了一个区域，宽231.93公里，高236.13公里。
bounds 属性的值来自一个更基本的属性：数据集的地理空间转换。
dataset.transform
Affine(30.0, 0.0, 1868454.913,
0.0, -30.0, 5353126.266)
数据集的 transform 是将（行、列）坐标中的像素位置映射到 (x, y) 空间位置的仿射变换矩阵。
以下的乘积 (0, 0) 数据集左上角的行和列坐标是左上角的空间位置。
dataset.transform * (0, 0)
(1868454.913, 5353126.266)
同样获得右下角的位置。
dataset.transform * (dataset.width, dataset.height)
(1913454.913, 5326126.266)
但这些数字是什么意思？离哪里4028985米？这些坐标值相对于数据集坐标参考系（CRS）的原点。
dataset.crs
CRS.from_wkt('PROJCS["Albers_Beijing54",GEOGCS["Unknown datum based upon the Krassowsky 1940 ellipsoid",DATUM["Not_specified_based_on_Krassowsky_1940_ellipsoid",SPHEROID["Krassowsky 1940",6378245,298.3,AUTHORITY["EPSG","7024"]],AUTHORITY["EPSG","6024"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["latitude_of_center",0],PARAMETER["longitude_of_center",105],PARAMETER["standard_parallel_1",25],PARAMETER["standard_parallel_2",47],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]')
dataset.indexes
(1, 2, 3)
band1 = dataset.read(1)
这个 read() 方法返回一个numpy n-d数组。
band1
array([[145, 146, 145, ..., 67, 43, 28],
[158, 157, 155, ..., 74, 53, 35],
[156, 157, 156, ..., 52, 42, 28],
...,
[223, 220, 219, ..., 24, 3, 0],
[222, 219, 219, ..., 76, 41, 56],
[221, 218, 218, ..., 70, 64, 86]], dtype=uint8)
数组中的值可以通过行、列索引寻址。
band1[dataset.height // 2, dataset.width // 2]
np.uint8(55)
x, y = (dataset.bounds.left + 100, dataset.bounds.top - 50)
row, col = dataset.index(x, y)
row, col
(1, 3)
band1[row, col]
np.uint8(156)
要获取像素的空间坐标，请使用数据集的 xy() 方法。图像中心的坐标可以这样计算。
dataset.xy(dataset.height // 2, dataset.width // 2)
(np.float64(1890969.913), np.float64(5339611.266))
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(np.float64(0.03333333333333333), np.float64(0.0), np.float64(-4.016666666666667),
np.float64(0.0), np.float64(-0.03333333333333333), np.float64(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,
)
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 部分。