#!/bin/env python#coding=utf-8import codecsimport numpy as npimport datetimeimport reimport pandas as pdimport cartopy.crs as ccrsimport cartopy.feature as cfeaturefrom metpy.calc import wind_componentsfrom metpy.plots import add_metpy_logo, current_weather, sky_cover, StationPlot, wx_code_mapfrom cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatterfrom metpy.units import unitsimport matplotlib.pyplot as pltfrom matplotlib.patches import Polygonimport cartopy.io.shapereader as shpreader
def func_return_self(x):    return x
def split(str_in, sep=None, convert_func=None):    ss = str_in.split(sep)    nn = len(ss)    rst = [None]*nn    if convert_func is None:        cf = [func_return_self] * nn    else:        cf = [_ for _ in convert_func]   # print(nn, len(cf))    for i in range(nn):        xcf = func_return_self if cf[i] is None else cf[i]        rst[i] =  xcf(ss[i])    return rst
def d1_parse(filename,encoding="gb2312"):    """    读micaps第1类数据文件到内存    :return:    """    with codecs.open(filename, mode='r', encoding=encoding) as fid:        contents = fid.readlines()    dataflag, style, title = split(contents[0],None,None)
    nn = len(contents)    nx = int((nn-2)/2)
    data = {        "stid":[None]*nx,        "lon":[None]*nx,        "lat":[None]*nx,        "height":[None]*nx,        "stlevel":[None]*nx,        "cloud":[None]*nx,        "wd":[None]*nx,        "ws":[None]*nx,        "psfc":[None]*nx,        "pd3":[None]*nx,        "pw1":[None]*nx,        "pw2":[None]*nx,        "prc6":[None]*nx,        "lc1":[None]*nx,        "lc2":[None]*nx,        "lc3":[None]*nx,        "dp":[None]*nx,        "njd":[None]*nx,        "cw":[None]*nx,        "tc":[None]*nx,        "mc":[None]*nx,        "hc":[None]*nx,        "flag1":[None]*nx,        "flag2":[None]*nx,        "dtc24":[None]*nx,        "dp24":[None]*nx    }    data = pd.DataFrame(data)
    j=0    for i in range(2,nn,2):        iline1 = contents[i]        iline2 = contents[i+1]        iline = " ".join([iline1,iline2])        #for iline  in contents[2:]:        data.stid[j], data.lon[j], data.lat[j], data.height[j], data.stlevel[j],data.cloud[j],data.wd[j],data.ws[j],data.psfc[j],data.pd3[j], data.pw1[j], data.pw2[j], \        data.prc6[j], data.lc1[j],data.lc2[j],data.lc3[j],       data.dp[j],   data.njd[j],  data.cw[j],   data.tc[j],   data.mc[j],    data.hc[j],   data.flag1[j], data.flag2[j], data.dtc24[j],data.dp24[j] \            =split(iline,None,        [int, float,float,float,int,float,float,float,float,float,int,int,\         float,float,float,float,float,float,int,  float,float, float,int, int,float,float]                   )        j+=1        #data.append([stid, lon, lat, height, stlevel,cloud,wd,ws,psfc,pd3, pw1, pw2,\        #prc6, lc1,lc2,lc3,       dp,   njd,  cw,   tc,   mc,    hc,   flag1, flag2, dtc24,dp24])    #print(data)    return dataflag, style, title,datadef d4_parse(filename,encoding="gb2312"):    """    读micaps第4类数据文件到内存    :return:    """    try:        file_object = codecs.open(filename, mode='r', encoding=encoding)        all_the_text = file_object.read().strip()        file_object.close()        contents = re.split('[\s]+', all_the_text)        if len(contents) < 23:            return        dataflag = contents[0].strip()        style = contents[1].strip()        title = contents[2].strip()        yy = int(contents[3].strip())        mm = int(contents[4].strip())        dd = int(contents[5].strip())        hh = int(contents[6].strip())        forehh = int(contents[7].strip())        level = contents[8].strip()        deltalon = float(contents[9].strip())        deltalat = float(contents[10].strip())        beginlon = float(contents[11].strip())        endlon = float(contents[12].strip())        beginlat = float(contents[13].strip())        endlat = float(contents[14].strip())        sumlon = int(contents[15].strip())        sumlat = int(contents[16].strip())        distance = float(contents[17].strip())        min = float(contents[18].strip())        max = float(contents[19].strip())        def1 = contents[20].strip()        def2 = contents[21].strip()        x = np.arange(beginlon, endlon + 0.9*deltalon, deltalon)        y = np.arange(beginlat, endlat + 0.9*deltalat, deltalat)        X, Y = np.meshgrid(x, y)        if dataflag == 'diamond' and style == '4':            begin = 22            Z = np.zeros((sumlat, sumlon))            for i in range(sumlon):                for j in range(sumlat):                    Z[j, i] = float(contents[begin + j * sumlon + i])        return dataflag, style,title,yy,mm,dd,hh,forehh,level,X,Y,Z,distance,min,max,def1,def2    except Exception as err:        print(u'【{0}】{1}-{2}'.format(filename, err, datetime.datetime.now()))
