Python如何处理卫星云图？netCDF4库指南

import netCDF4
import numpy as np

# 假设你有一个名为 'satellite_data.nc' 的NetCDF文件
# 实际工作中，文件路径可能更复杂，比如从某个数据中心下载的
file_path = 'satellite_data.nc'

try:
    # 打开NetCDF文件，'r' 表示只读模式
    with netCDF4.Dataset(file_path, 'r') as ds:
        print(f"成功打开文件: {file_path}\n")

        # 1. 查看全局属性 (Global Attributes)
        # 这些通常包含了关于整个数据集的元信息，比如数据来源、创建者、版本等
        print("--- 全局属性 (Global Attributes) ---")
        for attr_name in ds.ncattrs():
            print(f"  {attr_name}: {getattr(ds, attr_name)}")
        print("\n")

        # 2. 查看维度 (Dimensions)
        # 维度定义了数据的形状，比如时间、纬度、经度、高度等
        print("--- 维度 (Dimensions) ---")
        for dim_name, dim_obj in ds.dimensions.items():
            # dim_obj.isunlimited 表示这个维度是否是无限增长的，通常是时间维度
            print(f"  {dim_name}: size={len(dim_obj)}, is_unlimited={dim_obj.isunlimited()}")
        print("\n")

        # 3. 查看变量 (Variables)
        # 这是最核心的部分，包含了实际的数据和每个变量的属性
        print("--- 变量 (Variables) ---")
        for var_name, var_obj in ds.variables.items():
            print(f"  变量名: {var_name}")
            print(f"    数据类型: {var_obj.dtype}")
            print(f"    维度: {var_obj.dimensions}")
            print(f"    形状: {var_obj.shape}")

            # 遍历变量的属性
            print("    属性:")
            for attr in var_obj.ncattrs():
                print(f"      {attr}: {getattr(var_obj, attr)}")

            # 尝试读取一小部分数据或查看统计信息，避免加载大文件到内存
            if var_obj.ndim > 0: # 确保变量有数据维度
                # 对于非常大的变量，只看前几个或统计信息
                if np.prod(var_obj.shape) > 10000: # 假设超过1万个元素就算大
                    print(f"    数据预览 (前5个): {var_obj[:5]}")
                    print(f"    数据统计 (min/max): {var_obj[:].min()} / {var_obj[:].max()}")
                else:
                    print(f"    数据预览: {var_obj[:]}")
                    print(f"    数据统计 (min/max): {var_obj[:].min()} / {var_obj[:].max()}")
            else:
                print(f"    数据值: {var_obj[()]}") # 标量变量
            print("-" * 30)

except FileNotFoundError:
    print(f"错误：文件 '{file_path}' 未找到。请确保文件路径正确。")
except Exception as e:
    print(f"处理文件时发生错误: {e}")

import netCDF4
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np

# 假设文件路径和变量名
file_path = 'satellite_data.nc' # 替换为你的NetCDF文件路径
data_var_name = 'ir_band'      # 卫星云图数据变量名，例如红外波段
lat_var_name = 'lat'           # 纬度变量名
lon_var_name = 'lon'           # 经度变量名

try:
    with netCDF4.Dataset(file_path, 'r') as ds:
        # 1. 提取数据和坐标
        # 注意：这里假设ir_band是 (time, lat, lon) 结构，我们取第一个时间步
        # 实际数据结构可能不同，需要根据ds.variables[data_var_name].dimensions调整

        # 尝试获取经纬度数据
        try:
            lons = ds.variables[lon_var_name][:]
            lats = ds.variables[lat_var_name][:]
        except KeyError:
            print(f"警告: 未找到经纬度变量 '{lon_var_name}' 或 '{lat_var_name}'。尝试从数据维度推断。")
            # 如果经纬度是维度，可能直接就是维度变量
            if lon_var_name in ds.dimensions and lat_var_name in ds.dimensions:
                lons = ds.variables[lon_var_name][:]
                lats = ds.variables[lat_var_name][:]
            else:
                raise ValueError("无法找到有效的经纬度信息。")

        # 提取卫星云图数据
        # 假设数据是三维的 (time, lat, lon)，我们取第一个时间步
        # 如果是二维 (lat, lon)，直接取即可
        if ds.variables[data_var_name].ndim == 3:
            # 假设时间维度是第一个维度
            image_data = ds.variables[data_var_name][0, :, :] 
        elif ds.variables[data_var_name].ndim == 2:
            image_data = ds.variables[data_var_name][:, :]
        else:
            raise ValueError(f"数据变量 '{data_var_name}' 维度不符合预期，当前维度: {ds.variables[data_var_name].ndim}")

        # 处理可能存在的填充值（_FillValue）或无效值
        if '_FillValue' in ds.variables[data_var_name].ncattrs():
            fill_value = ds.variables[data_var_name]._FillValue
            image_data = np.ma.masked_equal(image_data, fill_value)

        # 2. 创建绘图
        fig = plt.figure(figsize=(10, 8))

        # 定义地图投影。对于全球或大区域，PlateCarree (等经纬度) 常用且简单
        # 对于特定卫星（如地球静止轨道卫星），可能需要更特殊的投影，如 Geostationary
        ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

        # 绘制图像数据
        # pcolormesh 适用于不规则网格或需要精细控制颜色映射的情况
        # imshow 适用于规则网格，且性能通常更好，但可能需要调整extent参数
        # 这里的 extent=[lons.min(), lons.max(), lats.min(), lats.max()] 是一个简化
        # 更准确的做法是使用 transform=ccrs.PlateCarree()

        # 自动调整颜色范围，或者根据数据特性手动设定 vmin, vmax
        # 例如，对于红外数据，可能希望反转颜色，让冷（高云）显示为白色或亮色
        # cmap = 'gray_r' 或 'viridis'

        # 确保数据是浮点类型，避免整数截断导致颜色不连续
        plot = ax.pcolormesh(lons, lats, image_data.astype(float), 
                             transform=ccrs.PlateCarree(), 
                             cmap='viridis', # 或者 'gray_r' 对于红外云图
                             shading='auto') # 'auto' 尝试智能选择 flat 或 gouraud

        # 添加地理要素
        ax.coastlines(resolution='50m', color='black', linewidth=0.8) # 分辨率可选 '10m', '50m', '110m'
        ax.add_feature(ccrs.feature.BORDERS, linestyle=':', edgecolor='gray')
        ax.add_feature(ccrs.feature.LAKES, alpha=0.5, facecolor='blue')
        ax.add_feature(ccrs.feature.RIVERS, edgecolor='blue')

        # 添加网格线和标签
        gl = ax.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False)
        gl.top_labels = False # 不显示顶部标签
        gl.right_labels = False # 不显示右侧标签

        # 添加颜色条
        cbar = fig.colorbar(plot, ax=ax, orientation='vertical', pad=0.05, shrink=0.7)
        cbar.set_label(f'{data_var_name} ({ds.variables[data_var_name].units if "units" in ds.variables[data_var_name].ncattrs() else "无单位"})')

        # 设置标题
        ax.set_title(f'卫星云图: {data_var_name} (第一时间步)', fontsize=14)

        plt.show()

except FileNotFoundError:
    print(f"错误：文件 '{file_path}' 未找到。")
except ValueError as ve:
    print(f"数据处理错误: {ve}")
except Exception as e:
    print(f"绘图时发生错误: {e}")