整合tif图与geojson数据

import numpy as np
import rasterio
import rasterio.plot
from rasterio.mask import mask
from rasterio.warp import reproject, calculate_default_transform, Resampling as WarpResampling
import geopandas as gpd
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from shapely.geometry import box
import pandas as pd
import matplotlib

# 设置中文字体（需要系统中有相应的字体）
matplotlib.rcParams['font.sans-serif'] = ['SimHei']  # 使用黑体
matplotlib.rcParams['axes.unicode_minus'] = False    # 正常显示负号

# 1. 定义函数1：读取 base.tif 文件
def read_base_tif(file_path):
    with rasterio.open(file_path) as src:
        data = src.read(1)  # 读取第一波段
        profile = src.profile.copy()  # 获取元数据
        transform = src.transform  # 地理变换信息
        crs = src.crs  # 坐标参考系
        bounds = src.bounds  # 图像范围
    return data, profile, transform, crs, bounds

# 2. 定义函数2：读取 intensity.tif 文件
def read_intensity_tif(file_path):
    with rasterio.open(file_path) as src:
        data = src.read(1)
        profile = src.profile.copy()
        transform = src.transform
        crs = src.crs
        bounds = src.bounds
    return data, profile, transform, crs, bounds

# 3. 读取并合并 GeoJSON 文件
def read_geojson_files(file_list):
    gdf_list = []
    for file in file_list:
        gdf = gpd.read_file(file)
        gdf_list.append(gdf)
    combined_gdf = pd.concat(gdf_list, ignore_index=True)
    combined_gdf = gpd.GeoDataFrame(combined_gdf)
    return combined_gdf

# 5. 截取 base.tif 与 intensity.tif 重叠的部分
def extract_overlapping_region(base_data, base_transform, base_crs, base_bounds,
                               intensity_data, intensity_transform, intensity_crs, intensity_bounds):
    # 确保坐标参考系一致
    if base_crs != intensity_crs:
        raise ValueError("base.tif 和 intensity.tif 的坐标参考系不一致，请进行坐标转换。")

    # 计算重叠区域
    base_geom = box(*base_bounds)
    intensity_geom = box(*intensity_bounds)
    overlap_geom = base_geom.intersection(intensity_geom)

    if overlap_geom.is_empty:
        raise ValueError("base.tif 和 intensity.tif 没有重叠区域。")

    # 将重叠区域转换为 GeoJSON 格式
    overlap_geojson = [overlap_geom.__geo_interface__]

    # 裁剪 base_data
    with rasterio.open('base.tif') as src:
        base_overlap_data, base_overlap_transform = mask(src, overlap_geojson, crop=True)
        base_overlap_meta = src.meta.copy()
        base_overlap_meta.update({
            "height": base_overlap_data.shape[1],
            "width": base_overlap_data.shape[2],
            "transform": base_overlap_transform
        })

    # 裁剪 intensity_data
    with rasterio.open('intensity.tif') as src:
        intensity_overlap_data, intensity_overlap_transform = mask(src, overlap_geojson, crop=True)
        intensity_overlap_meta = src.meta.copy()
        intensity_overlap_meta.update({
            "height": intensity_overlap_data.shape[1],
            "width": intensity_overlap_data.shape[2],
            "transform": intensity_overlap_transform
        })

    # 重采样 base_overlap_data 以匹配 intensity_overlap_data
    base_overlap_resampled = np.empty_like(intensity_overlap_data, dtype=base_overlap_data.dtype)

    reproject(
        source=base_overlap_data,
        destination=base_overlap_resampled,
        src_transform=base_overlap_transform,
        src_crs=base_crs,
        dst_transform=intensity_overlap_transform,
        dst_crs=intensity_crs,
        resampling=WarpResampling.nearest
    )

    return base_overlap_resampled[0], intensity_overlap_data[0], intensity_overlap_transform, intensity_crs

# 6. 根据给定公式计算风险概率
def calculate_risk_probability(intensity, base):
    # intensity 是仪器烈度，如 6、7、8 等
    # base 是底图层
    # 计算 A = 0.00250165038535077 * e^(0.6931 * intensity)
    A = 0.00250165038535077 * np.exp(0.6931 * intensity)

    # 计算指数 exponent = - (A + base)
    exponent = - (A + base)

    # 计算风险概率 result = 1 / (1 + e^(exponent))
    with np.errstate(over='ignore', invalid='ignore'):
        result = 1 / (1 + np.exp(exponent))

    result = np.nan_to_num(result, nan=0.0, posinf=0.0, neginf=0.0)
    return result

