- 积分
- 63453
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2011-7-23
- 最后登录
- 1970-1-1
|
发表于 2024-11-13 10:33:00
|
显示全部楼层
- """
- Spyder Editor
- OLDLee.
- """
- ###引用部分
- import numpy as np
- import pandas as pd
- from PIL import Image
- from metpy.calc import wind_components
- from metpy.plots import StationPlot
- from metpy.plots.wx_symbols import current_weather, sky_cover
- from metpy.units import units
- import cartopy.crs as ccrs
- import cartopy.feature as feat
- import matplotlib.pyplot as plt
- import cartopy.io.shapereader as shpreader
- plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文
- plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
- ## 数据读取 与 整理
- a=r"data/Metpy/Surf14090120.000"
- df=pd.read_csv(a,skiprows=3,header=None,encoding="GB2312",sep='\s+').drop([0])
- df1=df.iloc[::2].reset_index(drop=True)
- df2=df.iloc[1::2].reset_index(drop=True)
- df1=df1.drop(columns=[12,13])
- head=["stid","lon","lat","h","lev","cloud_fraction","wind_dir",
- "wind_speed","slp","p3","w1","w2","r6","lc","lcc","lch",
- "dew_point_temperature","vv","weather","air_temperature","mc","hc","s1","s2","T24","P24"]
- res=pd.concat([df1,df2],axis=1,ignore_index=True)
- res.columns=head
- res[['stid','cloud_fraction','weather','w1','w2','lcc','lc','lch','mc','hc','p3','r6']]=res[['stid','cloud_fraction','weather','w1','w2','lcc','lc','lch','mc','hc','p3','r6']].astype('int')
- res=res[res['stid']<80000]
- res['weather']=res['weather'].replace(9999,0)
- res['wind_speed']=res['wind_speed'].replace(9999,0)
- res['weather']=res['weather'].astype('int')
- res['cloud_fraction']=res['cloud_fraction'].replace(9999,10)
- ## 数据筛选 与 处理
- all_data = res[['stid','lat','lon','slp','air_temperature','cloud_fraction','dew_point_temperature','weather','wind_dir','wind_speed','p3','r6']]
- #挑选目标区域内的站点站号白名单
- with open ('data/Metpy/whitelist.txt') as f:
- whitelist=[int(id) for id in list(f)]
- # 找出数据和白名单交集的站号
- data=all_data[all_data['stid'].isin(whitelist)]
- ###########################################
- # 现在有了需要的数据, 需要做一些转换:
- #
- # - 获得所需站点站号的字符串列表 Get a list of strings for the station IDs
- # - 由风速、风向，计算U、V分量 Get wind components from speed and direction
- # - 转换总云量为电码** Convert cloud fraction values to integer codes [0 - 8]
- # - 映射天气现象电码，以匹配符号 Map METAR weather codes to WMO codes for weather symbols
- p3 = data['p3']
- p3.where(p3>999,np.nan,inplace=True)
- r6 = data['r6']
- r6.where(r6>150,np.nan,inplace=True)
- r6.where(r6<=0,np.nan,inplace=True)
- data['slp'][data['slp']>9000]=np.NaN
- stids =[str(id) for id in list(data['stid'].values)]
- u, v = wind_components((data['wind_speed'].values * units('m/s')),data['wind_dir'].values * units.degree)
- cloud_frac = (data['cloud_fraction']).astype(int)
- wx_text = [s for s in data['weather']]
- li=np.arange(0,100)
- wx_codes = dict(zip(list(map(str,li)),li))
- wx = [wx_codes[s.split()[0] if ' ' in str(s) else str(s)] for s in wx_text]
- ## 绘图部分
- proj = ccrs.PlateCarree()
- fig = plt.figure(figsize=(20,14),dpi=200)
- ax = fig.add_subplot(1, 1, 1, frameon=False,projection=ccrs.LambertConformal(central_longitude=80,central_latitude=35.))
- ax.set_extent((70, 100, 30, 50))
- ax.add_feature(feat.LAKES, zorder=1)
- ax.add_feature(feat.LAND, zorder=-1)
- ax.add_feature(feat.OCEAN, zorder=-1)
- ax.add_feature(feat.RIVERS, zorder=1)
- # ax.coastlines(resolution='110m', zorder=2, color='black')
- # ax.add_feature(feat.BORDERS, linewidth='2', edgecolor='black')
- ##### 声明地面填图绘图区系列
- stationplot = StationPlot(ax, data['lon'], data['lat'], transform=proj,fontsize=4.5)
- # 在左上角绘制温度（红色），左下角绘制露点（绿色）
- data['air_temperature'][data['air_temperature']>9000]=np.NaN
- data['dew_point_temperature'][data['dew_point_temperature']>9000]=np.NaN
- stationplot.plot_parameter('NW', data['air_temperature'], color='red',transform=proj)
- stationplot.plot_parameter('SW', data['dew_point_temperature'], color='darkgreen',transform=proj)
- # 按照天气图填图规范，填注三位数的海平面气压：
- # 1000hPa以下减900，再乘以10；1000以上减1000，再乘以10
- stationplot.plot_parameter('NE', data['slp'],
- formatter=lambda v: format(1 * v, '.0f')[-4:],transform=proj)
- # 在站点右上角填注3小时变压
- stationplot.plot_parameter('E', data['p3'],transform=proj)
- # 在站点右下角填注6小时降水
- stationplot.plot_parameter('SE', r6,transform=proj)
- # 在居中位置绘制总云量符号。将总云量成数，转为MetPy定义的sky_cover。
- stationplot.plot_symbol('C', cloud_frac, sky_cover,transform=proj)
- # 将站点中心视为原点，在(-1,0)位置，也就是正西方向，绘制现在天气符号；在(2,-1)右下角位置，填写站号
- stationplot.plot_symbol((-1,0), wx, current_weather,transform=proj,fontsize=5)
- stationplot.plot_text((2,-1), stids,transform=proj,fontsize=2.5)
- # 添加风羽
- stationplot.plot_barb(u, v,transform=proj,barb_increments={'half':2, 'full':4, 'flag':20},sizes={'emptybarb':0.0})
- ax.gridlines(linestyle='--',
- xlocs=[60,70,80,90,100,110,120,130,140,150],
- ylocs=[10,20,30,40,50,60,70])
- # 添加一些基本的地理信息
- shpborder='shp/NationalBorder.shp'
- adm1_shapes = list(shpreader.Reader(shpborder).geometries())
- ax.add_geometries(adm1_shapes[:],ccrs.PlateCarree(),edgecolor='k',facecolor='None',linewidth=0.8)
- ### 这一段的意义在于：
- ### 注意这个x和y的值很大，是最后一个“局地绘图区”的相对位置！！！
- ### 如后续例图所示
- # #添加左下角的文本、矩形框
- # rect = plt.Rectangle((-2650000,-1800000),1450000,850000,
- # facecolor ='w',edgecolor='k',alpha=0.9)
- # ax.add_patch(rect)
- # ax.text(-2600000, -1100000,'地 面 天 气 图',fontsize=20)
- # ax.text(-2600000, -1300000,year+' 年 '+month+'月 '+day+'日 '
- # +valid+'时',fontsize=20)
- # ax.text(-2600000, -1500000,'分析预报员：__________',fontsize=20)
- ## 预览输出
- plt.show()
|
|
窥视卡
雷达卡
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
楼主