2020-07-13

#coding=utf-8

import numpyas np

import pandasas pd

import os

import re

import datetimeas dt

from matplotlib.datesimport date2num, DateFormatter

import matplotlib.pyplotas plt

import matplotlib.datesas mdates

from matplotlib.font_managerimport FontProperties

font = FontProperties(fname=r"C:\\Windows\\Fonts\\simsun.ttc", size=25)

from mpl_toolkits.axes_grid1import make_axes_locatable

import sys

sys.path.insert(1,r'H:\microwave_radiometer/')

import variables_dictionaryas vd

# vardic=vd.vardic

def create_empty_mother_array_for_dates_and_data(varname):

dates_all=np.array([],dtype=int).reshape(0, 6)#### year,mon,day,hour,min,sec

    if vd.vardic[varname]['levels']!=None:

levels=vd.vardic[varname]['levels']

data_all=np.array([]).reshape(0,len(levels))

# print levels

    if vd.vardic[varname]['levels']==None:

data_all=np.array([]).reshape(0,1)

return dates_all,data_all

print ""

def draw_2D(varname,times,data_all):

levels=vd.vardic[varname]['levels']

xx,yy=np.meshgrid(times,levels)

# print data_all.shape

    fig=plt.figure(figsize=(10,6))

rect=0.1, 0.1, 0.7, 0.8

    ax=fig.add_axes(rect)### rect=l, b, w, h

    fig.autofmt_xdate()

pcolormesh=ax.pcolormesh(xx,yy,data_all.transpose());

plt.title(unicode(varname,'utf-8')+'(%s)'%unicode(vd.vardic[varname]['unit'],'utf-8'),fontproperties=font)

xfmt = mdates.DateFormatter('%d %H:%M')

ax.xaxis.set_major_formatter(xfmt)

ax.set_xlim(dt.datetime(2020,4,16,20,0,0),dt.datetime(2020,4,17,20,0,0))

#### colorbar

    ax_colorbar=fig.add_axes([0.85,0.1,0.05,0.8])

# ax_colorbar=fig.add_axes(rect=[0.85,0.1,0.1,0.8])

    fig.colorbar(mappable=pcolormesh,cax=ax_colorbar)

plt.savefig(unicode(varname,'utf-8')+'.png')

plt.close()

def draw_1D(varname,times,data_all,yaxis_range=None):

fig=plt.figure(figsize=(10,6))

rect=0.1, 0.1, 0.7, 0.8

    ax=fig.add_axes(rect)

fig.autofmt_xdate()

# print times.shape,data_all.shape

# ax.scatter(times,data_all.transpose().flatten())

    ax.plot(times,data_all.transpose().flatten())

#plt.title(unicode(varname,'utf-8')+'(%s)'%unicode(vd.vardic[varname]['unit'],'utf-8'),fontproperties=font)

    plt.title(varname,fontproperties=font)

xfmt = mdates.DateFormatter('%d %H:%M')

ax.xaxis.set_major_formatter(xfmt)

ax.set_xlim(dt.datetime(2020,4,16,20,0,0),dt.datetime(2020,4,17,20,0,0))

if yaxis_range==None:pass

    if yaxis_range!=None:ax.set_ylim(yaxis_range[0],yaxis_range[-1])

#plt.savefig(unicode(varname,'utf-8')+'.png')

    plt.savefig(varname+'.png')

# plt.show()

    plt.close()

# dates_list=[]

# data=pd.read_csv(r'H:\microwave_radiometer\2020_04_16_case_study\\20200415_20200417.txt',skiprows=0,header=1,sep=' ')

# dates=np.array(data[['Year','Mon','Day','Hour']],dtype=int)

# rain=np.array(data[['PRE_1h']],dtype=np.float32)

# for year,mon,day,hour in dates:

#    # print year,mon,day,hour

#    dates_list.append(dt.datetime(year,mon,day,hour)+dt.timedelta(hours=+8))

# print dates_list

# draw_1D(unicode('PRE_1h'),dates_list,rain)

def read_MZ_minute_txt(filepath):