# 7. 根据概率值对风险级别进行分类
def classify_risk_levels(probability):
    risk_levels = np.zeros_like(probability, dtype=np.uint8)

    # 风险级别从 1 到 5，对应五种颜色，1 表示低风险
    risk_levels[probability >= 0.1] = 5  # 高风险（红色）
    risk_levels[(probability >= 0.01) & (probability < 0.1)] = 4  # 中高风险（橙色）
    risk_levels[(probability >= 0.001) & (probability < 0.01)] = 3  # 中风险（黄色）
    risk_levels[(probability >= 0.0001) & (probability < 0.001)] = 2  # 中低风险（蓝色）
    risk_levels[probability < 0.0001] = 1  # 低风险（白色）

    return risk_levels

# 8. 叠加区县信息并可视化
def visualize_risk(risk_levels, transform, crs, combined_gdf):
    # 生成颜色映射
    colors = ['white', 'blue', 'yellow', 'orange', 'red']
    cmap = ListedColormap(colors)

    fig, ax = plt.subplots(figsize=(10, 8))

    # 显示风险等级图
    rasterio.plot.show(risk_levels, transform=transform, cmap=cmap, ax=ax)
    ax.set_title("次生地质灾害风险概率图")
    ax.set_xlabel('经度')
    ax.set_ylabel('纬度')

    # 添加图例
    labels = ['低风险', '中低风险', '中风险', '中高风险', '高风险']
    patches = [plt.plot([], [], marker="s", ms=10, ls="", mec=None, color=colors[i], label="{:s}".format(labels[i]))[0] for i in range(len(colors))]
    ax.legend(handles=patches, bbox_to_anchor=(1.05, 1), loc='upper left')

    # 确保坐标参考系一致
    if combined_gdf.crs != crs:
        combined_gdf = combined_gdf.to_crs(crs)

    # 裁剪区县数据到重叠区域
    width = risk_levels.shape[1]
    height = risk_levels.shape[0]
    left, top = transform * (0, 0)
    right, bottom = transform * (width, height)
    overlap_geom = box(left, bottom, right, top)
    overlap_gdf = gpd.GeoDataFrame(geometry=[overlap_geom], crs=crs)
    clipped_gdf = gpd.overlay(combined_gdf, overlap_gdf, how='intersection')

    # 绘制区县边界
    clipped_gdf.boundary.plot(ax=ax, edgecolor='black', linewidth=0.5)

    # 添加区县名称标签
    for idx, row in clipped_gdf.iterrows():
        centroid = row['geometry'].centroid
        ax.text(centroid.x, centroid.y, row.get('name', '未知'), fontsize=6, ha='center', va='center', color='black')

    plt.show()

    # 保存结果为 TIFF 文件
    plt.savefig('final_result.tif', dpi=300, bbox_inches='tight', format='tif')
    print("风险概率图已成功保存为 'final_result.tif'。")



# 9. 主程序
if __name__ == "__main__":
    # 读取 base.tif 文件
    base_data, base_profile, base_transform, base_crs, base_bounds = read_base_tif('base.tif')

    # 读取 intensity.tif 文件
    intensity_data, intensity_profile, intensity_transform, intensity_crs, intensity_bounds = read_intensity_tif('intensity.tif')

    # 读取 GeoJSON 文件
    sichuan_json_files = ['阿坝藏族羌族自治州.json', '成都市.json', '德阳市.json', '甘孜藏族自治州.json',
                          '乐山市.json', '凉山彝族自治州.json', '眉山市.json', '绵阳市.json',
                          '内江市.json', '遂宁市.json', '雅安市.json', '宜宾市.json',
                          '资阳市.json', '自贡市.json']
    xizang_json_files = ['昌都市.json', '那曲市.json', '日喀则市.json', '山南市.json', '林芝市.json']
    all_json_files = sichuan_json_files + xizang_json_files

    combined_gdf = read_geojson_files(all_json_files)

    # 截取重叠区域
    base_overlap_data, intensity_overlap_data, overlap_transform, overlap_crs = extract_overlapping_region(
        base_data, base_transform, base_crs, base_bounds,
        intensity_data, intensity_transform, intensity_crs, intensity_bounds
    )

    # 计算风险概率
    risk_probability = calculate_risk_probability(intensity_overlap_data, base_overlap_data)

    # 对结果进行分级
    risk_levels = classify_risk_levels(risk_probability)

    # 保存风险等级为 TIFF 文件
    out_meta = base_profile.copy()
    out_meta.update({
        "driver": "GTiff",
        "height": risk_levels.shape[0],
        "width": risk_levels.shape[1],
        "transform": overlap_transform,
        "crs": overlap_crs,
        "count": 1,
        "dtype": 'uint8',
        "nodata": 0
    })
    with rasterio.open('risk_level.tif', 'w', **out_meta) as dest:
        dest.write(risk_levels, 1)

    # 可视化结果并叠加区县信息
    visualize_risk(risk_levels, overlap_transform, overlap_crs, combined_gdf)