代码实战 | 如何使用地形数据去除ERA5低层虚假数据

用户11172986

发布于 2025-01-09 18:39:28

5050

引言

在气象学和气候研究中，理解不同高度层的温度变化对于准确预测天气模式和气候变化至关重要。特别是850 hPa这一高度层，它位于对流层下部，是监测天气系统发展的重要层次之一。然而，在某些地区，尤其是地形复杂的区域，如山脉或高原附近，ERA5再分析数据可能会因为模型分辨率限制而产生虚假的850 hPa数据——这些地方实际上可能是地形表面而非自由大气。

为了更精确地反映真实情况，并排除地形影响导致的数据偏差，我们决定绘制经过地形过滤后的850 hPa温度分布图。这不仅有助于提高数据分析的准确性，还能为气象预报提供更加可靠的支持。通过将ERA5全球再分析数据与中国高分辨率DEM（数字高程模型）相结合，我们可以创建出既科学又直观的地图，帮助研究人员更好地理解和解释气象现象。

方法概述

为了更好地组织代码并提高可读性，我们将整个流程拆分为多个模块。每个模块负责一个具体的功能或步骤，最后提供完整的代码清单。

模块1：加载和准备数据


import xarray as xr

def load_data(era5_path, dem_path):
    """加载ERA5和DEM数据"""
    # 加载ERA5数据
    era5_data = xr.open_dataset(era5_path)
    
    # 加载中国DEM数据并调整坐标名称
    dem_data = xr.open_dataset(dem_path)
    if'LAT'in dem_data.coords and'LON'in dem_data.coords:
        dem_data = dem_data.rename({'LAT': 'latitude', 'LON': 'longitude'})
    
    return era5_data, dem_data['elevation']

模块2：插值DEM到ERA5网格上

def interpolate_dem_to_era5(dem_elevation, era5_data):
    """将DEM数据插值到ERA5网格上"""
    dem_interpolated = dem_elevation.interp_like(era5_data)
    return dem_interpolated

这里使用了xarray的内部函数interp_like，这方便我们将较细的地形数据插值到era5的格点上，方便数据筛选

模块3：应用地形掩码

def apply_terrain_mask(t, z, dem_interpolated):
    """根据地形高度应用掩码"""
    condition = (z/9.8 < dem_interpolated)
    t_masked = t.where(~condition | dem_interpolated.isnull(), np.nan)
    return t_masked

这里我们使用了位势高度当作实际高度来筛选数据

模块4：绘制850 hPa温度分布图


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

def plot_temperature_distribution(t_masked, cmap='coolwarm'):
    """绘制850 hPa温度分布图"""
    temperature_850 = t_masked.sel(level=850, method='nearest').isel(time=[0])

    plt.figure(figsize=(12, 8))
    ax = plt.axes(projection=ccrs.PlateCarree())

    contourf = temperature_850[0].plot.contourf(
        ax=ax,
        transform=ccrs.PlateCarree(),
        levels=np.arange(270, 300, 2),
        cmap=cmap,
        extend='both'
    )

    ax.coastlines()
    plt.colorbar(contourf, orientation='horizontal', pad=0.07, label='Temperature (K)', shrink=0.8)

    plt.title('Filtered 850 hPa Temperature Distribution Over China on Sep 4, 2022', fontsize=16)
    plt.xlabel('Longitude', fontsize=14)
    plt.ylabel('Latitude', fontsize=14)

    # 设置字体大小
    for item in ([ax.title, ax.xaxis.label, ax.yaxis.label] +
                 ax.get_xticklabels() + ax.get_yticklabels()):
        item.set_fontsize(12)

    plt.tight_layout()
    plt.show()

完整代码

以下是所有模块组合在一起的完整代码：

import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np

# 模块1：加载和准备数据
def load_data(era5_path, dem_path):
    """加载ERA5和DEM数据"""
    # 加载ERA5数据
    era5_data = xr.open_dataset(era5_path)
    
    # 加载中国DEM数据并调整坐标名称
    dem_data = xr.open_dataset(dem_path)
    if'LAT'in dem_data.coords and'LON'in dem_data.coords:
        dem_data = dem_data.rename({'LAT': 'latitude', 'LON': 'longitude'})
    
    return era5_data, dem_data['elevation']

# 模块2：插值DEM到ERA5网格上
def interpolate_dem_to_era5(dem_elevation, era5_data):
    """将DEM数据插值到ERA5网格上"""
    dem_interpolated = dem_elevation.interp_like(era5_data)
    return dem_interpolated

# 模块3：应用地形掩码
def apply_terrain_mask(t, z, dem_interpolated):
    """根据地形高度应用掩码"""
    condition = (z/9.8 < dem_interpolated)
    t_masked = t.where(~condition | dem_interpolated.isnull(), np.nan)
    return t_masked