yearmonday_of_file=filepath[-12:-4]

# print yearmondate_of_file

    year,mon,day=yearmonday_of_file[0:4],yearmonday_of_file[4:6],yearmonday_of_file[6:8];year,mon,day=map(int,(year,mon,day))

# print year,mon,day

    rain_list=[]

time_list=[]

ff=open(filepath,'r')

ff.readline()

for iin range(1440):

line=ff.readline()

rain_list.append(float(line[56:60])/10)

# print dt.datetime(year,mon,day,hour=20,minute=1)+dt.timedelta(minutes=i)+dt.timedelta(days=-1)

        time_list.append(dt.datetime(year,mon,day,hour=20,minute=1)+dt.timedelta(minutes=i)+dt.timedelta(days=-1))

return rain_list,time_list

# rain_list1,time_list1=read_MZ_minute_txt(r'E:\microwave_radiometer\2020_04_16_case_study\AWS_M_Z_54237_20200416.txt')

# rain_list2,time_list2=read_MZ_minute_txt(r'E:\microwave_radiometer\2020_04_16_case_study\AWS_M_Z_54237_20200417.txt')

#

# rain_list=rain_list1+rain_list2

# time_list=time_list1+time_list2

# draw_1D(unicode('PRE_1min'),time_list,np.array(rain_list),yaxis_range=[0,1])

def get_data_by_varname_of_micro_and_daterange(varname='水汽廓线',startday='20200501',endday='20200509'):#####

    dates_all,data_all=create_empty_mother_array_for_dates_and_data(varname=varname)

year,mon,day=map(int,(startday[:4],startday[4:6],startday[6:]));startday_dt=dt.datetime(year,mon,day)

year,mon,day=map(int,(endday[:4],endday[4:6],endday[6:]));endday_dt=dt.datetime(year,mon,day)

while(startday_dt!=endday_dt+dt.timedelta(days=1)):

result=vd.get_dates_and_data(date=startday_dt.strftime('%Y%m%d'),varname=varname)

if result=="This date has no data":

print "%s has no data of %s"%(startday_dt,varname)

startday_dt=startday_dt+dt.timedelta(days=1)

continue

        print startday_dt

dates_array,data_array=result

dates_all=np.concatenate([dates_all,dates_array])

data_all=np.concatenate([data_all,data_array])

startday_dt=startday_dt+dt.timedelta(days=1)

return dates_all,data_all##### dates_all[0]=(year,mon,day,hour,min,sec)

# dates_all,data_all=get_data_by_varname_of_micro_and_daterange(varname='水汽廓线',startday='20200501',endday='20200510')

# print dates_all

# exit()

hourlyrainfaile=open(r'H:\microwave_radiometer\something_about_data_usage\2019june-2020june_hourly_rain.txt')

dates_rain=[];data_rain=[]

hourlyrainfaile.readline()

while True:

line=hourlyrainfaile.readline()

if line=='':

break

    ss=line.split(',')

# station=ss[5];year=ss[]

    year=ss[10];mon=ss[11];day=ss[12];hour=ss[13];rain_1h=ss[-2]

year,mon,day,hour=map(int,(year,mon,day,hour));rain_1h=float(rain_1h)

# print ss

# print "fasdfasdfasdf"

    dates_rain.append(dt.datetime(year,mon,day,hour))

data_rain.append(rain_1h)

data_rain=np.array(data_rain)

print data_rain.shape

print np.corrcoef(data_rain,data_rain)

exit()

font = FontProperties(fname=r"C:\\Windows\\Fonts\\simsun.ttc", size=45)

def draw_vardata_1D_and_rain_2yaxis(varname,vardata,vardata_times,raindata,rain_times,xaxis=['201907110000','202006100000']):#xaxis=['202004162130','202004162300']

    print vardata.shape,len(raindata)

# print raindata

    fig, ax1 = plt.subplots(figsize=(15,9),tight_layout=True)

ax2=ax1.twinx()

# a1=ax1.scatter(rain_times, raindata,s=1, c='blue') #绘制折线图像1,圆形点，标签，线宽

