"""pyrad.graph.plot_timeseries===========================Functions to plot Pyrad datasets.. autosummary:: :toctree: generated/ plot_timeseries plot_timeseries_comp plot_monitoring_ts plot_intercomp_scores_ts plot_ml_ts plot_sun_retrieval_ts"""importpyartimportmatplotlib.pyplotaspltimportmatplotlib.datesasmdatesfromwarningsimportwarnimportnumpyasnpimportdatetimeimportmatplotlibasmplmpl.use("Agg")# Increase a bit font sizempl.rcParams.update({"font.size":16})mpl.rcParams.update({"font.family":"sans-serif"})
[docs]defplot_timeseries(tvec,data_list,fname_list,labelx="Time [UTC]",labely="Value",labels=["Sensor"],title="Time Series",period=0,timeformat=None,colors=None,linestyles=None,markers=None,ymin=None,ymax=None,dpi=72,):""" plots a time series Parameters ---------- tvec : datetime object time of the time series data_list : list of float array values of the time series fname_list : list of str list of names of the files where to store the plot labelx : str The label of the X axis labely : str The label of the Y axis labels : array of str The label of the legend title : str The figure title period : float measurement period in seconds used to compute accumulation. If 0 no accumulation is computed timeformat : str Specifies the tvec and time format on the x axis colors : array of str Specifies the colors of each line linestyles : array of str Specifies the line style of each line markers: array of str Specify the markers to be used for each line ymin, ymax: float Lower/Upper limit of y axis dpi : int dots per inch Returns ------- fname_list : list of str list of names of the created plots History -------- 201?.??.?? -fvj- creation 2017.08.21 -jgr- modified margins and grid + minor graphical updates 2018.03.05 -jgr- added x-limit of x axis to avoid unwanted error messages """tvec=np.array(tvec)# convert from pandas if neededifperiod>0:fori,datainenumerate(data_list):data*=period/3600.0data_list[i]=np.ma.cumsum(data)fig,ax=plt.subplots(figsize=[10,6],dpi=dpi)lab=Nonecol=Nonelstyle="--"marker="o"fori,datainenumerate(data_list):iflabelsisnotNone:lab=labels[i]ifcolorsisnotNone:col=colors[i]iflinestylesisnotNone:lstyle=linestyles[i]ifmarkersisnotNone:marker=markers[i]ax.plot(tvec,data,label=lab,color=col,linestyle=lstyle,marker=marker)ax.set_title(title)ax.set_xlabel(labelx)ax.set_ylabel(labely)ax.set_ylim(bottom=ymin,top=ymax)ax.set_xlim([tvec[0],tvec[-1]])# Turn on the gridax.grid()iftimeformatisnotNone:ax.xaxis.set_major_formatter(mdates.DateFormatter(timeformat))# rotates and right aligns the x labels, and moves the bottom of the# axes up to make room for themfig.autofmt_xdate()# Make a tight layoutfig.tight_layout()forfnameinfname_list:fig.savefig(fname,dpi=dpi)plt.close(fig)returnfname_list
[docs]defplot_timeseries_comp(date1,value1,date2,value2,fname_list,labelx="Time [UTC]",labely="Value",label1="Sensor 1",label2="Sensor 2",titl="Time Series Comparison",period1=0,period2=0,ymin=None,ymax=None,dpi=72,):""" plots 2 time series in the same graph Parameters ---------- date1 : datetime object time of the first time series value1 : float array values of the first time series date2 : datetime object time of the second time series value2 : float array values of the second time series fname_list : list of str list of names of the files where to store the plot labelx : str The label of the X axis labely : str The label of the Y axis label1, label2 : str legend label for each time series titl : str The figure title period1, period2 : float measurement period in seconds used to compute accumulation. If 0 no accumulation is computed dpi : int dots per inch ymin, ymax : float The limits of the Y-axis. None will keep the default limit. Returns ------- fname_list : list of str list of names of the created plots History -------- 201?.??.?? -fvj- created 2017.08.21 -jgr- changed some graphical aspects """if(period1>0)and(period2>0):# TODO: document this and check (sometimes artefacts)value1*=period1/3600.0value1=np.ma.cumsum(value1)value2*=period2/3600.0value2=np.ma.cumsum(value2)fig,ax=plt.subplots(figsize=[10,6.5],dpi=dpi)ax.plot(date1,value1,"b",label=label1,linestyle="--",marker="o")ax.plot(date2,value2,"r",label=label2,linestyle="--",marker="s")ax.legend(loc="best")ax.set_xlabel(labelx)ax.set_ylabel(labely)ax.set_title(titl)ax.grid()ax.set_ylim(bottom=ymin,top=ymax)ax.set_xlim([date2[0],date2[-1]])# rotates and right aligns the x labels, and moves the bottom of the# axes up to make room for themfig.autofmt_xdate()# Make a tight layoutfig.tight_layout()forfnameinfname_list:fig.savefig(fname,dpi=dpi)plt.close(fig)returnfname_list
[docs]defplot_monitoring_ts(date,np_t,cquant,lquant,hquant,field_name,fname_list,ref_value=None,vmin=None,vmax=None,np_min=0,labelx="Time [UTC]",labely="Value",titl="Time Series",dpi=72,plot_until_year_end=False,):""" plots a time series of monitoring data Parameters ---------- date : datetime object time of the time series np_t : int array number of points cquant, lquant, hquant : float array values of the central, low and high quantiles field_name : str name of the field fname_list : list of str list of names of the files where to store the plot ref_value : float the reference value vmin, vmax : float The limits of the y axis np_min : int minimum number of points to consider the sample plotable labelx : str The label of the X axis labely : str The label of the Y axis titl : str The figure title dpi : int dots per inch plot_until_year_end : bool if true will always set the maximum xtick to the end of the year Returns ------- fname_list : list of str list of names of the created plots """vmin_pyart,vmax_pyart=pyart.config.get_field_limits(field_name)ifvminisNone:vmin=vmin_pyartifvmaxisNone:vmax=vmax_pyart# plot only valid data (but keep first and last date)date2=np.array(date)isvalid=np.logical_not(np.ma.getmaskarray(cquant))ifnp_min>0:has_np=np_t>np_minisvalid=np.logical_and(isvalid,has_np)cquant_plt=cquant[isvalid]lquant_plt=lquant[isvalid]hquant_plt=hquant[isvalid]date_plt=date2[isvalid]ifnotisvalid[0]:cquant_plt=np.ma.append(np.ma.masked,cquant_plt)lquant_plt=np.ma.append(np.ma.masked,lquant_plt)hquant_plt=np.ma.append(np.ma.masked,hquant_plt)date_plt=np.ma.append(date2[0],date_plt)ifnotisvalid[-1]:cquant_plt=np.ma.append(cquant_plt,np.ma.masked)lquant_plt=np.ma.append(lquant_plt,np.ma.masked)hquant_plt=np.ma.append(hquant_plt,np.ma.masked)date_plt=np.ma.append(date_plt,date2[-1])fig=plt.figure(figsize=[15,13],dpi=dpi)ax=fig.add_subplot(2,1,1)ax.plot(date_plt,cquant_plt,"x-")ax.plot(date_plt,lquant_plt,"rx-")ax.plot(date_plt,hquant_plt,"rx-")ifref_valueisnotNone:ax.plot(date_plt,np.zeros(len(date_plt))+ref_value,"k--")ax.set_ylabel(labely)ax.set_title(titl)ax.set_ylim([vmin,vmax])ifplot_until_year_end:t0=date_plt[0]tend=datetime.datetime(year=t0.year,month=12,day=31,hour=23,minute=59)ax.set_xlim([t0,tend])else:# tight x axisax.autoscale(enable=True,axis="x",tight=True)ax.grid(True)ax=fig.add_subplot(2,1,2)ax.plot(date,np_t,"x-")ifnp_minisnotNone:ax.plot(date,np.zeros(len(date))+np_min,"k--")ifplot_until_year_end:ax.set_xlim([t0,tend])else:# tight x axisax.autoscale(enable=True,axis="x",tight=True)ax.set_ylabel("Number of Samples")ax.set_xlabel(labelx)# rotates and right aligns the x labels, and moves the bottom of the# axes up to make room for themfig.autofmt_xdate()forfnameinfname_list:fig.savefig(fname,dpi=dpi)plt.close(fig)returnfname_list
[docs]defplot_intercomp_scores_ts(date_vec,np_vec,meanbias_vec,medianbias_vec,quant25bias_vec,quant75bias_vec,modebias_vec,corr_vec,slope_vec,intercep_vec,intercep_slope1_vec,fname_list,ref_value=0.0,np_min=0,corr_min=0.0,labelx="Time UTC",titl="RADAR001-RADAR002 intercomparison",dpi=72,):""" plots a time series of radar intercomparison scores Parameters ---------- date_vec : datetime object time of the time series np_vec : int array number of points meanbias_vec, medianbias_vec, modebias_vec : float array mean, median and mode bias quant25bias_vec, quant75bias_vec: 25th and 75th percentile of the bias corr_vec : float array correlation slope_vec, intercep_vec : float array slope and intercep of a linear regression intercep_slope1_vec : float the intercep point of a inear regression of slope 1 ref_value : float the reference value np_min : int The minimum number of points to consider the result valid corr_min : float The minimum correlation to consider the results valid labelx : str The label of the X axis titl : str The figure title Returns ------- fname_list : list of str list of names of the created plots """# plot only valid data (but keep first and last date)date2=np.array(date_vec)isvalid=np.logical_not(np.ma.getmaskarray(meanbias_vec))isvalid_corr=np.logical_not(np.ma.getmaskarray(corr_vec))ifnp_min>0:has_np=np_vec>np_minisvalid=np.logical_and(isvalid,has_np)ifcorr_min>0:has_corr_min=corr_vec>corr_minisvalid=np.logical_and(isvalid,has_corr_min)meanbias_plt=meanbias_vec[isvalid]medianbias_plt=medianbias_vec[isvalid]quant25bias_plt=quant25bias_vec[isvalid]quant75bias_plt=quant75bias_vec[isvalid]modebias_plt=modebias_vec[isvalid]intercep_plt=intercep_slope1_vec[isvalid]corr_plt=corr_vec[isvalid_corr]date_corr=date2[isvalid_corr]date_plt=date2[isvalid]ifnotisvalid[0]:meanbias_plt=np.ma.append(np.ma.masked,meanbias_plt)medianbias_plt=np.ma.append(np.ma.masked,medianbias_plt)quant25bias_plt=np.ma.append(np.ma.masked,quant25bias_plt)quant75bias_plt=np.ma.append(np.ma.masked,quant75bias_plt)modebias_plt=np.ma.append(np.ma.masked,modebias_plt)intercep_plt=np.ma.append(np.ma.masked,intercep_plt)date_plt=np.ma.append(date2[0],date_plt)ifnotisvalid[-1]:meanbias_plt=np.ma.append(meanbias_plt,np.ma.masked)medianbias_plt=np.ma.append(medianbias_plt,np.ma.masked)quant25bias_plt=np.ma.append(quant25bias_plt,np.ma.masked)quant75bias_plt=np.ma.append(quant75bias_plt,np.ma.masked)modebias_plt=np.ma.append(modebias_plt,np.ma.masked)intercep_plt=np.ma.append(intercep_plt,np.ma.masked)date_plt=np.ma.append(date_plt,date2[-1])ifnotisvalid_corr[0]:corr_plt=np.ma.append(np.ma.masked,corr_plt)date_corr=np.ma.append(date2[0],date_corr)ifnotisvalid_corr[-1]:corr_plt=np.ma.append(corr_plt,np.ma.masked)date_corr=np.ma.append(date_corr,date2[-1])fig=plt.figure(figsize=[10,20],dpi=dpi)ax=fig.add_subplot(4,1,1)ax.plot(date_plt,medianbias_plt,"bx-",label="median")ax.plot(date_plt,meanbias_plt,"rx-",label="mean")ax.plot(date_plt,modebias_plt,"gx-",label="mode")ax.plot(date_plt,intercep_plt,"yx-",label="intercep of slope 1 LR")ifref_valueisnotNone:ax.plot(date_plt,np.zeros(len(date_plt))+ref_value,"k--")# plt.legend(loc='best')ax.set_ylabel("bias [dB]")ax.set_title(titl)ax.set_ylim([-5.0,5.0])# tight x axisax.autoscale(enable=True,axis="x",tight=True)ax.grid(True)ax=fig.add_subplot(4,1,2)ax.plot(date_plt,medianbias_plt,"bx-",label="median")ax.plot(date_plt,quant25bias_plt,"rx-",label="25-percentile")ax.plot(date_plt,quant75bias_plt,"rx-",label="75-percentile")ifref_valueisnotNone:ax.plot(date_plt,np.zeros(len(date_plt))+ref_value,"k--")# plt.legend(loc='best')ax.set_ylabel("bias [dB]")ax.set_ylim([-5.0,5.0])# tight x axisax.autoscale(enable=True,axis="x",tight=True)ax.grid(True)ax=fig.add_subplot(4,1,3)ax.plot(date_corr,corr_plt,"bx-")ifcorr_min>0:ax.plot(date_corr,np.zeros(len(date_corr))+corr_min,"k--")ax.set_ylabel("correlation")ax.set_ylim([0.0,1.0])# tight x axisax.autoscale(enable=True,axis="x",tight=True)ax.grid(True)ax=fig.add_subplot(4,1,4)ax.plot(date2,np_vec,"bx-")ifnp_min>0:ax.plot(date2,np.zeros(len(date2))+np_min,"k--")ax.set_ylabel("Number of Samples")ax.set_xlabel(labelx)# tight x axisax.autoscale(enable=True,axis="x",tight=True)# rotates and right aligns the x labels, and moves the bottom of the# axes up to make room for themfig.autofmt_xdate()forfnameinfname_list:fig.savefig(fname,dpi=dpi)plt.close(fig)returnfname_list
[docs]defplot_ml_ts(dt_ml_arr,ml_top_avg_arr,ml_top_std_arr,thick_avg_arr,thick_std_arr,nrays_valid_arr,nrays_total_arr,fname_list,labelx="Time UTC",titl="Melting layer time series",dpi=72,):""" plots a time series of melting layer data Parameters ---------- dt_ml_arr : datetime object time of the time series np_vec : int array number of points meanbias_vec, medianbias_vec, modebias_vec : float array mean, median and mode bias quant25bias_vec, quant75bias_vec: 25th and 75th percentile of the bias corr_vec : float array correlation slope_vec, intercep_vec : float array slope and intercep of a linear regression intercep_slope1_vec : float the intercep point of a inear regression of slope 1 ref_value : float the reference value np_min : int The minimum number of points to consider the result valid corr_min : float The minimum correlation to consider the results valid labelx : str The label of the X axis titl : str The figure title Returns ------- fname_list : list of str list of names of the created plots """fig=plt.figure(figsize=[10,15],dpi=dpi)ax=fig.add_subplot(3,1,1)ax.plot(dt_ml_arr,ml_top_avg_arr,"bx-",label="avg")ax.plot(dt_ml_arr,ml_top_avg_arr+ml_top_std_arr,"rx-",label="avg+std")ax.plot(dt_ml_arr,ml_top_avg_arr-ml_top_std_arr,"rx-",label="avg-std")# plt.legend(loc='best')ax.set_ylabel("Top height [m MSL]")ax.set_title(titl)ax.set_ylim([0.0,6000.0])ax.set_xlim([dt_ml_arr[0],dt_ml_arr[-1]])# tight x axisax.autoscale(enable=True,axis="x",tight=True)ax.grid(True)ax=fig.add_subplot(3,1,2)ax.plot(dt_ml_arr,thick_avg_arr,"bx-",label="avg")ax.plot(dt_ml_arr,thick_avg_arr+thick_std_arr,"rx-",label="avg+std")ax.plot(dt_ml_arr,thick_avg_arr-thick_std_arr,"rx-",label="avg-std")# plt.legend(loc='best')ax.set_ylabel("Thickness [m]")ax.set_ylim([0.0,3000.0])ax.set_xlim([dt_ml_arr[0],dt_ml_arr[-1]])# tight x axisax.autoscale(enable=True,axis="x",tight=True)ax.grid(True)ax=fig.add_subplot(3,1,3)ax.plot(dt_ml_arr,nrays_valid_arr,"bx-",label="N valid rays")ax.plot(dt_ml_arr,nrays_total_arr,"rx-",label="rays total")# plt.legend(loc='best')ax.set_ylabel("Rays")ax.set_xlabel(labelx)ax.set_ylim([0,np.max(nrays_total_arr)+5])ax.set_xlim([dt_ml_arr[0],dt_ml_arr[-1]])# tight x axisax.autoscale(enable=True,axis="x",tight=True)ax.grid(True)# rotates and right aligns the x labels, and moves the bottom of the# axes up to make room for themfig.autofmt_xdate()forfnameinfname_list:fig.savefig(fname,dpi=dpi)plt.close(fig)returnfname_list
[docs]defplot_sun_retrieval_ts(sun_retrieval,data_type,fname_list,labelx="Date",titl="Sun retrieval Time Series",dpi=72,):""" plots sun retrieval time series series Parameters ---------- sun_retrieval : tuple tuple containing the retrieved parameters data_type : str parameter to be plotted fname_list : list of str list of names of the files where to store the plot labelx : str the x label titl : str the title of the plot dpi : int dots per inch Returns ------- fname_list : list of str list of names of the created plots """value_std=Noneref=Nonedate=sun_retrieval[1]ifdata_type=="nhits_h":value=sun_retrieval[2]labely="Number of sun hits H channel"vmin=0vmax=np.max(sun_retrieval[2])+1elifdata_type=="el_width_h":value=sun_retrieval[3]labely="Elevation beamwidth H channel (Deg)"vmin=0.0vmax=4.0elifdata_type=="az_width_h":value=sun_retrieval[4]labely="Azimuth beamwidth H channel (Deg)"vmin=0.0vmax=4.0elifdata_type=="el_bias_h":value=sun_retrieval[5]ref=np.zeros(len(value))labely="Elevation pointing bias H channel (Deg)"vmin=-2.0vmax=2.0elifdata_type=="az_bias_h":value=sun_retrieval[6]ref=np.zeros(len(value))labely="Azimuth pointing bias H channel (Deg)"vmin=-2.0vmax=2.0elifdata_type=="dBm_sun_est":value=sun_retrieval[7]value_std=sun_retrieval[8]labely="Sun Power H channel (dBm)"vmin=-110.0vmax=-90.0elifdata_type=="rx_bias_h":value=10.0*np.ma.log10(sun_retrieval[9])-10.0*np.ma.log10(sun_retrieval[21])value_std=sun_retrieval[8]ref=np.zeros(len(value))labely="Receiver bias H channel (dB)"vmin=-5.0vmax=5.0elifdata_type=="sf_h":value=10.0*np.ma.log10(sun_retrieval[9])# value_std = sun_retrieval[8]ref=10.0*np.ma.log10(sun_retrieval[21])labely="Observed solar flux H channel (dB(sfu))"vmin=15.0vmax=30.0elifdata_type=="nhits_v":value=sun_retrieval[10]labely="Number of sun hits V channel"vmin=0vmax=np.max(sun_retrieval[10])+1elifdata_type=="el_width_v":value=sun_retrieval[11]labely="Elevation beamwidth V channel (Deg)"vmin=0.0vmax=4.0elifdata_type=="az_width_v":value=sun_retrieval[12]labely="Azimuth beamwidth V channel (Deg)"vmin=0.0vmax=4.0elifdata_type=="el_bias_v":value=sun_retrieval[13]ref=np.zeros(len(value))labely="Elevation pointing bias V channel (Deg)"vmin=-2.0vmax=2.0elifdata_type=="az_bias_v":value=sun_retrieval[14]ref=np.zeros(len(value))labely="Azimuth pointing bias V channel (Deg)"vmin=-2.0vmax=2.0elifdata_type=="dBmv_sun_est":value=sun_retrieval[15]value_std=sun_retrieval[16]labely="Sun Power V channel (dBm)"vmin=-110.0vmax=-90.0elifdata_type=="rx_bias_v":value=10.0*np.ma.log10(sun_retrieval[17])-10.0*np.ma.log10(sun_retrieval[21])value_std=sun_retrieval[16]ref=np.zeros(len(value))labely="Receiver bias V channel (dB)"vmin=-5.0vmax=5.0elifdata_type=="sf_v":value=10.0*np.ma.log10(sun_retrieval[17])# value_std = sun_retrieval[16]ref=10.0*np.ma.log10(sun_retrieval[21])labely="Observed solar flux V channel (dB(sfu))"vmin=15.0vmax=30.0elifdata_type=="nhits_zdr":value=sun_retrieval[18]labely="Number of sun hits ZDR"vmin=0vmax=np.max(sun_retrieval[18])+1elifdata_type=="ZDR_sun_est":value=sun_retrieval[19]value_std=sun_retrieval[20]ref=np.zeros(len(value))labely="Sun ZDR (dB)"vmin=-2.0vmax=2.0mask=np.ma.getmaskarray(value)ifmask.all():warn("Unable to create figure "+" ".join(fname_list)+". No valid data")returnNone# plot only valid data (but keep first and last date)isvalid=np.logical_not(mask)date2=np.array(date)value_plt=value[isvalid]date_plt=date2[isvalid]ifnotisvalid[0]:value_plt=np.ma.append(np.ma.masked,value_plt)date_plt=np.ma.append(date2[0],date_plt)ifnotisvalid[-1]:value_plt=np.ma.append(value_plt,np.ma.masked)date_plt=np.ma.append(date_plt,date2[-1])fig,ax=plt.subplots(figsize=[10,6],dpi=dpi)ax.plot(date_plt,value_plt,"x-")ifvalue_stdisnotNone:value_std_plt=value_std[isvalid]ifnotisvalid[0]:value_std_plt=np.ma.append(np.ma.masked,value_std_plt)ifnotisvalid[-1]:value_std_plt=np.ma.append(value_std_plt,np.ma.masked)ax.plot(date_plt,value_plt+value_std_plt,"rx-")ax.plot(date_plt,value_plt-value_std_plt,"rx-")ifrefisnotNone:ref_plt=ref[isvalid]ifnotisvalid[0]:ref_plt=np.ma.append(ref[0],ref_plt)ifnotisvalid[-1]:ref_plt=np.ma.append(ref_plt,ref[-1])ax.plot(date_plt,ref_plt,"k--")ax.set_xlabel(labelx)ax.set_ylabel(labely)ax.set_title(titl)ax.set_ylim([vmin,vmax])ax.set_xlim([date_plt[0],date_plt[-1]])# tight x axisax.autoscale(enable=True,axis="x",tight=True)ax.grid(True)# rotates and right aligns the x labels, and moves the bottom of the# axes up to make room for themfig.autofmt_xdate()forfnameinfname_list:fig.savefig(fname,dpi=dpi)plt.close(fig)returnfname_list