From 693494a38a8ade6c79ebdcb7edaa7d8e69d8dcde Mon Sep 17 00:00:00 2001 From: jtang <26022410@qq.com> Date: Fri, 17 Mar 2023 15:30:16 +0800 Subject: [PATCH] add era5 support MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 增加era5数据 修改实况数据，增加20-20降水，使得08-08和20-20的降水量的中间时刻为环流场时刻 --- apps/reanalysis_maps/app.py | 18 +- apps/reanalysis_maps/app_era5.py | 414 +++++++++++++++++++ apps/reanalysis_maps/draw_maps.py | 18 +- apps/reanalysis_maps/draw_maps_era5.py | 544 +++++++++++++++++++++++++ 4 files changed, 987 insertions(+), 7 deletions(-) create mode 100644 apps/reanalysis_maps/app_era5.py create mode 100644 apps/reanalysis_maps/draw_maps_era5.py diff --git a/apps/reanalysis_maps/app.py b/apps/reanalysis_maps/app.py index 8c58e6d..910265e 100644 --- a/apps/reanalysis_maps/app.py +++ b/apps/reanalysis_maps/app.py @@ -96,9 +96,20 @@ def main(): # date_obj.strftime('%Y%m%d000000'), data_code="SURF_CHN_MUL_DAY", sta_levels="011,012,013", # elements="Station_Id_C,Lat,Lon,Alti,TEM_Max,TEM_Min,VIS_Min,PRE_Time_0808,WIN_S_Max") # Also support cmadaas server - obs_data = cmadaas_obs_by_time( - date_obj.strftime('%Y%m%d000000'), data_code="SURF_CHN_MUL_DAY", sta_levels="011,012,013", - elements="Station_Id_C,Lat,Lon,Alti,TEM_Max,TEM_Min,VIS_Min,PRE_Time_0808,WIN_S_Max") + # SURF_CHN_MUL_DAY(_N) 接口里 + # 20220405000000 + # PRE_Time_0808 指5日08-6日08 + # PRE_Time_2020，是指4日20-5日20 + if data_time == '00': + # 20 + obs_data = cmadaas_obs_by_time( + date_obj.strftime('%Y%m%d000000'), data_code="SURF_CHN_MUL_DAY", sta_levels="011,012,013", + elements="Station_Id_C,Lat,Lon,Alti,TEM_Max,TEM_Min,VIS_Min,PRE_Time_2020,WIN_S_Max") + else: + # 0520--- 05 get 0508-0608 + obs_data = cmadaas_obs_by_time( + date_obj.strftime('%Y%m%d000000'), data_code="SURF_CHN_MUL_DAY", sta_levels="011,012,013", + elements="Station_Id_C,Lat,Lon,Alti,TEM_Max,TEM_Min,VIS_Min,PRE_Time_0808,WIN_S_Max") if obs_data is not None: state.img3 = draw_maps.draw_observation(obs_data, date_obj, map_region) else: @@ -163,6 +174,7 @@ def read_variable(varname, date_obj): result = requests.get(filepath_html%(varname)) if result.status_code == 200: +# if result.status_code == 200 and yyyy != 2014: data = xr.open_dataset(filepath_local%(varname)) else: data = xr.open_dataset(filepath_remote) diff --git a/apps/reanalysis_maps/app_era5.py b/apps/reanalysis_maps/app_era5.py new file mode 100644 index 0000000..4b5d678 --- /dev/null +++ b/apps/reanalysis_maps/app_era5.py @@ -0,0 +1,414 @@ +# _*_ coding: utf-8 _*_ + +# Copyright (c) 2020 NMC Developers. +# Distributed under the terms of the GPL V3 License. + +""" +检索特定日期和范围的CFSR再分析数据, 并绘制相关的天气图. +""" + +import sys +import requests +import datetime +from multiprocessing import Process, Manager +import numpy as np +import xarray as xr +import streamlit as st +from metpy.units import units + +# set page title +st.set_page_config(page_title="历史天气图分析", layout="wide") + +from nmc_met_io.retrieve_cmadaas import cmadaas_obs_by_time +from nmc_met_io.retrieve_cimiss_server import cimiss_obs_by_time +from nmc_met_graphics.util import get_map_regions +from nmc_met_graphics.web import SessionState, ipyplot + +sys.path.append('.') +import draw_maps_era5 as draw_maps + +# set session state +state = SessionState.get(img1=None, img2=None, img3=None) + + +def main(): + # application title + st.title("历史天气图分析") + st.header('——检索1959年1月1日以来的地面观测日值和天气图') + + # application information + st.sidebar.image('http://image.nmc.cn/assets/img/index/nmc_logo_3.png', width=300) + st.sidebar.markdown('**天气预报技术研发室**') + + # Input data source + data_source = st.sidebar.selectbox('数据源:', ('CFSR', 'ERA5'), index=1) + + # Input data date + data_date = st.sidebar.date_input( + "输入日期:", datetime.date(2020, 7, 1), + min_value=datetime.date(1959, 1, 1), + max_value=datetime.date.today() - datetime.timedelta(days=2)) + + # Input data time + data_time = st.sidebar.selectbox('选择时刻(UTC):', ('00', '06', '12', '18'), index=2) + + # construct datetime string + datetime_str = data_date.strftime('%Y%m%d') + data_time + date_obj = datetime.datetime.strptime(datetime_str,'%Y%m%d%H') + + # subset data + map_regions = get_map_regions() + select_items = list(map_regions.keys()) + select_items.append('自定义') + select_item = st.sidebar.selectbox('选择空间范围:', select_items) + if select_item == '自定义': + map_region_str = st.sidebar.text_input('输入空间范围[West, East, South, North]:', '70, 140, 10, 65') + map_region = list(map(float, map_region_str.split(", "))) + else: + map_region = map_regions[select_item] + + # draw the figures + if st.sidebar.button('绘制天气图'): + # load data + data = load_variables(date_obj, map_region=map_region, data_source=data_source) + + # draw the synoptic compsite + fig = draw_maps.draw_composite_map( + date_obj, data['t850'], data['u200'], data['v200'], data['u500'], + data['v500'], data['mslp'], data['gh500'], data['u850'], data['v850'], data['pwat']) + state.img1 = fig + + # draw weather analysis maps + manager = Manager() + return_dict = manager.dict() + p = Process(target=draw_maps.draw_weather_analysis, args=(date_obj, data, map_region, return_dict, data_source)) + p.start() + p.join() + if return_dict[0] is not None: + state.img2 = return_dict[0] + else: + st.info('制作各层天气图失效!') + state.img2 = None + return_dict.clear() + manager.shutdown() + p.close() + del data + + # load observation data + # obs_data = cimiss_obs_by_time( + # date_obj.strftime('%Y%m%d000000'), data_code="SURF_CHN_MUL_DAY", sta_levels="011,012,013", + # elements="Station_Id_C,Lat,Lon,Alti,TEM_Max,TEM_Min,VIS_Min,PRE_Time_0808,WIN_S_Max") + # Also support cmadaas server + # SURF_CHN_MUL_DAY(_N) 接口里 + # 20220405000000 + # PRE_Time_0808 指5日08-6日08 + # PRE_Time_2020，是指4日20-5日20 + if data_time == '00': + # 20 + obs_data = cmadaas_obs_by_time( + date_obj.strftime('%Y%m%d000000'), data_code="SURF_CHN_MUL_DAY", sta_levels="011,012,013", + elements="Station_Id_C,Lat,Lon,Alti,TEM_Max,TEM_Min,VIS_Min,PRE_Time_2020,WIN_S_Max") + else: + # 0520--- 05 get 0508-0608 + obs_data = cmadaas_obs_by_time( + date_obj.strftime('%Y%m%d000000'), data_code="SURF_CHN_MUL_DAY", sta_levels="011,012,013", + elements="Station_Id_C,Lat,Lon,Alti,TEM_Max,TEM_Min,VIS_Min,PRE_Time_0808,WIN_S_Max") + if obs_data is not None: + state.img3 = draw_maps.draw_observation(obs_data, date_obj, map_region) + else: + state.img3 = None + + if state.img1 is None or state.img2 is None: + st.info('请点击左侧**绘制天气图**按钮生成或更新图像.') + else: + # display observation + if state.img3 is not None: + st.markdown( + ''' + ------ + ### 选择站点实况观测''') + options = [key for key in state.img3.keys()] + option = st.selectbox('', options, 0) + st.plotly_chart(state.img3[option], use_container_width=True) + + # display synoptic composite + st.markdown( + ''' + ------ + ### 环流形势综合图''') + st.pyplot(state.img1) + + # display weather analysis maps + st.markdown( + ''' + ------ + ### 点击天气图弹出放大''') + images = np.asarray([*state.img2.values()], dtype=object) + labels = np.asarray([*state.img2.keys()]) + html = ipyplot.display_image_gallery(images, labels, img_width=200) + st.markdown(html, unsafe_allow_html=True) + # options = [key for key in state.img2.keys()] + #options = st.multiselect( + # '', options, ['500hPa_height', 'precipitable_water', 'mean_sea_level_pressure']) + #for option in options: + # st.image(state.img2[option], use_column_width=True) + + st.sidebar.markdown( + ''' + ------ + [CFSR](https://climatedataguide.ucar.edu/climate-data/climate-forecast-system-reanalysis-cfsr)(CLIMATE FORECAST SYSTEM REANALYSIS) + 再分析资料是NCEP提供的第三代再分析产品. 本程序从Albany大学Thredds服务器检索指定日期时刻及范围的CFSR数据（从1979年1月1日以来每日4次, 全球范围）, + 并绘制高空环流形势天气图. + ''') + + +def read_variable(varname, date_obj, data_source='CFSR'): + """ + read varaible from load or remote thredds server. + """ + # To make parsing the date easier, convert it into a datetime object + # and get it into various formats + yyyy = date_obj.year + yyyymmdd = date_obj.strftime('%Y%m%d') + + filepath_html = "http://10.28.49.118:8080/thredds/dodsC/CFSR/%s/%s.%s.0p5.anl.nc.html"%(yyyy,'%s',yyyy) + filepath_local = "http://10.28.49.118:8080/thredds/dodsC/CFSR/%s/%s.%s.0p5.anl.nc"%(yyyy,'%s',yyyy) + filepath_remote = "https://tangjian%%40cma.gov.cn:051087908986@rda.ucar.edu/thredds/dodsC/files/g/ds094.0/%s/cdas1.%s.pgrbanl.tar"%(yyyy,yyyymmdd) + filepath_era5 = "http://10.28.49.118:8080/thredds/dodsC/ERA5/china/%s/era5.%s.%s.nc"%('%s','%s',yyyymmdd) + + if data_source == 'CFSR': + result = requests.get(filepath_html%(varname)) + if result.status_code == 200: + data = xr.open_dataset(filepath_local%(varname)) + else: + data = xr.open_dataset(filepath_remote) + data = data.sel(time=date_obj) + else: + print(filepath_era5) + data = xr.open_dataset(filepath_era5%(varname,varname)) + data = data.sel(time=date_obj) + + return data + + +def load_variables(date_obj, map_region=[50, 160, 6, 60], data_source='CFSR'): + """ + Load the variables from UAlbany's opendap server + + Args: + date_obj (datetime): a datetime object + """ + if data_source=='CFSR': + if int(date_obj.year) >= 2020: # add rda.ucar.edu data source by tangjian. + st.info('从rda.ucar.edu载入CFSR数据(1x1度, Waiting ~30s.)') + + # construct sub region + sub_region = {'lon':slice(map_region[0], map_region[1]), + 'lat':slice(map_region[3], map_region[2])} + + # Subset and load data + subdata = {} + + # wind field + my_bar = st.progress(0) + data = read_variable('u-component_of_wind_isobaric', date_obj) + subdata['u200'] = data['u-component_of_wind_isobaric'].sel(isobaric1=20000.0, **sub_region).load() ; my_bar.progress(4) + subdata['u500'] = data['u-component_of_wind_isobaric'].sel(isobaric1=50000.0, **sub_region).load() ; my_bar.progress(8) + subdata['u700'] = data['u-component_of_wind_isobaric'].sel(isobaric1=70000.0, **sub_region).load() ; my_bar.progress(12) + subdata['u850'] = data['u-component_of_wind_isobaric'].sel(isobaric1=85000.0, **sub_region).load() ; my_bar.progress(16) + subdata['u925'] = data['u-component_of_wind_isobaric'].sel(isobaric1=92500.0, **sub_region).load() ; my_bar.progress(20) + data = read_variable('v-component_of_wind_isobaric', date_obj) + subdata['v200'] = data['v-component_of_wind_isobaric'].sel(isobaric1=20000.0, **sub_region).load() ; my_bar.progress(24) + subdata['v500'] = data['v-component_of_wind_isobaric'].sel(isobaric1=50000.0, **sub_region).load() ; my_bar.progress(28) + subdata['v700'] = data['v-component_of_wind_isobaric'].sel(isobaric1=70000.0, **sub_region).load() ; my_bar.progress(32) + subdata['v850'] = data['v-component_of_wind_isobaric'].sel(isobaric1=85000.0, **sub_region).load() ; my_bar.progress(36) + subdata['v925'] = data['v-component_of_wind_isobaric'].sel(isobaric1=92500.0, **sub_region).load() ; my_bar.progress(40) + + # vertical velocity field + data = read_variable('Vertical_velocity_pressure_isobaric', date_obj) + subdata['w700'] = data['Vertical_velocity_pressure_isobaric'].sel(isobaric1=70000.0, **sub_region).load() ; my_bar.progress(45) + + # pressure on pv surface + data = read_variable('Pressure_potential_vorticity_surface', date_obj) + subdata['pres_pv2'] = data['Pressure_potential_vorticity_surface'].sel(potential_vorticity_surface=2.0E-6, **sub_region).load() ; my_bar.progress(50) + subdata['pres_pv2'].metpy.convert_units('hPa') + + # geopotential height + data = read_variable('Geopotential_height_isobaric', date_obj) + subdata['gh200'] = data['Geopotential_height_isobaric'].sel(isobaric1=20000.0, **sub_region).load() ; my_bar.progress(55) + subdata['gh500'] = data['Geopotential_height_isobaric'].sel(isobaric1=50000.0, **sub_region).load() ; my_bar.progress(60) + subdata['gh700'] = data['Geopotential_height_isobaric'].sel(isobaric1=70000.0, **sub_region).load() ; my_bar.progress(62) + + # high temperature + data = read_variable('Temperature_isobaric', date_obj) + subdata['t500'] = data['Temperature_isobaric'].sel(isobaric1=50000.0, **sub_region).load() ; my_bar.progress(64) + subdata['t700'] = data['Temperature_isobaric'].sel(isobaric1=70000.0, **sub_region).load() ; my_bar.progress(66) + subdata['t850'] = data['Temperature_isobaric'].sel(isobaric1=85000.0, **sub_region).load() ; my_bar.progress(68) + subdata['t925'] = data['Temperature_isobaric'].sel(isobaric1=92500.0, **sub_region).load() ; my_bar.progress(70) + subdata['t500'] = subdata['t500'].metpy.convert_units('degC') + subdata['t700'] = subdata['t700'].metpy.convert_units('degC') + subdata['t850'] = subdata['t850'].metpy.convert_units('degC') + subdata['t925'] = subdata['t925'].metpy.convert_units('degC') + + # high moisture field + data = read_variable('Specific_humidity_isobaric', date_obj) + subdata['q700'] = data['Specific_humidity_isobaric'].sel(isobaric1=70000.0, **sub_region).load() ; my_bar.progress(75) + subdata['q850'] = data['Specific_humidity_isobaric'].sel(isobaric1=85000.0, **sub_region).load() ; my_bar.progress(80) + subdata['q925'] = data['Specific_humidity_isobaric'].sel(isobaric1=92500.0, **sub_region).load() ; my_bar.progress(85) + + # mean sea level pressure + data = read_variable('Pressure_msl', date_obj) + subdata['mslp'] = data['Pressure_msl'].sel(**sub_region).load() ; my_bar.progress(90) + subdata['mslp'] = subdata['mslp'].metpy.convert_units('hPa') + #subdata['mslp'] = subdata['mslp'].sortby('lat', ascending=True) ; 会造成mslp与其它变量的纬度信息排序不一致, 绘图出现矛盾 + + # precipitable water + data = read_variable('Precipitable_water_entire_atmosphere_single_layer', date_obj) + subdata['pwat'] = data['Precipitable_water_entire_atmosphere_single_layer'].sel(**sub_region).load() ; my_bar.progress(100) + + # surface temperature, 2018-, this data is wrong. + #data = read_variable('tsfc', date_obj) + #subdata['tsfc'] = data['tsfc'].sel(**sub_region).load() ; my_bar.progress(100) + #subdata['tsfc'].metpy.convert_units('degC') + else: + st.info('从10.28.49.118本地服务载入CFSR数据(0.5x0.5度, Waiting ~30s.)') + + # construct sub region + sub_region = {'lon':slice(map_region[0], map_region[1]), + 'lat':slice(map_region[2], map_region[3])} + + # Subset and load data + subdata = {} + + # wind field + my_bar = st.progress(0) + data = read_variable('u', date_obj) + subdata['u200'] = data['u'].sel(lev=200, **sub_region).load() ; my_bar.progress(4) + subdata['u500'] = data['u'].sel(lev=500, **sub_region).load() ; my_bar.progress(8) + subdata['u700'] = data['u'].sel(lev=700, **sub_region).load() ; my_bar.progress(12) + subdata['u850'] = data['u'].sel(lev=850, **sub_region).load() ; my_bar.progress(16) + subdata['u925'] = data['u'].sel(lev=925, **sub_region).load() ; my_bar.progress(20) + data = read_variable('v', date_obj) + subdata['v200'] = data['v'].sel(lev=200, **sub_region).load() ; my_bar.progress(24) + subdata['v500'] = data['v'].sel(lev=500, **sub_region).load() ; my_bar.progress(28) + subdata['v700'] = data['v'].sel(lev=700, **sub_region).load() ; my_bar.progress(32) + subdata['v850'] = data['v'].sel(lev=850, **sub_region).load() ; my_bar.progress(36) + subdata['v925'] = data['v'].sel(lev=925, **sub_region).load() ; my_bar.progress(40) + + # vertical velocity field + data = read_variable('w', date_obj) + subdata['w700'] = data['w'].sel(lev=700, **sub_region).load() ; my_bar.progress(45) + + # pressure on pv surface + data = read_variable('pres_pv', date_obj) + subdata['pres_pv2'] = data['pres_pv'].sel(lev=2.0E-6, **sub_region).load() ; my_bar.progress(50) + subdata['pres_pv2'].metpy.convert_units('hPa') + + # geopotential height + data = read_variable('g', date_obj) + subdata['gh200'] = data['g'].sel(lev=200, **sub_region).load() ; my_bar.progress(55) + subdata['gh500'] = data['g'].sel(lev=500, **sub_region).load() ; my_bar.progress(60) + subdata['gh700'] = data['g'].sel(lev=700, **sub_region).load() ; my_bar.progress(62) + + # high temperature + data = read_variable('t', date_obj) + subdata['t500'] = data['t'].sel(lev=500, **sub_region).load() ; my_bar.progress(64) + subdata['t700'] = data['t'].sel(lev=700, **sub_region).load() ; my_bar.progress(66) + subdata['t850'] = data['t'].sel(lev=850, **sub_region).load() ; my_bar.progress(68) + subdata['t925'] = data['t'].sel(lev=925, **sub_region).load() ; my_bar.progress(70) + subdata['t500'] = subdata['t500'].metpy.convert_units('degC') + subdata['t700'] = subdata['t700'].metpy.convert_units('degC') + subdata['t850'] = subdata['t850'].metpy.convert_units('degC') + subdata['t925'] = subdata['t925'].metpy.convert_units('degC') + + # high moisture field + data = read_variable('q', date_obj) + subdata['q700'] = data['q'].sel(lev=700, **sub_region).load() ; my_bar.progress(75) + subdata['q850'] = data['q'].sel(lev=850, **sub_region).load() ; my_bar.progress(80) + subdata['q925'] = data['q'].sel(lev=925, **sub_region).load() ; my_bar.progress(85) + + # mean sea level pressure + data = read_variable('pmsl', date_obj) + subdata['mslp'] = data['pmsl'].sel(**sub_region).load() ; my_bar.progress(90) + subdata['mslp'] = subdata['mslp'].metpy.convert_units('hPa') + #subdata['mslp'] = subdata['mslp'].sortby('lat', ascending=True) ; 会造成mslp与其它变量的纬度信息排序不一致, 绘图出现矛盾 + + # precipitable water + data = read_variable('pwat', date_obj) + subdata['pwat'] = data['pwat'].sel(**sub_region).load() ; my_bar.progress(100) + subdata['pwat'].metpy.convert_units('mm') + + # surface temperature, 2018-, this data is wrong. + #data = read_variable('tsfc', date_obj) + #subdata['tsfc'] = data['tsfc'].sel(**sub_region).load() ; my_bar.progress(100) + #subdata['tsfc'].metpy.convert_units('degC') + else: + st.info('从10.28.49.118本地服务载入ERA5数据(0.25x0.25度, Waiting ~30s.)') + + # construct sub region + sub_region = {'longitude':slice(map_region[0], map_region[1]), + 'latitude':slice(map_region[3], map_region[2])} + + # Subset and load data + subdata = {} + + # wind field + my_bar = st.progress(0) + data = read_variable('u_component_of_wind', date_obj, data_source) + subdata['u200'] = data['u'].sel(level=200, **sub_region).load().rename({'latitude':'lat', 'longitude':'lon'}) ; my_bar.progress(4) + subdata['u500'] = data['u'].sel(level=500, **sub_region).load().rename({'latitude':'lat', 'longitude':'lon'}) ; my_bar.progress(8) + subdata['u700'] = data['u'].sel(level=700, **sub_region).load().rename({'latitude':'lat', 'longitude':'lon'}) ; my_bar.progress(12) + subdata['u850'] = data['u'].sel(level=850, **sub_region).load().rename({'latitude':'lat', 'longitude':'lon'}) ; my_bar.progress(16) + subdata['u925'] = data['u'].sel(level=950, **sub_region).load().rename({'latitude':'lat', 'longitude':'lon'}) ; my_bar.progress(20) + data = read_variable('v_component_of_wind', date_obj, data_source) + subdata['v200'] = data['v'].sel(level=200, **sub_region).load().rename({'latitude':'lat', 'longitude':'lon'}) ; my_bar.progress(24) + subdata['v500'] = data['v'].sel(level=500, **sub_region).load().rename({'latitude':'lat', 'longitude':'lon'}) ; my_bar.progress(28) + subdata['v700'] = data['v'].sel(level=700, **sub_region).load().rename({'latitude':'lat', 'longitude':'lon'}) ; my_bar.progress(32) + subdata['v850'] = data['v'].sel(level=850, **sub_region).load().rename({'latitude':'lat', 'longitude':'lon'}) ; my_bar.progress(36) + subdata['v925'] = data['v'].sel(level=950, **sub_region).load().rename({'latitude':'lat', 'longitude':'lon'}) ; my_bar.progress(40) + + # vertical velocity field + data = read_variable('vertical_velocity', date_obj, data_source) + subdata['w700'] = data['w'].sel(level=700, **sub_region).load().rename({'latitude':'lat', 'longitude':'lon'}) ; my_bar.progress(45) + + + # geopotential height + data = read_variable('geopotential', date_obj, data_source) + subdata['gh200'] = data['z'].sel(level=200, **sub_region).load().rename({'latitude':'lat', 'longitude':'lon'}) ; my_bar.progress(55) + subdata['gh500'] = data['z'].sel(level=500, **sub_region).load().rename({'latitude':'lat', 'longitude':'lon'}) ; my_bar.progress(60) + subdata['gh700'] = data['z'].sel(level=700, **sub_region).load().rename({'latitude':'lat', 'longitude':'lon'}) ; my_bar.progress(62) + + # high temperature + data = read_variable('temperature', date_obj, data_source) + subdata['t500'] = data['t'].sel(level=500, **sub_region).load().rename({'latitude':'lat', 'longitude':'lon'}) ; my_bar.progress(64) + subdata['t700'] = data['t'].sel(level=700, **sub_region).load().rename({'latitude':'lat', 'longitude':'lon'}) ; my_bar.progress(66) + subdata['t850'] = data['t'].sel(level=850, **sub_region).load().rename({'latitude':'lat', 'longitude':'lon'}) ; my_bar.progress(68) + subdata['t925'] = data['t'].sel(level=950, **sub_region).load().rename({'latitude':'lat', 'longitude':'lon'}) ; my_bar.progress(70) + subdata['t500'] = subdata['t500'].metpy.convert_units('degC') + subdata['t700'] = subdata['t700'].metpy.convert_units('degC') + subdata['t850'] = subdata['t850'].metpy.convert_units('degC') + subdata['t925'] = subdata['t925'].metpy.convert_units('degC') + + # high moisture field + data = read_variable('specific_humidity', date_obj, data_source) + subdata['q700'] = data['q'].sel(level=700, **sub_region).load().rename({'latitude':'lat', 'longitude':'lon'}) ; my_bar.progress(75) + subdata['q850'] = data['q'].sel(level=850, **sub_region).load().rename({'latitude':'lat', 'longitude':'lon'}) ; my_bar.progress(80) + subdata['q925'] = data['q'].sel(level=950, **sub_region).load().rename({'latitude':'lat', 'longitude':'lon'}) ; my_bar.progress(85) + + # mean sea level pressure + data = read_variable('mslp', date_obj, data_source) + subdata['mslp'] = data['msl'].sel(**sub_region).load().rename({'latitude':'lat', 'longitude':'lon'}) ; my_bar.progress(90) + subdata['mslp'] = subdata['mslp'].metpy.convert_units('hPa') + + # precipitable water + data = read_variable('total_column_water', date_obj, data_source) + subdata['pwat'] = data['tcw'].sel(**sub_region).load().rename({'latitude':'lat', 'longitude':'lon'}) ; my_bar.progress(100) + #subdata['pwat'].metpy.convert_units('mm') + subdata['pwat'].attrs['units']='mm' + + return subdata + + +if __name__ == "__main__": + main() + diff --git a/apps/reanalysis_maps/draw_maps.py b/apps/reanalysis_maps/draw_maps.py index e642d89..4584658 100644 --- a/apps/reanalysis_maps/draw_maps.py +++ b/apps/reanalysis_maps/draw_maps.py @@ -56,17 +56,27 @@ def draw_observation(data, date_obj, map_region): keys = ['0.1~10', '10~25', '25~50', '50~100', '100~250', '>=250'] cols = ['lightgreen', 'yellow', 'lightskyblue', 'blue', 'magenta','maroon'] cols_map = dict(zip(keys, cols)) - data['rain'] = pd.cut(data['PRE_Time_0808'], bins=bins, labels=keys) - data['Rainfall'] = '['+data['Lon'].round(2).astype(str) + ',' + data['Lat'].round(2).astype(str) + ']: ' + \ + if 'PRE_Time_0808' in data: + data['rain'] = pd.cut(data['PRE_Time_0808'], bins=bins, labels=keys) + data['Rainfall'] = '['+data['Lon'].round(2).astype(str) + ',' + data['Lat'].round(2).astype(str) + ']: ' + \ data['PRE_Time_0808'].astype(str) - data['rain_size'] = data['PRE_Time_0808'] + data['PRE_Time_0808'].mean() + data['rain_size'] = data['PRE_Time_0808'] + data['PRE_Time_0808'].mean() + title_tp = 'Accumulated precipitation ({}'.format(date_obj.strftime("%Y%m%d08 - "))+'{})'.format((date_obj+dt.timedelta(days=1)).strftime("%Y%m%d08")) + if 'PRE_Time_2020' in data: + data['rain'] = pd.cut(data['PRE_Time_2020'], bins=bins, labels=keys) + data['Rainfall'] = '['+data['Lon'].round(2).astype(str) + ',' + data['Lat'].round(2).astype(str) + ']: ' + \ + data['PRE_Time_2020'].astype(str) + data['rain_size'] = data['PRE_Time_2020'] + data['PRE_Time_2020'].mean() + title_tp = 'Accumulated precipitation ({}'.format((date_obj-dt.timedelta(days=1)).strftime("%Y%m%d20 - "))+'{})'.format(date_obj.strftime("%Y%m%d20")) + df = data[data['rain'].notna()] if df.shape[0] >= 2: figs['Rainfall'] = px.scatter_mapbox( df, lat="Lat", lon="Lon", color="rain", category_orders={'rain': keys}, color_discrete_map = cols_map, hover_data={'Rainfall':True, 'Lon':False, 'Lat':False, 'rain':False, 'rain_size':False}, mapbox_style='satellite-streets', size="rain_size", center=map_center, size_max=10, zoom=4, - title = 'Accumulated precipitation ({})'.format(date_obj.strftime("%Y%m%d 08-08")), +# title = 'Accumulated precipitation ({})'.format(date_obj.strftime("%Y%m%d 08-08")), + title = title_tp, width=900, height=700) # draw maximum temperature diff --git a/apps/reanalysis_maps/draw_maps_era5.py b/apps/reanalysis_maps/draw_maps_era5.py new file mode 100644 index 0000000..b37feff --- /dev/null +++ b/apps/reanalysis_maps/draw_maps_era5.py @@ -0,0 +1,544 @@ +# _*_ coding: utf-8 _*_ + +# Copyright (c) 2020 NMC Developers. +# Distributed under the terms of the GPL V3 License. + + +""" +Synoptic Composite Ploting script + +refer to https://github.com/tomerburg/python_gallery +""" + + +#Import the necessary libraries +import collections +import numpy as np +import xarray as xr +import pandas as pd +import datetime as dt +import scipy.ndimage as ndimage +import matplotlib.pyplot as plt +import matplotlib.colors as col +import matplotlib.gridspec as gridspec + +from cartopy import crs as ccrs +import cartopy.feature as cfeature +from cartopy import util as cu + +import nmc_met_io.config as CONFIG +import plotly.express as px + +import metpy.calc as calc +from metpy.units import units + +from nmc_met_graphics.cmap.cm import gradient +from nmc_met_graphics.plot.util import add_mslp_label +from nmc_met_graphics.plot.mapview import add_china_map_2cartopy +from nmc_met_graphics.magics import dynamics, thermal, pv, moisture + + +def draw_observation(data, date_obj, map_region): + """ + Draw observation map with plotly + """ + + # set mapbox token + px.set_mapbox_access_token(CONFIG.CONFIG['MAPBOX']['token']) + + # create figures + map_center = {'lat':(map_region[2] + map_region[3]) * 0.5, + 'lon':(map_region[0] + map_region[1]) * 0.5} + figs = collections.OrderedDict() + + # draw precipitation + bins = [0.1, 10, 25, 50, 100, 250, 1200] + keys = ['0.1~10', '10~25', '25~50', '50~100', '100~250', '>=250'] + cols = ['lightgreen', 'yellow', 'lightskyblue', 'blue', 'magenta','maroon'] + cols_map = dict(zip(keys, cols)) + if 'PRE_Time_0808' in data: + data['rain'] = pd.cut(data['PRE_Time_0808'], bins=bins, labels=keys) + data['Rainfall'] = '['+data['Lon'].round(2).astype(str) + ',' + data['Lat'].round(2).astype(str) + ']: ' + \ + data['PRE_Time_0808'].astype(str) + data['rain_size'] = data['PRE_Time_0808'] + data['PRE_Time_0808'].mean() + title_tp = 'Accumulated precipitation ({}'.format(date_obj.strftime("%Y%m%d08 - "))+'{})'.format((date_obj+dt.timedelta(days=1)).strftime("%Y%m%d08")) + if 'PRE_Time_2020' in data: + data['rain'] = pd.cut(data['PRE_Time_2020'], bins=bins, labels=keys) + data['Rainfall'] = '['+data['Lon'].round(2).astype(str) + ',' + data['Lat'].round(2).astype(str) + ']: ' + \ + data['PRE_Time_2020'].astype(str) + data['rain_size'] = data['PRE_Time_2020'] + data['PRE_Time_2020'].mean() + title_tp = 'Accumulated precipitation ({}'.format((date_obj-dt.timedelta(days=1)).strftime("%Y%m%d20 - "))+'{})'.format(date_obj.strftime("%Y%m%d20")) + + df = data[data['rain'].notna()] + if df.shape[0] >= 2: + figs['Rainfall'] = px.scatter_mapbox( + df, lat="Lat", lon="Lon", color="rain", category_orders={'rain': keys}, color_discrete_map = cols_map, + hover_data={'Rainfall':True, 'Lon':False, 'Lat':False, 'rain':False, 'rain_size':False}, + mapbox_style='satellite-streets', size="rain_size", center=map_center, size_max=10, zoom=4, +# title = 'Accumulated precipitation ({})'.format(date_obj.strftime("%Y%m%d 08-08")), + title = title_tp, + width=900, height=700) + + # draw maximum temperature + bins = [35, 37, 40, 60] + keys = ['35~37', '37~40', '>=40'] + cols = ['rgb(255,191,187)', 'rgb(250,89,0)', 'rgb(230,0,8)'] + cols_map = dict(zip(keys, cols)) + data['max_temp_warning'] = pd.cut(data['TEM_Max'], bins=bins, labels=keys) + data['max_temp'] = '['+data['Lon'].round(2).astype(str) + ',' + data['Lat'].round(2).astype(str) + ']: ' + \ + data['TEM_Max'].astype(str) + df = data[data['max_temp_warning'].notna()] + if df.shape[0] >= 2: + figs['Max_temperature'] = px.scatter_mapbox( + df, lat="Lat", lon="Lon", color="max_temp_warning", category_orders={'max_temp_warning': keys}, + color_discrete_map = cols_map, + hover_data={'max_temp':True, 'Lon':False, 'Lat':False, 'max_temp_warning':False, 'TEM_Max':False}, + mapbox_style='satellite-streets', size="TEM_Max", center=map_center, size_max=10, zoom=4, + title = 'Maximum temperature ({})'.format(date_obj.strftime("%Y%m%d 08-08")), + width=900, height=700) + + # draw minimum temperature + bins = [-120, -40, -30, -20, -10, 0] + keys = ['<=-40','-40~-30', '-30~-20', '-20~-10', '-10~0'] + cols = ['rgb(178,1,223)', 'rgb(8,7,249)', 'rgb(5,71,162)', 'rgb(5,109,250)', 'rgb(111,176,248)'] + cols_map = dict(zip(keys, cols)) + data['min_temp_warning'] = pd.cut(data['TEM_Min'], bins=bins, labels=keys) + data['min_temp'] = '['+data['Lon'].round(2).astype(str) + ',' + data['Lat'].round(2).astype(str) + ']: ' + \ + data['TEM_Min'].astype(str) + df = data[data['min_temp_warning'].notna()] + if df.shape[0] >= 2: + figs['Min_temprature'] = px.scatter_mapbox( + df, lat="Lat", lon="Lon", color="min_temp_warning", category_orders={'min_temp_warning': keys}, + color_discrete_map = cols_map, + hover_data={'min_temp':True, 'Lon':False, 'Lat':False, 'min_temp_warning':False, 'TEM_Min':False}, + mapbox_style='satellite-streets', size=-1.0*df["TEM_Min"], center=map_center, size_max=10, zoom=4, + title = 'Minimum temperature ({})'.format(date_obj.strftime("%Y%m%d 08-08")), + width=900, height=700) + + # draw low visibility + data['VIS_Min'] /= 1000.0 + bins = [0, 0.05, 0.2, 0.5, 1] + keys = ['<=0.05','0.05~0.2', '0.2~0.5', '0.5~1'] + cols = ['rgb(0,82,77)', 'rgb(0,153,160)', 'rgb(0,210,204)', 'rgb(95,255,252)'] + cols_map = dict(zip(keys, cols)) + data['min_vis_warning'] = pd.cut(data['VIS_Min'], bins=bins, labels=keys) + data['VIS_Min_size'] = 2.0-data["VIS_Min"] + data['min_vis'] = '['+data['Lon'].round(2).astype(str) + ',' + data['Lat'].round(2).astype(str) + ']: ' + \ + data['VIS_Min'].astype(str) + df = data[data['min_vis_warning'].notna()] + if df.shape[0] >= 2: + figs['Low_visibility'] = px.scatter_mapbox( + df, lat="Lat", lon="Lon", color="min_vis_warning", category_orders={'min_vis_warning': keys}, + color_discrete_map = cols_map, + hover_data={'min_vis':True, 'Lon':False, 'Lat':False, 'min_vis_warning':False, 'VIS_Min_size':False}, + mapbox_style='satellite-streets', size="VIS_Min_size", center=map_center, size_max=10, zoom=4, + title = 'Low visibility ({})'.format(date_obj.strftime("%Y%m%d 08-08")), + width=900, height=700) + + # draw high wind + bins = [10.8, 13.9, 17.2, 20.8, 24.5, 28.5, 32.7, 37.0, 120] + keys = ['10.8~13.8','13.9~17.1', '17.2~20.7', '20.8~24.4', '24.5~28.4', '28.5~32.6', '32.7~36.9', '>=37.0'] + cols = ['rgb(0,210,244)', 'rgb(0,125,255)', 'rgb(253,255,0)', 'rgb(247,213,0)', + 'rgb(255,141,0)', 'rgb(251,89,91)', 'rgb(255,3,0)', 'rgb(178,1,223)'] + cols_map = dict(zip(keys, cols)) + data['max_win_warning'] = pd.cut(data['WIN_S_Max'], bins=bins, labels=keys) + data['max_win'] = '['+data['Lon'].round(2).astype(str) + ',' + data['Lat'].round(2).astype(str) + ']: ' + \ + data['WIN_S_Max'].astype(str) + df = data[data['max_win_warning'].notna()] + if df.shape[0] >= 2: + figs['High_wind'] = px.scatter_mapbox( + df, lat="Lat", lon="Lon", color="max_win_warning", category_orders={'max_win_warning': keys}, + color_discrete_map = cols_map, + hover_data={'max_win':True, 'Lon':False, 'Lat':False, 'max_win_warning':False, 'WIN_S_Max':False}, + mapbox_style='satellite-streets', size="WIN_S_Max", center=map_center, size_max=10, zoom=4, + title = 'Maximum wind speed ({})'.format(date_obj.strftime("%Y%m%d 08-08")), + width=1000, height=800) + + return figs + + +def draw_weather_analysis(date_obj, data, map_region, return_dict, data_source='CFSR'): + """ + Draw weather analysis map. + """ + + # image dictionary + images = collections.OrderedDict() + return_dict[0] = None + + if data_source == 'CFSR': + # draw 2PVU surface pressure + image = pv.draw_pres_pv2( + data['pres_pv2'].values, data['pres_pv2']['lon'].values, data['pres_pv2']['lat'].values, + map_region=map_region, title_kwargs={'name':data_source, 'time': date_obj}) + images['2PVU_Surface_Pressure'] = image + + # draw 200hPa wind field + image = dynamics.draw_wind_upper( + data['u200'].values, data['v200'].values, + data['u200']['lon'].values, data['u200']['lat'].values, + gh=data['gh200'].values, map_region=map_region, + title_kwargs={'name':data_source, 'head': "200hPa Wind | GH", 'time': date_obj}) + images['200hPa_Wind'] = image + + # draw 500hPa height and temperature + image = dynamics.draw_height_temp( + data['gh500'].values, data['t500'].values, + data['gh500']['lon'].values, data['gh500']['lat'].values, map_region=map_region, + title_kwargs={'name':data_source, 'head': "500hPa GH | T", 'time': date_obj}) + images['500hPa_Height'] = image + + # draw 500hPa vorticity + image = dynamics.draw_vort_high( + data['u500'].values, data['v500'].values, + data['u500']['lon'].values, data['u500']['lat'].values, + gh=data['gh500'].values, map_region=map_region, + title_kwargs={'name':data_source, 'head': "500hPa Wind | Vorticity | GH", 'time': date_obj}) + images['500hPa_Vorticity'] = image + + # draw 700hPa vertical velocity + image = dynamics.draw_vvel_high( + data['u700'].values, data['v700'].values, data['w700'].values, + data['w700']['lon'].values, data['w700']['lat'].values, + gh=data['gh700'].values, map_region=map_region, + title_kwargs={'name':data_source, 'head': "700hPa Vertical Velocity | Wind | GH", 'time': date_obj}) + images['700hPa_Vertical_Velocity'] = image + + # draw 700hPa wind field + image = dynamics.draw_wind_high( + data['u700'].values, data['v700'].values, + data['u700']['lon'].values, data['u700']['lat'].values, + gh=data['gh500'].values, map_region=map_region, + title_kwargs={'name':data_source, 'head': "700hPa Wind | 500hPa GH", 'time': date_obj}) + images['700hPa_Wind'] = image + + # draw 700hPa temperature field + image = thermal.draw_temp_high( + data['t700'].values, data['t700']['lon'].values, data['t700']['lat'].values, + gh=data['gh500'].values, map_region=map_region, + title_kwargs={'name':data_source, 'head': "700hPa T | 500hPa GH", 'time': date_obj}) + images['700hPa_Temperature'] = image + + # draw 700hPa relative humidity + rh = calc.relative_humidity_from_specific_humidity(700 * units.hPa, data['t700'], data['q700']) * 100 + image = moisture.draw_rh_high( + data['u700'].values, data['v700'].values, rh.values, + data['u700']['lon'].values, data['u700']['lat'].values, + gh=data['gh500'].values, map_region=map_region, + title_kwargs={'name':data_source, 'head': "700hPa RH | Wind | 500hPa GH", 'time': date_obj}) + images['700hPa_Relative_Humidity'] = image + + # draw 850hPa wind field + image = dynamics.draw_wind_high( + data['u850'].values, data['v850'].values, + data['u850']['lon'].values, data['u850']['lat'].values, + gh=data['gh500'].values, map_region=map_region, + title_kwargs={'name':data_source, 'head': "850hPa Wind | 500hPa GH", 'time': date_obj}) + images['850hPa_Wind'] = image + + # draw 850hPa temperature field + image = thermal.draw_temp_high( + data['t850'].values, data['t850']['lon'].values, data['t850']['lat'].values, + gh=data['gh500'].values, map_region=map_region, + title_kwargs={'name':data_source, 'head': "850hPa T | 500hPa GH", 'time': date_obj}) + images['850hPa_Temperature'] = image + + # draw 850hPa relative humidity + rh = calc.relative_humidity_from_specific_humidity(850 * units.hPa, data['t850'], data['q850']) * 100 + image = moisture.draw_rh_high( + data['u850'].values, data['v850'].values, rh.values, + data['u850']['lon'].values, data['u850']['lat'].values, + gh=data['gh500'].values, map_region=map_region, + title_kwargs={'name':data_source, 'head': "850hPa RH | Wind | 500hPa GH", 'time': date_obj}) + images['850hPa_Relative_Humidity'] = image + + # draw 850hPa specific field + image = moisture.draw_sp_high( + data['u850'].values, data['v850'].values, data['q850'].values*1000., + data['q850']['lon'].values, data['q850']['lat'].values, + gh=data['gh500'].values, map_region=map_region, + title_kwargs={'name':data_source, 'head': "850hPa SP | Wind | 500hPa GH", 'time': date_obj}) + images['850hPa_Specific_Humidity'] = image + + # draw 925hPa temperature field + image = thermal.draw_temp_high( + data['t925'].values, data['t925']['lon'].values, data['t925']['lat'].values, + gh=data['gh500'].values, map_region=map_region, + title_kwargs={'name':data_source, 'head': "925hPa T | 500hPa GH", 'time': date_obj}) + images['925hPa_Temperature'] = image + + # draw 925hPa wind field + image = dynamics.draw_wind_high( + data['u925'].values, data['v925'].values, + data['u925']['lon'].values, data['u925']['lat'].values, + gh=data['gh500'].values, map_region=map_region, + title_kwargs={'name':data_source, 'head': "925hPa Wind | 500hPa GH", 'time': date_obj}) + images['925hPa_Wind'] = image + + # draw 925hPa relative humidity + rh = calc.relative_humidity_from_specific_humidity(925 * units.hPa, data['t925'], data['q925']) * 100 + image = moisture.draw_rh_high( + data['u925'].values, data['v925'].values, rh.values, + data['u925']['lon'].values, data['u925']['lat'].values, + gh=data['gh500'].values, map_region=map_region, + title_kwargs={'name':data_source, 'head': "925hPa RH | Wind | 500hPa GH", 'time': date_obj}) + images['925hPa_Relative_Humdity'] = image + + # draw 925hPa specific field + image = moisture.draw_sp_high( + data['u925'].values, data['v925'].values, data['q925'].values*1000., + data['q925']['lon'].values, data['q925']['lat'].values, + gh=data['gh500'].values, map_region=map_region, + title_kwargs={'name':data_source, 'head': "925hPa SP | Wind | 500hPa GH", 'time': date_obj}) + images['925hPa_Specific_Humidity'] = image + + # draw precipitable water field + image = moisture.draw_pwat( + data['pwat'].values, data['pwat']['lon'].values, data['pwat']['lat'].values, + gh=data['gh500'].values, map_region=map_region, + title_kwargs={'name':data_source, 'head': "Precipitable Water | 500hPa GH", 'time': date_obj}) + images['Precipitable_Water'] = image + + # draw mean sea level pressure field + image = dynamics.draw_mslp( + data['mslp'].values, data['mslp']['lon'].values, data['mslp']['lat'].values, + gh=data['gh500'].values, map_region=map_region, + title_kwargs={'name':data_source, 'head': "MSLP | 500hPa GH", 'time': date_obj}) + images['Mean_Sea_Level_Pressure'] = image + + return_dict[0] = images + + +#Spatially smooth a 2D variable +def smooth(prod,sig): + + #Check if variable is an xarray dataarray + try: + lats = prod.lat.values + lons = prod.lon.values + prod = ndimage.gaussian_filter(prod,sigma=sig,order=0) + prod = xr.DataArray(prod, coords=[lats, lons], dims=['lat', 'lon']) + except: + prod = ndimage.gaussian_filter(prod,sigma=sig,order=0) + + return prod + + +def draw_composite_map(date_obj, t850, u200, v200, u500, v500, mslp, gh500, u850, v850, pwat): + """ + Draw synoptic composite map. + All variables must have the same region. + + Args: + date_obj (datetime) : datetime ojbect. + t850 (xarray): 850hPa temperature, must have lat and lon coordinates. + u200, v200 (xarray): 200hPa u and v wind component. + u500, v500 (xarray): 500hPa u and v wind component. + mslp (xarray): mean sea level pressure. + gh500 (xarray): 500hPa geopotential height. + u850, v850 (xarray): 850hPa u and v wind component. + pwat (xarray): precipitable water. + """ + + #Get lat and lon arrays for this dataset: + lat = t850.lat.values + lon = t850.lon.values + + #======================================================================================================== + # Create a Basemap plotting figure and add geography + #======================================================================================================== + + #Create a Plate Carree projection object + proj_ccrs = ccrs.Miller(central_longitude=0.0) + + #Create figure and axes for main plot and colorbars + fig = plt.figure(figsize=(18,12),dpi=125) + gs = gridspec.GridSpec(12, 36, figure=fig) #[ytop:ybot, xleft:xright] + ax = plt.subplot(gs[:, :-1],projection=proj_ccrs) #main plot + ax.set_xticklabels([]) + ax.set_yticklabels([]) + ax2 = plt.subplot(gs[:4, -1]) #top plot + ax2.set_xticklabels([]) + ax2.set_yticklabels([]) + ax3 = plt.subplot(gs[4:8, -1]) #bottom plot + ax3.set_xticklabels([]) + ax3.set_yticklabels([]) + ax4 = plt.subplot(gs[8:, -1]) #bottom plot + ax4.set_xticklabels([]) + ax4.set_yticklabels([]) + + #Add political boundaries and coastlines + ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidths=1.2) + ax.add_feature(cfeature.BORDERS.with_scale('50m'), linewidths=1.2) + ax.add_feature(cfeature.STATES.with_scale('50m'), linewidths=0.5) + + #Add land/lake/ocean masking + land_mask = cfeature.NaturalEarthFeature('physical', 'land', '50m', + edgecolor='face', facecolor='#e6e6e6') + sea_mask = cfeature.NaturalEarthFeature('physical', 'ocean', '50m', + edgecolor='face', facecolor='#ffffff') + lake_mask = cfeature.NaturalEarthFeature('physical', 'lakes', '50m', + edgecolor='face', facecolor='#ffffff') + ax.add_feature(sea_mask,zorder=0) + ax.add_feature(land_mask,zorder=0) + ax.add_feature(lake_mask,zorder=0) + + #======================================================================================================== + # Fill contours + #======================================================================================================== + + #-------------------------------------------------------------------------------------------------------- + # 850-hPa temperature + #-------------------------------------------------------------------------------------------------------- + + #Specify contour settings + clevs = np.arange(-40,40,1) + cmap = plt.get_cmap('jet') + extend = "both" + + #Contour fill this variable + norm = col.BoundaryNorm(clevs,cmap.N) + cs = ax.contourf(lon,lat,t850,clevs,cmap=cmap,norm=norm,extend=extend,transform=proj_ccrs,alpha=0.1) + + #-------------------------------------------------------------------------------------------------------- + # PWAT + #-------------------------------------------------------------------------------------------------------- + + #Specify contour settings + clevs = np.arange(20,71,0.5) + + #Define a color gradient for PWAT + pwat_colors = gradient([[(255,255,255),0.0],[(255,255,255),20.0]], + [[(205,255,205),20.0],[(0,255,0),34.0]], + [[(0,255,0),34.0],[(0,115,0),67.0]]) + cmap = pwat_colors.get_cmap(clevs) + extend = "max" + + #Contour fill this variable + norm = col.BoundaryNorm(clevs,cmap.N) + cs = ax.contourf(lon,lat,pwat,clevs,cmap=cmap,norm=norm,extend=extend,transform=proj_ccrs,alpha=0.9) + + #Add a color bar + _ = plt.colorbar(cs,cax=ax2,shrink=0.75,pad=0.01,ticks=[20,30,40,50,60,70]) + + #-------------------------------------------------------------------------------------------------------- + # 250-hPa wind + #-------------------------------------------------------------------------------------------------------- + + #Get the data for this variable + wind = calc.wind_speed(u200, v200) + + #Specify contour settings + clevs = [40,50,60,70,80,90,100,110] + cmap = col.ListedColormap(['#99E3FB','#47B6FB','#0F77F7','#AC97F5','#A267F4','#9126F5','#E118F3','#E118F3']) + extend = "max" + + #Contour fill this variable + norm = col.BoundaryNorm(clevs,cmap.N) + cs = ax.contourf(lon,lat,wind,clevs,cmap=cmap,norm=norm,extend=extend,transform=proj_ccrs) + + #Add a color bar + _ = plt.colorbar(cs,cax=ax3,shrink=0.75,pad=0.01,ticks=clevs) + + #-------------------------------------------------------------------------------------------------------- + # 500-hPa smoothed vorticity + #-------------------------------------------------------------------------------------------------------- + + #Get the data for this variable + dx,dy = calc.lat_lon_grid_deltas(lon,lat) + vort = calc.vorticity(u500, v500, dx=dx, dy=dy) + smooth_vort = smooth(vort, 5.0) * 10**5 + + #Specify contour settings + clevs = np.arange(2,20,1) + cmap = plt.get_cmap('autumn_r') + extend = "max" + + #Contour fill this variable + norm = col.BoundaryNorm(clevs,cmap.N) + cs = ax.contourf(lon,lat,smooth_vort,clevs,cmap=cmap,norm=norm,extend=extend,transform=proj_ccrs,alpha=0.3) + + #Add a color bar + _ = plt.colorbar(cs,cax=ax4,shrink=0.75,pad=0.01,ticks=clevs[::2]) + + #======================================================================================================== + # Contours + #======================================================================================================== + + #-------------------------------------------------------------------------------------------------------- + # MSLP + #-------------------------------------------------------------------------------------------------------- + + #Specify contour settings + clevs = np.arange(960,1040+4,4) + style = 'solid' #Plot solid lines + color = 'red' #Plot lines as gray + width = 0.8 #Width of contours 0.25 + + #Contour this variable + cs = ax.contour(lon,lat,mslp,clevs,colors=color,linewidths=width,linestyles=style,transform=proj_ccrs,alpha=0.9) + + #Include value labels + ax.clabel(cs, inline=1, fontsize=9, fmt='%d') + + #-------------------------------------------------------------------------------------------------------- + # Geopotential heights + #-------------------------------------------------------------------------------------------------------- + + #Get the data for this variable + gh500 = gh500 / 10.0 + + #Specify contour settings + clevs = np.arange(480,612,4) + style = 'solid' #Plot solid lines + color = 'black' #Plot lines as gray + width = 2.0 #Width of contours + + #Contour this variable + cs = ax.contour(lon,lat,gh500,clevs,colors=color,linewidths=width,linestyles=style,transform=proj_ccrs) + + #Include value labels + ax.clabel(cs, inline=1, fontsize=12, fmt='%d') + + #-------------------------------------------------------------------------------------------------------- + # Surface barbs + #-------------------------------------------------------------------------------------------------------- + + #Plot wind barbs + _ = ax.quiver(lon, lat, u850.values, v850.values, transform=proj_ccrs, regrid_shape=(38,30), scale=820, alpha=0.5) + + #-------------------------------------------------------------------------------------------------------- + # Label highs & lows + #-------------------------------------------------------------------------------------------------------- + + #Label highs and lows + add_mslp_label(ax, proj_ccrs, mslp, lat, lon) + + #======================================================================================================== + # Step 6. Add map boundary, legend, plot title, then save image and close + #======================================================================================================== + + #Add china province boundary + add_china_map_2cartopy(ax, name='province') + + #Add custom legend + from matplotlib.lines import Line2D + custom_lines = [Line2D([0], [0], color='#00A123', lw=5), + Line2D([0], [0], color='#0F77F7', lw=5), + Line2D([0], [0], color='#FFC000', lw=5), + Line2D([0], [0], color='k', lw=2), + Line2D([0], [0], color='k', lw=0.1, marker=r'$\rightarrow$', ms=20), + Line2D([0], [0], color='r', lw=0.8),] + + ax.legend(custom_lines, ['PWAT (mm)', '200-hPa Wind (m/s)', '500-hPa Vorticity', '500-hPa Height (dam)', '850-hPa Wind (m/s)', 'MSLP (hPa)'], loc=2, prop={'size':12}) + + #Format plot title + title = "Synoptic Composite \nValid: " + dt.datetime.strftime(date_obj,'%Y-%m-%d %H%M UTC') + st = plt.suptitle(title,fontweight='bold',fontsize=16) + st.set_y(0.92) + + #Return figuration + return(fig) + + +