# a2=ax2.scatter(vardata_times, vardata,s=1, c='red') #同上

# plt.legend((a1,a2),('rain',unicode(varname,'utf-8')),scatterpoints=10,loc=1,bbox_to_anchor=(1,1), bbox_transform=ax1.transAxes,prop=font)

    a1=ax1.plot(rain_times, raindata, 'b-',label='rain',alpha=0.5)#绘制折线图像1,圆形点，标签，线宽

    a2=ax2.plot(vardata_times, vardata, 'r.',label=unicode(varname,'utf-8'),markersize=3,alpha=0.5)#同上

    labs = [l.get_label()for lin a1+a2]

ax1.legend(a1+a2,labs,loc=1,bbox_to_anchor=(1,1), bbox_transform=ax1.transAxes,prop=font)

fig.autofmt_xdate()

# plt.xlim(xaxis[0],xaxis[1])

    if xaxis!=None:plt.xlim(dt.datetime.strptime(xaxis[0],'%Y%m%d%H%M'),dt.datetime.strptime(xaxis[1],'%Y%m%d%H%M'))

ax1.tick_params(axis='x',labelsize=18,length=10,width=5,direction='out')

ax1.tick_params(axis='y',labelsize=30,colors='blue')

ax2.tick_params(axis='y',labelsize=30,colors='red')

ax1.set_yticks([2, 4, 8,20,50], minor=True);

ax1.grid(True,axis='x');ax1.grid(True, which='minor',color='blue',linestyle='--',linewidth=1)

# plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%H:%M'))

    plt.savefig(unicode(varname,'utf-8')+'_plus_rain'+'.png')

# plt.savefig(unicode(varname,'utf-8')+'_plus_rain'+'.eps')

# plt.savefig('water_plus_rain'+'.eps')

    plt.close()

def draw_vardata_2D_and_rain_2yaxis(varname,vardata,vardata_times,raindata,rain_times,xaxis=['201907110000','202006100000']):

# fig=plt.figure(figsize=(10,6))

    fig, ax1 = plt.subplots(figsize=(15,9))

ax2=ax1.twinx()

levels=vd.vardic[varname]['levels'];

# if levels[-1]==10000:amplification_factor=20000

    if levels[-1]==10000:amplification_factor=10000.0/np.amax(raindata)

# if levels[-1]==2000:amplification_factor=20000/5

    if levels[-1]==2000:amplification_factor=2000/np.amax(raindata)

# levels=np.array(levels)/1000.0

# print levels

# exit()

    xx,yy=np.meshgrid(vardata_times,levels)

font = FontProperties(fname=r"C:\\Windows\\Fonts\\simsun.ttc", size=45)

# plt.legend((a1,a2),('rain',unicode(varname,'utf-8')),prop=font)

    vardata = np.ma.masked_where(np.isnan(vardata),vardata)

a2=ax1.pcolormesh(xx,yy, vardata.transpose(),label=unicode(varname,'utf-8'),cmap='coolwarm')#同上

# a2=ax2.contourf(xx,yy, vardata.transpose()) #同上

    a1=ax2.plot(rain_times, np.array(raindata)*amplification_factor, c='k',label='rain',linewidth=1,alpha=0.5)#绘制折线图像1,圆形点，标签，线宽

# ax2.set_ylim(0,np.amax(raindata)*amplification_factor)

# a1=ax1.scatter(rain_times, raindata,s=10, c='blue',label='rain') #绘制折线图像1,圆形点，标签，线宽

    ax2.get_yaxis().set_visible(False)

ax2.invert_yaxis()

fig.legend(loc=1,bbox_to_anchor=(1,1), bbox_transform=ax1.transAxes,prop=font)

fig.autofmt_xdate()

plt.xlim(dt.datetime.strptime(xaxis[0],'%Y%m%d%H%M'),

            dt.datetime.strptime(xaxis[1],'%Y%m%d%H%M'))

ax1.tick_params(axis='x',labelsize=20,length=10,width=5,labelrotation=90)

