用Python画的ECMF10日四要素垂直剖面图

betasy

用Python画的ECMF10日四要素垂直剖面图部分实现代码：
def plot_profiles(rhu, tem, vvp, u, v, rhu_levels=None, save_path=None):
    m_rhu = rhu['matrix']
    m_tem = tem['matrix']
    m_vvp = vvp['matrix']
    m_u = u['matrix']
    m_v = v['matrix']

    start_date = datetime.strptime(rhu['start_time'], "%Y-%m-%d %H:%M:%S")
    day_times = [0, 24, 48, 72, 120, 168, 216]
    xlocators = [i for i, t in enumerate(rhu['forecast_times']) if t in day_times]
    ylocators = [i for i in list(range(len(rhu['levels'])))]
    f_dates = [start_date + timedelta(hours=i) for i in day_times]
    fig, ax = plt.subplots(figsize=(12, 6), dpi=180, facecolor='white')

    ## RHU plot
if rhu_levels is None:
        rhu_levels = [0, 50, 60, 70, 80, 90, 100]
    rhu_colors = ["#ffffff", "#C8FFBE", "#85CC7F", "#64B35F", "#43993F", "#218020", "#006600"]
    cf = ax.contourf(gaussian_filter(m_rhu, sigma=0.9),
                     levels=rhu_levels, colors=rhu_colors, extend='max')
    c = ax.contour(gaussian_filter(m_rhu, sigma=0.9),
                   levels=rhu_levels, colors="black", extend='max', linewidths=0.4)
    ax.clabel(c, inline=True, fontsize=6, fmt="%1.0f", inline_spacing=10)
    colorbar_ax = fig.add_axes([0.904, 0.24, 0.006, 0.5])
    fig.colorbar(cf, cax=colorbar_ax, ax=cf, shrink=0.6, fraction=0.1, drawedges=True, orientation='vertical')

    ## VVP plot
c2 = ax.contour(gaussian_filter(m_vvp, sigma=0.9),
                   levels=[float(i) / 100 for i in list(range(-100, 100, 10))], colors="blue", extend='max',
                   linewidths=0.4)
    ax.clabel(c2, inline=True, fontsize=6, fmt="%1.2f", inline_spacing=10)

    ## TEM plot
c3 = ax.contour(gaussian_filter(m_tem - 273.15, sigma=0.9),
                   levels=list(range(-80, 37, 2)), colors="red", extend='max', linewidths=0.4)
    ax.clabel(c3, inline=True, fontsize=6, fmt="%1.0f", inline_spacing=10)
    ## WINUV plot
sizes = dict(emptybarb=0.06, spacing=0.1, height=0.6)
    b = ax.barbs(m_u, m_v, sizes=sizes, length=4,
                 linewidth=0.5, pivot="middle")
    ax.xaxis.set_major_locator(FixedLocator(xlocators))
    ax.xaxis.set_minor_locator(LinearLocator(numticks=len(rhu['forecast_times'])))
    ax.set_xticklabels([f.strftime("%m-%d") for f in f_dates], fontsize=8)
    ax.yaxis.set_major_locator(FixedLocator(ylocators))
    ax.set_yticklabels(rhu["levels"][::-1], fontsize=8)
    ax.set_ylabel("Level(hpa)")
    ax.set_title("ECMF\nRHU(shade,%), TEM(red line, $℃$), VVP(black line, Pa/s), WIND(wind vector)",
                horizontalalignment="left", loc ='left')
    if save_path:
        plt.savefig(save_path, bbox_inches='tight', pad_inches=1)
    else:
        plt.show()


欢迎同好交流

betasy

csv文件是站点信息和示例数据，下载后放到一个目录下，然后在profile_plot_eg.py文件中修改data_dir为这个目录。
requirements:
pandas, matplotlib, datetime, scipy.ndimage（等值线平滑）

还有别的绘图问题也可以交流哈
TQECProfiles.py
"""
NAFP_FOR_FTM_HIGH_EC_ANEA
垂直速度  VVP
相对湿度  RHU
气温 TEM
风 WIV WIU
56187  成都市    103.870000 30.750000
56196  绵阳市    104.730000 31.450000
56198  德阳市    104.500000 31.320000
56287  雅安市    103.000000 29.980000
56298  资阳市    104.600000 30.130000
56386  乐山市    103.750000 29.570000
56391  眉山市    103.820000 30.080000
57405  遂宁市    105.550000 30.500000
"""
import json
# import neccessary modules
from datetime import datetime, timedelta

