请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 1499|回复: 14

[经验总结] 用Python画的ECMF10日四要素垂直剖面图

[复制链接]

新浪微博达人勋

发表于 2022-1-7 13:27:59 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 betasy 于 2022-1-7 13:27 编辑

用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()


欢迎同好交流

四要素垂直剖面图

四要素垂直剖面图
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2022-1-10 14:54:56 | 显示全部楼层
本帖最后由 betasy 于 2022-1-10 15:00 编辑

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...")



EC_PROFILE_20220117200000_56386.png
EC_PROFILE_20220117200000_56391.png
EC_PROFILE_20220117200000_57405.png
EC_PROFILE_20220117200000_56198.png
EC_PROFILE_20220117200000_56196.png
EC_PROFILE_20220117200000_56187.png
EC_PROFILE_20220117200000_56298.png
EC_PROFILE_20220117200000_56287.png

示例.zip

8.24 MB, 下载次数: 841, 下载积分: 金钱 -5

密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

发表于 2022-1-7 16:47:46 | 显示全部楼层
有完整数据代码吗?
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2022-1-9 20:16:59 | 显示全部楼层
本帖最后由 betasy 于 2022-1-9 20:18 编辑

  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}

数据是从天擎接口上下来的json数组,就是附件中的这种形式,然后用上面这个函数把站点数据提取出来
Snipaste_2022-01-09_20-18-31.png

TEM2.json

153.23 KB, 下载次数: 12, 下载积分: 金钱 -5

密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2022-1-10 11:20:56 | 显示全部楼层
图很漂亮,收藏先
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2022-1-10 12:32:21 | 显示全部楼层
可否提供完整的示例?
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2022-1-10 15:32:31 | 显示全部楼层
楼主高人啊。这玩的溜得很。
我还在学,学出一身汗。
谢谢楼主分享。
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2022-1-11 10:20:45 | 显示全部楼层
楼主威武,赞
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2022-1-11 10:32:38 | 显示全部楼层
厉害啊就,要好好学python
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2022-1-11 16:59:02 | 显示全部楼层
感谢感谢,好好学一学!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表