【数据集】2000-2020年中国4 km日尺度气象数据集 (CDMet)

    import os
    import numpy as np
    import geopandas as gpd
    from netCDF4 import Dataset
    import matplotlib.pyplot as plt
    
    # 设置全局字体为 Times New Roman
    plt.rcParams["font.family"] = "Times New Roman"
    
    # 数据文件夹路径
    data_dir = r"C:\Database\3_CDMet\Relative humidity"  # 替换为您的文件夹路径
    variable_name = "rhu"  # 相对湿度变量名
    
    # SHP 文件路径
    Chinashp_file_path = "C:/PycharmProjects/GBA_Data_Process/shpFile/China/中国.shp"
    
    # 主程序：加载站点数据和 SHP 文件
    gdf_China = gpd.read_file(Chinashp_file_path)
    
    def read_daily_rhu(data_dir, year, day, lat_range=None, lon_range=None):
    """
    从 NetCDF 文件中读取指定年份和某一天的湿度数据，并应用缩放和无效值处理。
    可根据给定的经纬度范围裁剪数据。
    :param data_dir: 包含 NetCDF 文件的文件夹路径
    :param year: 年份 (如 2020)
    :param day: 天索引 (0 表示 1 月 1 日)
    :param lat_range: 纬度范围 (min_lat, max_lat)，默认 None 表示全范围
    :param lon_range: 经度范围 (min_lon, max_lon)，默认 None 表示全范围
    :return: 湿度数据二维数组 (lat, lon) 和裁剪后的经纬度数组
    """
    # 构造文件路径
    file_path = os.path.join(data_dir, f"CDMet_rhu_{year}.nc")  # 文件名为 CDMet_rhu_{year}.nc
    if not os.path.exists(file_path):
        raise FileNotFoundError(f"文件未找到: {file_path}")
    
    # 打开 NetCDF 文件
    with Dataset(file_path, mode="r") as nc_file:
        # 获取湿度变量
        rhu_var = nc_file.variables[variable_name]
    
        # 获取缩放因子和偏移量
        scale_factor = getattr(rhu_var, "scale_factor", 1.0)
        add_offset = getattr(rhu_var, "add_offset", 0.0)
        fill_value = getattr(rhu_var, "_FillValue", -32767.0)
    
        # 获取经纬度
        lat = nc_file.variables["lat"][:]
        lon = nc_file.variables["lon"][:]
    
        # 读取指定天的数据 (二维数组: lat × lon)
        daily_data = rhu_var[day, :, :].astype(np.float32)
    
        # 应用缩放因子和偏移量
        # daily_data = daily_data * scale_factor + add_offset
    
        # 处理无效值 (_FillValue / missing_value)
        daily_data[daily_data == fill_value] = np.nan
    
        # 裁剪经纬度范围（如果指定）
        if lat_range:
            lat_min, lat_max = lat_range
            lat_idx = np.where((lat >= lat_min) & (lat <= lat_max))[0]
            lat = lat[lat_idx]  # 更新纬度数组
            daily_data = daily_data[lat_idx, :]  # 裁剪数据的纬度维度
    
        if lon_range:
            lon_min, lon_max = lon_range
            lon_idx = np.where((lon >= lon_min) & (lon <= lon_max))[0]
            lon = lon[lon_idx]  # 更新经度数组
            daily_data = daily_data[:, lon_idx]  # 裁剪数据的经度维度
    
    return daily_data, lat, lon
    
    
    def plot_rhu_daily_with_inset(data, lat, lon, title, lat_range, lon_range, china_lat_range=None, china_lon_range=None):
    """
    绘制某一天的湿度数据，并在右下角添加中国边界的小窗口。
    :param data: 二维湿度数据 (lat, lon)
    :param lat: 纬度数组
    :param lon: 经度数组
    :param title: 图像标题
    :param lat_range: 主图纬度范围 (min_lat, max_lat)
    :param lon_range: 主图经度范围 (min_lon, max_lon)
    :param china_lat_range: 小窗口纬度范围 (min_lat, max_lat)
    :param china_lon_range: 小窗口经度范围 (min_lon, max_lon)
    """
    # 创建主图
    fig, ax_main = plt.subplots(figsize=(10, 6))
    
    # 绘制湿度数据
    rhu_plot = ax_main.pcolormesh(lon, lat, data, cmap="coolwarm", shading="auto")  # 使用经纬度作为坐标
    cbar = plt.colorbar(rhu_plot, ax=ax_main, fraction=0.05, pad=0.04)
    cbar.set_label("Relative Humidity (%)", fontsize=14)
    
    # 设置标题和坐标轴标签
    ax_main.set_title(title, fontsize=16)
    ax_main.set_xlabel("Longitude (°E)", fontsize=14)
    ax_main.set_ylabel("Latitude (°N)", fontsize=14)
    ax_main.set_xlim(lon_range)
    ax_main.set_ylim(lat_range)
    
    # 设置经纬度刻度字体大小
    ax_main.tick_params(axis="both", which="major", labelsize=13)
    
    # 绘制网格
    ax_main.grid(True, linestyle="--", linewidth=0.5, alpha=0.7)
    
    # 绘制中国边界
    gdf_China.boundary.plot(ax=ax_main, color="black", linewidth=1, label="China Boundary")
    # ax_main.legend()
    
    # 创建小窗口 (inset) 位于主图右下角，与主图边框完全契合
    inset_width = 0.16  # 小窗口宽度占主图的比例
    inset_height = 0.26  # 小窗口高度占主图的比例
    # ax_inset = fig.add_axes([1 - inset_width - 0.05, 0.05, inset_width, inset_height])  # [x, y, width, height]
    ax_inset = fig.add_axes([0.682, 0.11, inset_width, inset_height])  # [x, y, width, height]
    
    # 设置小窗口经纬度范围
    if china_lat_range:
        ax_inset.set_ylim(china_lat_range)
    if china_lon_range:
        ax_inset.set_xlim(china_lon_range)
    
    # 绘制中国边界在小窗口中
    gdf_China.boundary.plot(ax=ax_inset, color="black", linewidth=1)
    
    # 在小窗口中标注主图的经纬度范围
    """
    rect = plt.Rectangle(
        (lon_range[0], lat_range[0]),  # 起点 (左下角)
        lon_range[1] - lon_range[0],  # 矩形宽度
        lat_range[1] - lat_range[0],  # 矩形高度
        linewidth=1.5, edgecolor="red", facecolor="none"
    )
    ax_inset.add_patch(rect)
    
    # 设置小窗口坐标标签
    ax_inset.set_xlabel("Lon (°E)", fontsize=10)
    ax_inset.set_ylabel("Lat (°N)", fontsize=10)
    ax_inset.tick_params(labelsize=8)    
    """
    
    # 隐藏小窗口的经纬度标签，仅保留方框
    ax_inset.tick_params(left=False, bottom=False, labelleft=False, labelbottom=False)
    
    # 显示图像
    plt.show()
    
    # 主程序
    if __name__ == "__main__":
    # 指定年份和天数
    year = 2020
    day = 0  # 0 表示 1 月 1 日
    
    # 指定绘图的经纬度范围 (例如中国区域)
    lat_range = (16, 54)  # 主图纬度范围
    lon_range = (73, 137)  # 主图经度范围
    
    # 指定小窗口的经纬度范围 (覆盖整个中国)
    china_lat_range = (2, 25)  # 小窗口纬度范围
    china_lon_range = (102, 125)  # 小窗口经度范围
    
    # 读取 2020 年 1 月 1 日的湿度数据
    daily_data, lat, lon = read_daily_rhu(data_dir, year, day, lat_range, lon_range)
    
    # 打印数据信息
    print(f"2020 年 1 月 1 日数据形状: {daily_data.shape}")
    print(f"纬度范围: {lat.min()} 至 {lat.max()}")
    print(f"经度范围: {lon.min()} 至 {lon.max()}")
    print(f"数据范围: 最小值={np.nanmin(daily_data)}, 最大值={np.nanmax(daily_data)}")  # 排除 NaN
    
    # 绘制湿度分布图，带小窗口
    plot_rhu_daily_with_inset(
        daily_data, lat, lon,
        title="RHU in Jan.1 2020",
        lat_range=lat_range,
        lon_range=lon_range,
        china_lat_range=china_lat_range,
        china_lon_range=china_lon_range
    )

    import os
    import numpy as np
    from netCDF4 import Dataset
    import rasterio
    from rasterio.transform import from_origin
    
    
    def read_nc_file(file_path):
    """
    读取 .nc 文件，提取变量信息和经纬度数据。
    """
    if not os.path.exists(file_path):
        raise FileNotFoundError(f"文件未找到: {file_path}")
    
    # 打开 NetCDF 文件
    nc_file = Dataset(file_path, mode="r")
    return nc_file
    
    
    def extract_data_in_bbox(nc_file, variable_name, bbox, day):
    """
    提取指定经纬度范围的数据，并应用 scale_factor 和 add_offset。
    同时将缺测值 (_FillValue) 替换为 np.nan。
    """
    # 获取变量和经纬度
    rhu_var = nc_file.variables[variable_name]
    lat = nc_file.variables["lat"][:]
    lon = nc_file.variables["lon"][:]
    
    # 获取变量属性
    """
    scale_factor = getattr(rhu_var, "scale_factor", 1.0)
    print(f"缩放因子scale_factor为{scale_factor}")
    add_offset = getattr(rhu_var, "add_offset", 0.0)
    print(f"偏移值add_offset为{add_offset}")    
    """
    
    fill_value = getattr(rhu_var, "_FillValue", np.nan)
    # print(f"缺失值fill_value为{fill_value}")
    
    # 经纬度范围索引
    lat_min, lat_max, lon_min, lon_max = bbox
    lat_indices = np.where((lat >= lat_min) & (lat <= lat_max))[0]
    lon_indices = np.where((lon >= lon_min) & (lon <= lon_max))[0]
    
    if len(lat_indices) == 0 or len(lon_indices) == 0:
        raise ValueError("指定的经纬度范围无数据")
    
    # 提取当天数据
    daily_data = rhu_var[day, :, :].astype(np.float32)
    
    # 替换缺测值为 np.nan
    daily_data[daily_data == fill_value] = np.nan
    
    # 应用 scale_factor 和 add_offset
    # daily_data = daily_data * scale_factor + add_offset
    
    # 提取指定经纬度范围的数据子集
    subset_data = daily_data[np.min(lat_indices):np.max(lat_indices) + 1,
                             np.min(lon_indices):np.max(lon_indices) + 1]
    subset_lat = lat[np.min(lat_indices):np.max(lat_indices) + 1]
    subset_lon = lon[np.min(lon_indices):np.max(lon_indices) + 1]
    
    return subset_data, subset_lat, subset_lon
    
    
    def save_as_tif(output_dir, year, day, data, lat, lon, file_prefix="CDMet_rhu_GBA"):
    """
    将数据保存为 .tif 文件，并设置 NoData 值为 np.nan。
    """
    # 创建年份目录
    year_dir = os.path.join(output_dir, f"RHU_GBA_{year}")
    if not os.path.exists(year_dir):
        os.makedirs(year_dir)
    
    # 计算栅格分辨率和仿射变换参数
    lon_res = lon[1] - lon[0]  # 经度分辨率
    lat_res = lat[1] - lat[0]  # 纬度分辨率
    
    # 仿射变换 (左上角原点)
    transform = from_origin(west=lon[0], north=lat[-1], xsize=lon_res, ysize=-lat_res)
    
    # 设置保存路径
    output_file = os.path.join(year_dir, f"{file_prefix}_{year}_{day + 1:03d}.tif")
    
    # 保存为 GeoTIFF 文件
    with rasterio.open(
        output_file,
        "w",
        driver="GTiff",
        height=data.shape[0],
        width=data.shape[1],
        count=1,
        dtype=data.dtype,
        crs="EPSG:4326",
        transform=transform,
        nodata=np.nan  # 设置 NoData 值为 np.nan
    ) as dst:
        dst.write(data, 1)
    
    print(f"文件已保存: {output_file}")
    
    
    def process_nc_to_tif(file_path, output_dir, bbox, year, variable_name="rhu"):
    """
    主函数：读取 .nc 文件，提取指定经纬度范围的数据并保存为 .tif 文件。
    """
    # 打开 NetCDF 文件（在主函数中打开，避免上下文管理器关闭文件）
    nc_file = read_nc_file(file_path)
    
    # 获取时间维度大小
    num_days = nc_file.variables[variable_name].shape[0]
    print(f"检测到 {year} 年有 {num_days} 天的数据")
    
    # 遍历每一天的数据
    for day in range(num_days):
        print(f"正在处理 {year} 年第 {day + 1} 天的数据...")
        # 提取指定范围的数据
        subset_data, subset_lat, subset_lon = extract_data_in_bbox(
            nc_file, variable_name, bbox, day
        )
    
        # 保存为 .tif 文件
        save_as_tif(output_dir, year, day, subset_data, subset_lat, subset_lon)
    
    # 关闭 NetCDF 文件
    nc_file.close()
    print(f"{year} 年的数据处理完成！")
    
    
    def batch_process_nc_files(input_dir, output_dir, bbox, start_year, end_year, variable_name="rhu"):
    """
    批量处理指定年份范围内的 .nc 文件。
    """
    for year in range(start_year, end_year + 1):
        print(f"开始处理 {year} 年的数据...")
        # 构造 .nc 文件路径
        file_path = os.path.join(input_dir, f"CDMet_rhu_{year}.nc")
        if not os.path.exists(file_path):
            print(f"警告: {file_path} 不存在，跳过该年份！")
            continue
    
        # 处理当前年份的 .nc 文件
        process_nc_to_tif(file_path, output_dir, bbox, year, variable_name)
    
    print("所有年份的数据处理完成！")
    
    
    # 批量处理调用示例
    if __name__ == "__main__":
    # 输入目录和输出目录
    input_directory = r"C:\Database\3_CDMet\Relative humidity"
    output_directory = r"C:\Database\3_CDMet\RHU_GBA"
    bounding_box = (22, 24, 111, 115)  # 纬度范围 (22-24)，经度范围 (111-115)
    
    # 年份范围
    start_year = 2000
    end_year = 2020
    
    # 批量处理
    batch_process_nc_files(input_directory, output_directory, bounding_box, start_year, end_year)