# plotting modules configure
import matplotlib.pyplot as plt
from matplotlib.ticker import FixedLocator, LinearLocator
from scipy.ndimage.filters import gaussian_filter



def site_data_extract(sitedf, siteID, feature):
    onesite = sitedf[(sitedf['siteID'] == siteID) & (sitedf['feature'] == feature)]
    dt = onesite['dt'].drop_duplicates().tolist()[0]
    onesite = onesite.set_index(["times", "level"])["val"].unstack()
    times = onesite.index.tolist()
    levels = onesite.columns.tolist()
    mat = onesite.values.transpose()
    mat = mat[::-1]
    return {"matrix": mat,
            "start_time": dt,
            "forecast_times": [int(t) for t in times],
            "levels": [int(l) for l in levels],
            "site_id": siteID,
            "feature": feature}

def plot_profiles(rhu, tem, vvp, u, v, rhu_levels=None, save_path=None):
    m_rhu = rhu['matrix']
    m_tem = tem['matrix']
    m_vvp = vvp['matrix']
    m_u = u['matrix']
    m_v = v['matrix']

    start_date = datetime.strptime(rhu['start_time'], "%Y-%m-%d %H:%M:%S")
    day_times = [0, 24, 48, 72, 120, 168, 216]
    xlocators = [i for i, t in enumerate(rhu['forecast_times']) if t in day_times]
    ylocators = [i for i in list(range(len(rhu['levels'])))]
    f_dates = [start_date + timedelta(hours=i) for i in day_times]
    fig, ax = plt.subplots(figsize=(12, 6), dpi=180, facecolor='white')

    ## RHU plot
    if rhu_levels is None:
        rhu_levels = [0, 50, 60, 70, 80, 90, 100]
    rhu_colors = ["#ffffff", "#C8FFBE", "#85CC7F", "#64B35F", "#43993F", "#218020", "#006600"]
    cf = ax.contourf(gaussian_filter(m_rhu, sigma=0.9),
                     levels=rhu_levels, colors=rhu_colors, extend='max')
    c = ax.contour(gaussian_filter(m_rhu, sigma=0.9),
                   levels=rhu_levels, colors="black", extend='max', linewidths=0.4)
    ax.clabel(c, inline=True, fontsize=6, fmt="%1.0f", inline_spacing=10)
    colorbar_ax = fig.add_axes([0.904, 0.24, 0.006, 0.5])
    fig.colorbar(cf, cax=colorbar_ax, ax=cf, shrink=0.6, fraction=0.1, drawedges=True, orientation='vertical')

    ## VVP plot
    c2 = ax.contour(gaussian_filter(m_vvp, sigma=0.9),
                   levels=[float(i) / 100 for i in list(range(-100, 100, 10))], colors="blue", extend='max',
                   linewidths=0.4)
    ax.clabel(c2, inline=True, fontsize=6, fmt="%1.2f", inline_spacing=10)

    ## TEM plot
    c3 = ax.contour(gaussian_filter(m_tem - 273.15, sigma=0.9),
                   levels=list(range(-80, 37, 2)), colors="red", extend='max', linewidths=0.4)
    ax.clabel(c3, inline=True, fontsize=6, fmt="%1.0f", inline_spacing=10)
    ## WINUV plot
    sizes = dict(emptybarb=0.06, spacing=0.1, height=0.6)
    b = ax.barbs(m_u, m_v, sizes=sizes, length=4,
                 linewidth=0.5, pivot="middle")
    ax.xaxis.set_major_locator(FixedLocator(xlocators))
    ax.xaxis.set_minor_locator(LinearLocator(numticks=len(rhu['forecast_times'])))
    ax.set_xticklabels([f.strftime("%m-%d") for f in f_dates], fontsize=8)
    ax.yaxis.set_major_locator(FixedLocator(ylocators))
    ax.set_yticklabels(rhu["levels"][::-1], fontsize=8)
    ax.set_ylabel("Level(hpa)")
    ax.set_title("ECMF\nRHU(shade,%), TEM(red line, $℃$), VVP(black line, Pa/s), WIND(wind vector)",
                horizontalalignment="left", loc ='left')
    if save_path:
        plt.savefig(save_path, bbox_inches='tight', pad_inches=1)
    else:
        plt.show()