if __name__ == "__main__":    dataflag, style, title, yy, mm, dd, hh, forehh, level, X, Y, p0_data, distance, dmin, dmax, def1, def2 = d4_parse("/home/hliang/share/share/plotest/20190126/p0/19012602.000")    station_dataflag, staiton_style, station_title, station_data = d1_parse("/home/hliang/share/share/plotest/20190126/plot/19012108.000")    station_data = station_data.replace(9999,np.nan)    station_data = station_data[station_data.stlevel==2]#.dropna()    lat_max=50    lat_min=15    lon_max=130    lon_min=100    proj = ccrs.LambertCylindrical(central_longitude=(lon_max+lon_min)/2.0)#ccrs.LambertConformal(central_longitude=(lon_max+lon_min)/2.0, central_latitude=(lat_max+lat_min)/2.0,standard_parallels=[(lat_max+lat_min)/2.0])    fig = plt.figure(figsize=(10, 10))    ax = fig.add_subplot(1, 1, 1, projection=proj)    ax.set_extent([lon_min,lon_max,lat_min,lat_max],proj)    ax.set_xticks(range(lon_min,lon_max+2,10), crs=proj)    ax.set_yticks(range(lat_min,lat_max+2,10), crs=proj)    shp = shpreader.Reader('/home/hliang/share/share/plotest/shp_resources/CHN_adm/CHN_adm1.shp')    ax.add_geometries(shp.geometries(), crs=proj, facecolor="none", edgecolor="grey", lw=0.8)
    contours = ax.contour(X, Y, p0_data,levels=np.arange(dmin,dmax+0.9*distance,distance), colors="black",transform=proj,linewidths=0.5)    plt.clabel(contours, fontsize=8,fmt='%.1f')    stationplot = StationPlot(ax, station_data['lon'].values, station_data['lat'].values, clip_on=True,                              transform=proj, fontsize=6)    u,v =  wind_components(station_data['ws'].values * units('m/s'),station_data['wd'].values * units.degree)    stationplot.plot_parameter('NW', station_data['tc'], color='red')    stationplot.plot_barb(u, v,sizes=dict(emptybarb=0.15))
    cw = station_data.cw.replace(np.nan,0).astype(np.int)    pw1 = station_data.pw1.replace(np.nan,0).astype(np.int)    stationplot.plot_symbol('W', cw, current_weather,color='fuchsia')    stationplot.plot_symbol('SE', pw1, current_weather,color='y')
    '''    for record,state in zip(shp.records(), shp.geometries()):        name = record.attributes["NAME_1"]        if name=="Jiangsu":            ax.add_geometries(state,crs=proj,facecolor="none",edgecolor="grey",lw=0.5)    '''    plt.savefig('./fig1.png',dpi=300)    #plt.show()