# 模块4：绘制850 hPa温度分布图
def plot_temperature_distribution(t_masked, cmap='coolwarm'):
    """绘制850 hPa温度分布图"""
    temperature_850 = t_masked.sel(level=850, method='nearest').isel(time=[0])

    plt.figure(figsize=(12, 8),dpi=200)
    ax = plt.axes(projection=ccrs.PlateCarree())

    contourf = temperature_850[0].plot.contourf(
        ax=ax,
        transform=ccrs.PlateCarree(),
        levels=np.arange(270, 300, 2),
        cmap=cmap,
        extend='both'
    )

    ax.coastlines()

    plt.title('Filtered 850 hPa Temperature Distribution Over China 202308', fontsize=16)
    plt.xlabel('Longitude', fontsize=14)
    plt.ylabel('Latitude', fontsize=14)

    # 设置字体大小
    for item in ([ax.title, ax.xaxis.label, ax.yaxis.label] +
                 ax.get_xticklabels() + ax.get_yticklabels()):
        item.set_fontsize(12)

    plt.tight_layout()
    plt.show()

# 主程序
if __name__ == "__main__":
    era5_filename = "/home/mw/input/era58091/ERA5-2023-08_pl.nc"
    dem_filename = "/home/mw/input/china_dem3276/cldasgrid_dem.nc"

    # 加载数据
    era5_data, dem_elevation = load_data(era5_filename, dem_filename)

    # 插值DEM到ERA5网格上
    dem_interpolated = interpolate_dem_to_era5(dem_elevation, era5_data)

    # 提取温度和位势高度数据，并应用掩码
    t = era5_data.t
    z = era5_data.z
    t_masked = apply_terrain_mask(t, z, dem_interpolated)

    # 绘制850 hPa温度分布图
    plot_temperature_distribution(t_masked)

/opt/conda/lib/python3.9/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/110m_physical/ne_110m_coastline.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

打点

下面用sccater做了下打点，但是效果不佳

有更好的办法也可评论区提出

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
import numpy.ma as ma
import xarray as xr

def plot_temperature_with_nan(t_masked, cmap='coolwarm'):
    """绘制850 hPa温度分布图并在NaN值位置打点"""
    # 选择850 hPa水平的温度数据并检查是否为NaN
    temperature_850 = t_masked.sel(level=850, method='nearest').isel(time=0)
    
    # 创建图形和坐标轴
    fig, ax = plt.subplots(figsize=(12, 8), dpi=200, subplot_kw={'projection': ccrs.PlateCarree()})

    # 绘制温度等值线图
    temperature_850.plot.contourf(
        ax=ax,
        levels=np.arange(270, 300, 2),
        cmap=cmap,
        extend='both',
        transform=ccrs.PlateCarree()
    )

    # 添加海岸线
    ax.coastlines()

    # 找到NaN值的位置，并转换成适合matplotlib scatter的格式
    nan_mask = np.isnan(temperature_850)

    lon, lat = np.meshgrid(temperature_850['longitude'], temperature_850['latitude'])
    lon_nan= lon[nan_mask]
    lat_nan= lat[nan_mask]
    if len(lon_nan) > 0and len(lat_nan) > 0:
        # 使用matplotlib直接在轴上打点显示NaN值的位置
        step=5
        ax.scatter(
            lon_nan[::step], lat_nan[::step],
            s=5, c='black', marker='.', alpha=0.5,
            transform=ccrs.PlateCarree(),
            zorder=1# 提高绘图层级确保点在最上层
        )

    # 设置标题和标签
    ax.set_title('Filtered 850 hPa Temperature Distribution Over China on Sep 4, 2022', fontsize=16)
    ax.set_xlabel('Longitude', fontsize=14)
    ax.set_ylabel('Latitude', fontsize=14)

    # 调整布局以防止标签被裁剪
    plt.tight_layout()

    # 显示图形
    plt.show()

# 主程序
if __name__ == "__main__":
    era5_filename = "/home/mw/input/era58091/ERA5-2023-08_pl.nc"
    dem_filename = "/home/mw/input/china_dem3276/cldasgrid_dem.nc"

    # 加载数据
    era5_data, dem_elevation = load_data(era5_filename, dem_filename)

    # 插值DEM到ERA5网格上
    dem_interpolated = interpolate_dem_to_era5(dem_elevation, era5_data)

    # 提取温度和位势高度数据，并应用掩码
    t = era5_data.t
    z = era5_data.z
    t_masked = apply_terrain_mask(t, z, dem_interpolated)

    # 绘制850 hPa温度分布图并打点显示被掩码区域
    plot_temperature_with_nan(t_masked)