ax1.tick_params(axis='y',labelsize=30,colors='k')

ax2.tick_params(axis='y',labelsize=20,colors='blue')

ax1.set_yticks([2, 4, 8,20,50], minor=True);

# ax1.grid(True,axis='x')

# plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%H:%M'))

    plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%y-%m-%d'))

# divider = make_axes_locatable(ax1)

# cax = divider.append_axes('bottom', size='5%', pad=0.1)

    plt.colorbar(mappable=a2)

plt.tight_layout()

plt.savefig(unicode(varname,'utf-8')+'_plus_rain'+'.png')

plt.close()

for varnamein vd.vardic:

# for varname in ['水汽总含量']:

# for varname in ['AI','KI','TTI','LI','SI','LCL','LFC','EC','BLH ']:

# if vd.vardic[varname]['filemark']=='L36000AMT':continue

# if vd.vardic[varname]['levels']==None:continue

# if '廓线' not in varname:continue        #### 画廓线

    if varnamein ['K通道接收机电压1B', 'V通道接收机电压1B', 'K通道模块噪声源电压1B', 'V通道模块噪声源电压1B', '方位1B', '俯仰1B','V通道亮温1D',

              'K通道接收机电压1D', 'V通道接收机电压1D', 'K通道模块噪声源电压1D', 'V通道模块噪声源电压1D', 'K通道模块噪声源温度', 'V通道模块噪声源温度',

              'K通道腔体温度', 'V通道腔体温度', 'K通道模块温度', 'V通道模块温度', '方位1D', '俯仰1D', '方位2B', '俯仰2B', '方位2D', '俯仰2D',

              '方位3D', '俯仰3D']:continue

    # varname

    if vd.vardic[varname]['levels']==None:

#### 所有值均为90.0;  多个分号

        print varname

# continue

        dates_all,data_all=get_data_by_varname_of_micro_and_daterange(varname=varname,startday='20190701',endday='20200610')

times=np.array([dt.datetime(i[0],i[1],i[2],i[3],i[4],i[5])for iin dates_all])#### xaxis

# draw_vardata_2D_and_rain_2yaxis(varname=varname,vardata=data_all,vardata_times=times,raindata=rain_list,rain_times=time_list)

        draw_vardata_1D_and_rain_2yaxis(varname=varname,vardata=data_all,vardata_times=times,raindata=data_rain,rain_times=dates_rain)

# print

# if vd.vardic[varname]['levels'] != None and '廓线' in varname:

#    ##### 所有值均为90.0;  多个分号

#    dates_all, data_all = get_data_by_varname_of_micro_and_daterange(varname=varname, startday='20190701',endday='20200610')

#    # dates_all, data_all = get_data_by_varname_of_micro_and_daterange(varname=varname, startday='20190701',endday='20191001')

#

#    times = np.array([dt.datetime(i[0], i[1], i[2], i[3], i[4], i[5]) for i in dates_all])  #### xaxis

#    # draw_vardata_2D_and_rain_2yaxis(varname=varname,vardata=data_all,vardata_times=times,raindata=rain_list,rain_times=time_list)

#    draw_vardata_2D_and_rain_2yaxis(varname=varname, vardata=data_all, vardata_times=times, raindata=data_rain,

#                                    rain_times=dates_rain)

# exit()

# if vd.vardic[varname]['levels']==None:

#    draw_1D(varname,times,data_all)

# if vd.vardic[varname]['levels']!=None:

#    draw_2D(varname,times,data_all)

©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 219,366评论 6赞 508
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 93,521评论 3赞 395
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 165,689评论 0赞 356
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 58,925评论 1赞 295
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 67,942评论 6赞 392
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 51,727评论 1赞 305
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 40,447评论 3赞 420
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 39,349评论 0赞 276
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 45,820评论 1赞 317
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,990评论 3赞 337
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 40,127评论 1赞 351
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 35,812评论 5赞 346
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 41,471评论 3赞 331
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 32,017评论 0赞 22
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 33,142评论 1赞 272
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 48,388评论 3赞 373
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 45,066评论 2赞 355