读mdfs数据并画图

注：解决的几项技术细节问题：

1、cartopy叠加shp地图文件；
2、scipy.rbf 从站点插值到格点 一定要注意不能出现一个经纬度点出现多个值，否则会报错；
3、cartopy画图调整colorbar、调整线条颜色、阴影颜色、等值线label等细节问题。

mdfs(Ming distributed file system)是一个负载均衡的分布式文件系统，目前主要用于中国气象局MICAPS系统。本文使用了https://github.com/CyanideCN/micaps_mdfs中mdfs的库，解析并绘制了空中实况观测，绘制了空中等高线、等温线、露点温度差、风向风速等。效果如图

50020200914080000.png 70020200914080000.png 85020200914080000.png

话不多说，直接上代码
'''

import sys
sys.path.append('/home/shortterm/micaps_mdfs/micaps_mdfs-master/')
import datetime
from mdfs import Station
import numpy as np
import pandas as pd
import os
import matplotlib.pyplot as plt
import matplotlib.colors as ccolor
import matplotlib.cm as cmx
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
import matplotlib.ticker as mticker
from cartopy.io.shapereader import Reader
from matplotlib.ticker import MaxNLocator
from matplotlib.colors import BoundaryNorm
SHP_china = '/home/shortterm/micaps_mdfs/china_shp'
SHP_world = '/home/shortterm/.local/share/cartopy/shapefiles/natural_earth/physical'
from scipy.interpolate import griddata#引入插值函数
from scipy.interpolate import Rbf#引入径向基函数
import matplotlib.patches as patches

# plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
# plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
def griddata_data(var,lat,lon,extent):
    data=pd.DataFrame([lon,lat,var]).T
    data=data.dropna()
    lon=data.iloc[:,0].values.reshape(-1,1)
    lat=data.iloc[:,1].values.reshape(-1,1)
    var=data.iloc[:,2].values.reshape(-1,1)

    points = np.concatenate([lon,lat],axis = 1)
    lon_min = extent[0]
    lon_max = extent[1]
    lat_min = extent[2]
    lat_max = extent[3]

    det_grid=0.5
    lon_grid, lat_grid = np.meshgrid(np.arange(lon_min,lon_max+det_grid,det_grid), 
                        np.arange(lat_min,lat_max+det_grid,det_grid))

    grid_data = griddata(points,var,(lon_grid,lat_grid),method = 'cubic')
    grid_data = grid_data[:,:,0]


    if lat_grid[0,0]<lat_grid[1,0]:
        lat_grid = lat_grid[-1::-1]
        grid_data = grid_data[-1::-1]
    return grid_data,lat_grid,lon_grid

def rbf_data(var,lat,lon,extent):
    data=pd.DataFrame([lon,lat,var]).T
    data=data.dropna()
    data=data[data.iloc[:,1]>5]
    lon=data.iloc[:,0].values.flatten()
    lat=data.iloc[:,1].values.flatten()
    var=data.iloc[:,2].values.flatten()
    olon=np.linspace(70,140,141)#设置网格经度
    olat=np.linspace(10,60,101)#设置网格纬度
    olon,olat=np.meshgrid(olon,olat)
    func=Rbf(lon,lat,var,function='linear')
    grid_data=func(olon,olat)
    return grid_data,olat,olon

def wsd2uv(ws, wd):     #风向风速转U、V
    wd = 270 - wd
    wd = wd /180 *np.pi
    x = ws * np.cos(wd)
    y = ws * np.sin(wd)
    return(x, y)

if __name__ == '__main__':

    Filetime=datetime.datetime.now()

    file_hour=''
    if (datetime.datetime.now().hour>=10) and (datetime.datetime.now().hour<22): 
        file_hour='080000'
    if (datetime.datetime.now().hour>=22) and (datetime.datetime.now().hour<24): 
        file_hour='200000'
    if (datetime.datetime.now().hour>=0) and (datetime.datetime.now().hour<10): 
        Filetime=Filetime-datetime.timedelta(days=1)
        file_hour='200000'
    FiletimeStr=Filetime.strftime('%Y%m%d')

    print(FiletimeStr)
    lev=['100','150','200','250','300','400','500','700','850','925','1000']
    extent=[70, 135, 15, 55]   #经纬度范围

    for i in lev:
#        file_name_dir='/mnt/z/highnew/plot/'+i+'/'+NowTimeStr+file_hour+'.000'
        file_name_dir=r'/mnt/micaps/highnew/plot/'+i+'/'+FiletimeStr+file_hour+'.000'
        print(file_name_dir)


    #    file_data = Station('./20200907080000.000')
        file_data =Station(file_name_dir)
        lon = file_data.data['Lon'] 
        lat = file_data.data['Lat']

        #  3   测站高度
        #  421 高度
        #  803 温度露点差
        #  601 温度
        #  201 风向
        #  203 风速

        hgt = file_data.data[421] 
        temp=  file_data.data[601]
        wind_dir=file_data.data[201]
        wind_speed=file_data.data[203]
        dewt=file_data.data[803]

        #将风向风速转变为UV
        U=[]
        V=[]
        for k,j in zip(wind_speed,wind_dir):
            u,v=wsd2uv(k,j)
            U.append(u)
            V.append(v)

        grid_hgt,lat_grid,lon_grid=rbf_data(hgt,lat,lon,extent)
        grid_temp,lat_grid,lon_grid=rbf_data(temp,lat,lon,extent)
        grid_dewt,lat_grid,lon_grid=rbf_data(dewt,lat,lon,extent)
        grid_dewt[grid_dewt<0]=0.0 # 因为插值的原因，grid_dewt部分值为小于0，将小于0的全部设置为0
        grid_U,lat_grid,lon_grid=rbf_data(U,lat,lon,extent)
        grid_V,lat_grid,lon_grid=rbf_data(V,lat,lon,extent)

        #proj= ccrs.LambertConformal(central_longitude=110, central_latitude=35) 
        proj= ccrs.PlateCarree() 
        fig = plt.figure(figsize=(10,8),dpi=300)  
        ax = fig.subplots(1, 1, subplot_kw={'projection': proj})

        #kedu 
        ax.xlabel_style={'size':33.5}
        ax.ylabel_style={'size':33.5}
        ax.set_extent(extent, ccrs.PlateCarree())
#        ax.coastlines(resolution='10m', linewidth=0.3)
        ax.add_geometries(Reader(os.path.join(SHP_china, 'cnhimap.shp')).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k', linewidth=0.3)
        ax.add_geometries(Reader(os.path.join(SHP_world, '10m_coastline.shp')).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k', linewidth=0.3)
        ax.set_xticks([70,75,80,85,90,95,100,105,110,115,120,125,130,135,140])#需要显示的经度，一般可用np.arange
        ax.set_yticks([10,15,20,25,30,35,40,45,50,55,60])#需要显示的纬度
        ax.xaxis.set_major_formatter(LongitudeFormatter())#将横坐标转换为经度格式
        ax.yaxis.set_major_formatter(LatitudeFormatter())#将纵坐标转换为纬度格式
        ax.tick_params(axis='both',labelsize=8,direction='in',length=2.75,width=0.55,right=True,top=True)#修改刻度样式
        ax.grid(linewidth=0.4, color='k', alpha=0.45, linestyle='--')#开启网格线

        #画等高线
        levels=np.arange(pd.Series(hgt).dropna().round(0).min(),pd.Series(hgt).dropna().round(0).max(),4)
        hgt_line = ax.contour(lon_grid,lat_grid,grid_hgt,levels=levels,colors=['blue'],linewidths=1.0,transform=ccrs.PlateCarree())
    #    levelsf = MaxNLocator(nbins=18).tick_values(564,600)
        #cf = ax.contourf(lon_grid,lat_grid,grid_data,levels=levelsf,transform=ccrs.PlateCarree())
        ax.clabel(
                hgt_line,  # Typically best results when labelling line contours.
                colors=['black'],
                manual=False,  # Automatic placement vs manual placement.
                inline=True,  # Cut the line where the label will be placed.
                fmt=' {:.0f} '.format,  # Labes as integers, with some extra space.
                fontsize=4.5,
                inline_spacing=0,
                )

        #画等温线
        levels=np.arange(pd.Series(temp).dropna().round(0).min(),pd.Series(temp).dropna().round(0).max(),4)
        temp_line = ax.contour(lon_grid,lat_grid,grid_temp,levels=levels,colors=['red'],linestyles='--',linewidths=1.0,transform=ccrs.PlateCarree())
        ax.clabel(
                temp_line,  # Typically best results when labelling line contours.
                colors=['black'],
                manual=False,  # Automatic placement vs manual placement.
                inline=True,  # Cut the line where the label will be placed.
                fmt=' {:.0f} '.format,  # Labes as integers, with some extra space.
                fontsize=4.5,
                inline_spacing=0,
                )
        #画温度露点差
        levels=np.arange(pd.Series(dewt).dropna().round(0).min(),pd.Series(dewt).dropna().round(0).max(),4)
        dewt_fill = ax.contourf(lon_grid,lat_grid,grid_dewt,levels=levels,cmap='Blues_r',transform=ccrs.PlateCarree())
        position=fig.add_axes([0.92,0.16,0.02,0.68])
        cb=fig.colorbar(dewt_fill,cax=position,extend='both',shrink=0.3,pad=0.01)

        #画风场
        barbs_inspace=3
        ax.barbs(lon_grid[::barbs_inspace,::barbs_inspace],lat_grid[::barbs_inspace,::barbs_inspace],
        grid_U[::barbs_inspace,::barbs_inspace],grid_V[::barbs_inspace,::barbs_inspace],linewidth=0.35,
        barb_increments={'half':2,'full':4,'flag':20},
        sizes=dict(spacing=0.1,width=0.05,emptybarb=0.02,height=0.35),length=3.5,zorder=5)
#        ax.text(0.02, 1.02, i,horizontalalignment='center',verticalalignment='center',transform=ax.transAxes)
        ax.text(0.09, 1.01, FiletimeStr+file_hour,horizontalalignment='center',verticalalignment='center',transform=ax.transAxes)
        ax.text(0.9, 1.01, i+'hPa,H,T,Wind,Dewt',horizontalalignment='center',verticalalignment='center',transform=ax.transAxes)


        plt.savefig('/home/shortterm/micaps_mdfs/sounding/'+str(i)+FiletimeStr+file_hour+'.png') 

'''