根据烈度图生成危险概率分布图

import numpy as np
import rasterio
from rasterio.enums import Resampling
import matplotlib.pyplot as plt
import matplotlib

# 定义文件路径
intensity_file = 'intensity.tif'
base_file = 'base.tif'
output_file = 'risk_probability.tif'

# 风险概率计算函数
def calculate_risk_probability(intensity, base):
    return 1 / (1 + np.exp(0 - 0.00250165038535077 * np.exp(0.6931 * intensity) - base))

try:
    # 读取基础数据
    with rasterio.open(base_file) as src:
        base_data = src.read(1)  # 读取第一层数据
        base_meta = src.meta

    # 读取并重采样强度数据
    with rasterio.open(intensity_file) as src:
        intensity_data_resampled = src.read(
            out_shape=(1, base_data.shape[0], base_data.shape[1]),
            resampling=Resampling.bilinear
        )[0]  # 选择第一层数据

    print("Base data shape:", base_data.shape)
    print("Intensity data shape (after resampling):", intensity_data_resampled.shape)

    # 计算风险概率
    risk_probability = calculate_risk_probability(intensity_data_resampled, base_data)

    # 处理无效值
    risk_probability[np.isinf(risk_probability)] = 0
    risk_probability = np.clip(risk_probability, 0, 1)

    # 分级处理
    risk_classification = np.zeros_like(risk_probability)
    risk_classification[risk_probability >= 0.1] = 4  # 高风险
    risk_classification[(risk_probability >= 0.01) & (risk_probability < 0.1)] = 3  # 中高风险
    risk_classification[(risk_probability >= 0.001) & (risk_probability < 0.01)] = 2  # 中风险
    risk_classification[(risk_probability >= 0.0001) & (risk_probability < 0.001)] = 1  # 中低风险
    risk_classification[risk_probability < 0.0001] = 0  # 低风险

    # 可视化风险概率
    plt.figure(figsize=(10, 6))
    cmap = matplotlib.colormaps['RdYlBu_r']  # 使用反转的颜色映射
    plt.imshow(risk_classification, cmap=cmap, interpolation='nearest', extent=(base_meta['transform'][2], 
                                                                             base_meta['transform'][2] + base_meta['width'] * base_meta['transform'][0], 
                                                                             base_meta['transform'][5] + base_meta['height'] * base_meta['transform'][4], 
                                                                             base_meta['transform'][5]))

    # 添加经纬度刻度
    plt.colorbar(ticks=[0, 1, 2, 3, 4], label='风险等级')
    plt.title('次生地质灾害风险概率分级图')
    plt.xlabel('经度')
    plt.ylabel('纬度')
    plt.axis('on')  # 显示坐标轴
    plt.show()

    # 将风险概率写入新的 TIFF 文件
    risk_probability = (risk_probability * 255).astype(np.uint8)  # 归一化到 0-255
    with rasterio.open(output_file, 'w', driver='GTiff',
                       height=risk_probability.shape[0],
                       width=risk_probability.shape[1],
                       count=1,
                       dtype='uint8',
                       crs=base_meta['crs'],
                       transform=base_meta['transform']) as dst:
        dst.write(risk_probability, 1)

    print(f"风险概率地图已保存至 {output_file}")

except Exception as e:
    print(f"发生错误: {e}")

base_meta
以下是对从 GeoTIFF 文件中提取的属性的解释：g
'driver': 'GTiff':
解释：这表明文件格式是 GeoTIFF（GTiff 是 GeoTIFF 的缩写）。GeoTIFF 是一种广泛使用的地理空间栅格数据格式，它扩展了 TIFF（标签图像文件格式）以支持地理坐标系统信息。
'dtype': 'float32':
解释：这表示栅格数据中的每个像素值的数据类型是 32 位浮点数。这意味着每个像素可以存储一个浮点数值，其范围大约在 -3.4e+38 到 3.4e+38 之间。
'nodata': -3.4028230607370965e+38:
解释：这是定义的无数据值，用于表示栅格数据中缺失或无效的像素。在这个例子中，nodata 值设置为一个非常小的负数，这是 float32 类型所能表示的最小数值之一。
'width': 1822, 'height': 588:
解释：这些属性表示栅格数据的尺寸。width 是图像的宽度（以像素为单位），而 height 是图像的高度（以像素为单位）。所以这个 GeoTIFF 文件包含一个 1822 像素宽和 588 像素高的图像。
'count': 1:
解释：这表示栅格数据集包含一个波段（band）。如果是彩色图像，count 可能是 3（红、绿、蓝波段）或者更多。
'crs': CRS.from_wkt('...'):
解释：这是栅格数据的坐标参考系统（CRS）。CRS 定义了数据在地球上的位置和如何将像素坐标转换为地理坐标。这里的 WKT（Well-Known Text）字符串定义了一个地理坐标系 WGS 84，这是全球定位系统（GPS）使用的坐标系。
'transform': Affine(0.008333333300000006, 0.0, 89.82031276588475, 0.0, -0.0083333333, 32.3978756035457):
解释：这是一个仿射变换矩阵，用于将像素坐标转换为地理坐标。这个 Affine 对象包含六个参数，分别是：
a (x-scale): 0.008333333300000006，表示像素在 x 方向上的大小（经度方向上的分辨率）。
b (x-skew): 0.0，表示像素在 x 方向上的倾斜度（通常为 0）。
c (x-translate): 89.82031276588475，表示图像左上角像素的 x 坐标（经度）。
d (y-skew): 0.0，表示像素在 y 方向上的倾斜度（通常为 0）。
e (y-scale): -0.0083333333，表示像素在 y 方向上的大小（纬度方向上的分辨率，通常是负值）。
f (y-translate): 32.3978756035457，表示图像左上角像素的 y 坐标（纬度）。
这个仿射变换矩阵用于计算给定像素坐标 (x, y) 的地理坐标 (longitude, latitude)。