rasterio栅格数据分析代码

下面是使用Rasterio进行栅格数据分析的示例代码：

import rasterio
import numpy as np
import matplotlib.pyplot as plt

# 打开栅格数据
with rasterio.open('path/to/raster.tif') as src:
    # 读取数据和元数据
    data = src.read()
    meta = src.meta

# 查看数据的维度、数据类型和范围
print('Data shape:', data.shape)
print('Data type:', data.dtype)
print('Data range:', np.min(data), np.max(data))

# 可视化数据
plt.imshow(data[0], cmap='gray')
plt.show()

# 重投影数据
from rasterio.warp import calculate_default_transform, reproject, Resampling

dst_crs = 'EPSG:4326'
transform, width, height = calculate_default_transform(meta['crs'], dst_crs, meta['width'], meta['height'], *meta['bounds'])
meta.update({'crs': dst_crs, 'transform': transform, 'width': width, 'height': height})

with rasterio.open('path/to/reprojected.tif', 'w', **meta) as dst:
    reproject(
        source=rasterio.band(data, 1),
        destination=rasterio.band(dst, 1),
        src_transform=meta['transform'],
        src_crs=meta['crs'],
        dst_transform=transform,
        dst_crs=dst_crs,
        resampling=Resampling.nearest
    )

# 剪切数据
from rasterio.mask import mask

# 定义裁剪区域
geometry = {'type': 'Polygon', 'coordinates': [[(x1, y1), (x2, y2), (x3, y3), (x4, y4)]]}

# 根据裁剪区域获取掩膜
mask = rasterio.features.geometry_mask([geometry], out_shape=data.shape, transform=meta['transform'], invert=True)

# 剪切数据
clipped_data = data[:, mask]

# 可视化剪切后的数据
plt.imshow(clipped_data[0], cmap='gray')
plt.show()

这段代码首先使用`rasterio.open()`函数打开一个栅格数据，然后读取数据和元数据。接着，我们查看了数据的维度、数据类型和范围，并可视化了数据。然后，我们使用`rasterio.warp.reproject()`函数将数据重投影到一个新的坐标系，并保存为一个新的栅格数据。最后，我们使用`rasterio.mask.mask()`函数根据一个裁剪区域获取掩膜，并用掩膜剪切了数据。