profile_plot_eg.pyimport pandas as pd
import os
from datetime import datetime, timedelta

from TQECProfiles import site_data_extract, plot_profiles

dt = "20220109200000"

# 加载数据
data_dir = r"E:\Projects\pythonProject\readearth-plots\Data\tqecprofiles"
sites = pd.read_csv(os.path.join(data_dir, "sites.csv"), dtype={"siteID": str})
rhu = pd.read_csv(os.path.join(data_dir, "rhu.csv"), dtype={"siteID": str})
vvp = pd.read_csv(os.path.join(data_dir, "vvp.csv"), dtype={"siteID": str})
tem = pd.read_csv(os.path.join(data_dir, "tem.csv"), dtype={"siteID": str})
u = pd.read_csv(os.path.join(data_dir, "u.csv"), dtype={"siteID": str})
v = pd.read_csv(os.path.join(data_dir, "v.csv"), dtype={"siteID": str})

# 将数据和站点信息合并
# rhu_df = rhu.merge(sites, how="left", on=["lng", "lat"])
# vvp_df = vvp.merge(sites, how="left", on=["lng", "lat"])
# tem_df = tem.merge(sites, how="left", on=["lng", "lat"])
# udf = u.merge(sites, how="left", on=["lng", "lat"])
# vdf = v.merge(sites, how="left", on=["lng", "lat"])

siteIDs = sites["siteID"].drop_duplicates().tolist()
for site in siteIDs:
    dt_plus8 = (datetime.strptime(dt, "%Y%m%d%H%M%S") + timedelta(8)).strftime("%Y%m%d%H%M%S")
    output_path = f"EC_PROFILE_{dt_plus8}_{site}.png"
    if os.path.exists(output_path):
        print(os.path.split(output_path)[1] + "already exists...")
    else:
        r_data = site_data_extract(rhu, site, "RHU")
        vvp_data = site_data_extract(vvp, site, "VVP")
        t_data = site_data_extract(tem, site, "TEM")
        u_data = site_data_extract(u, site, "WIU")
        v_data = site_data_extract(v, site, "WIV")
        plot_profiles(r_data, t_data, vvp_data, u_data, v_data, rhu_levels=None, save_path=output_path)
        print(os.path.split(output_path)[1] + " plotted...")

木人

有完整数据代码吗?

betasy

1. def site_data_extract(sitedf, siteID, feature):
   
2.     onesite = sitedf[(sitedf['siteID'] == siteID) & (sitedf['feature'] == feature)]
   
3.     dt = onesite['dt'].drop_duplicates().tolist()[0]
   
4.     onesite = onesite.set_index(["times", "level"])["val"].unstack()
   
5.     times = onesite.index.tolist()
   
6.     levels = onesite.columns.tolist()
   
7.     mat = onesite.values.transpose()
   
8.     mat = mat[::-1]
   
9.     return {"matrix": mat,
   
10.             "start_time": dt,
    
11.             "forecast_times": [int(t) for t in times],
    
12.             "levels": [int(l) for l in levels],
    
13.             "site_id": siteID,
    
14.             "feature": feature}

def site_data_extract(sitedf, siteID, feature):<br />     onesite = sitedf[(sitedf['siteID'] == siteID) & (sitedf['feature'] == feature)]<br />     dt = onesite['dt'].drop_duplicates().tolist()[0]<br />     onesite = onesite.set_index(["times", "level"])["val"].unstack()<br />     times = onesite.index.tolist()<br />     levels = onesite.columns.tolist()<br />     mat = onesite.values.transpose()<br />     mat = mat[::-1]<br />     return {"matrix": mat,<br />             "start_time": dt,<br />             "forecast_times": [int(t) for t in times],<br />             "levels": [int(l) for l in levels],<br />             "site_id": siteID,<br />             "feature": feature}

数据是从天擎接口上下来的json数组，就是附件中的这种形式，然后用上面这个函数把站点数据提取出来

yqcxm

图很漂亮，收藏先

jonny_521

可否提供完整的示例？

独孤酒见

楼主高人啊。这玩的溜得很。
我还在学，学出一身汗。
谢谢楼主分享。

mico2008

楼主威武，赞

Crystal_qiong

厉害啊就，要好好学python

jonny_521

感谢感谢，好好学一学！