# -*- coding: utf-8 -*- """ Grafica ceilómetro geoportal con estimaciones de Capa Límite Lauri :) """ from Funciones_Ceilometro import * from Funciones_Estadistica import * from Funciones_Radiometro import * from Funciones_RWP import * import traceback import matplotlib.font_manager as fm prop = fm.FontProperties(fname='/home/lherrera/Pagina/fonts/AvenirLTStd-Heavy.otf') prop_1 = fm.FontProperties(fname='/home/lherrera/Pagina/fonts/AvenirLTStd-Book.otf') prop_2 = fm.FontProperties(fname='/home/lherrera/Pagina/fonts/AvenirLTStd-Black.otf') def Richardson_comp(Alt, Theta_v, Zonal, Meridional): Ri = [] for i in range(len(Zonal)): cte = 9.8 / Theta_v[0] num = Theta_v[i] - Theta_v[0] den = (Zonal[i])**2 + (Meridional[i])**2 Ri.append(cte * (num/den) * Alt[i]) Ri = np.array(Ri) return Ri def lee_datos(ceilometro, Ventana): Fecha_1 = datetime.strftime (datetime.utcnow(), "%Y-%m-%d") Fecha_2 = datetime.strftime (datetime.utcnow() - timedelta(hours=Ventana), "%Y-%m-%d") # Genera la lista de fechas que se deben leer de acuerdo con la ventana ingresada File_List = pd.date_range(Fecha_2, Fecha_1, freq='1D') Backs = [] Fechas = [] for idd, Fecha in enumerate(File_List): year = datetime.strftime(Fecha, '%Y') month = datetime.strftime(Fecha, '%b') day = datetime.strftime(Fecha, '%d') fname_miel = '/mnt/ALMACENAMIENTO/ceilometro/datos/ceilometro'+str(ceilometro) \ +'/'+year+'/' + month+'/CEILOMETER_1_LEVEL_2_' +day+'.his' fname = fname_miel[-27:] try: BIN_fname = np.genfromtxt(fname_miel,delimiter=', ',dtype=object,usecols=(0,4),skip_header=2,\ converters = {0: lambda s: datetime.strptime(s, "%Y-%m-%d %H:%M:%S")}) DATA = np.array([decode_hex_string(BIN_fname[i,1]) for i in range(len(BIN_fname))]).T File_Dates = BIN_fname[:,0].tolist() except: DATA = np.ones([450,5397]) * np.NaN File_Dates = np.ones([5397]) * np.NaN # Concatena las matrices de los días seleccionados if idd == 0: Backs = DATA Fechas = File_Dates else: Backs = np.concatenate([Backs,DATA],axis=1) Fechas = np.concatenate([Fechas, File_Dates]) Backs = promedio(Backs, 3, 15) # Funcion que promedia en X minutos y Y metros Backs = Backs.astype(np.int) Backs[Backs < 100] = 100 Backs = pd.DataFrame(Backs.T, Fechas - timedelta(hours=5)) return Backs # Graficación por ventanas def Grafica_Ventanas(Backs, Ventana, NombrePNG, CLA=None, Nombre_CLA=None): Fecha_1 = datetime.strftime (datetime.now(), "%Y-%m-%d %H:%M") Fecha_2 = datetime.strftime (datetime.now() - timedelta(hours=Ventana), "%Y-%m-%d %H:%M") Backs = Backs_Ceil[Fecha_2:Fecha_1] plt.close ('all') # Fechas, ticks x_lims = (Backs.index[0], Backs.index[-1]) x_lims = mdates.date2num(x_lims) y_lims = [0, 4500] if Ventana < 7: minute =mdates.MinuteLocator(interval = 25) elif 7 < Ventana < 20: minute =mdates.MinuteLocator(interval = 90) elif 20 < Ventana < 48: minute =mdates.MinuteLocator(interval = 180) elif 48 <= Ventana: minute =mdates.MinuteLocator(interval = 400) # Cmap, colores orig_cmap = matplotlib.cm.jet shrunk_cmap = shiftedColorMap(orig_cmap, start=0.0, midpoint=0.65, stop=0.85, name='shrunk') fig = plt.figure(facecolor='w', edgecolor='w', figsize=[20,6]) plt.subplots_adjust(hspace=0.5) gs = gridspec.GridSpec(1, 2, width_ratios=[1.5, 3.5]) gs.update(wspace = 0.08, hspace = 0) ax = fig.add_subplot(gs[1]) ax.spines['bottom'].set_color((0.45, 0.45, 0.45)) ax.spines['top'].set_color((0.45, 0.45, 0.45)) ax.xaxis.label.set_color((0.45, 0.45, 0.45)) ax.tick_params(axis='x', colors=(0.45, 0.45, 0.45)) ax.spines['right'].set_color((0.45, 0.45, 0.45)) ax.spines['left'].set_color((0.45, 0.45, 0.45)) ax.yaxis.label.set_color((0.45, 0.45, 0.45)) ax.tick_params(axis='y', colors=(0.45, 0.45, 0.45)) imresobj = ax.imshow(Backs.T,cmap = shrunk_cmap, norm=LogNorm(), extent =[x_lims[0], x_lims[-1], y_lims[0], y_lims[1]],vmin=10, vmax=13000, aspect = 'auto', origin='lower') if CLA is not None: CLA = CLA[Fecha_2:Fecha_1] ax.scatter(CLA.index, CLA, color='r',s=33,lw=0.5,edgecolor='k', label=u'Altura estimada de la CLA') legend = ax.legend(bbox_to_anchor=(-0.01, -0.66, 0.3, .3), loc=9,ncol=1, mode="expand", borderaxespad=0., prop={'size':13}, shadow=False) frame = legend.get_frame() frame.set_edgecolor('0.90') for text in legend.get_texts(): plt.setp(text, color = (0.45, 0.45, 0.45)) ax.set_ylabel('Altura [km]', fontsize = 17, fontproperties=prop_1) plt.xticks(fontsize=17, fontproperties=prop_1) plt.yticks(fontsize=17, fontproperties=prop_1) plt.ylim([0,4500]) plt.xlim([Backs.index[0], Backs.index[-1]]) ax.annotate(' ', xy=(0, -0.11), xycoords='axes fraction', xytext=(1, -0.11), arrowprops=dict(arrowstyle="-", color=(.57, 0.45, 0.45))) ax.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M\n\n%d-%b')) ax.xaxis.set_major_locator(minute) ax.xaxis_date() plt.yticks(range(0,5000,500),['0.0','0.5','1.0','1.5','2.0','2.5','3.0','3.5','4.0','4.5'],fontsize=18, fontproperties=prop_1) cbar = plt.colorbar(imresobj,ticks=[10,10**2,10**3,10**4], orientation='horizontal', format="%.1f", cax=fig.add_axes([0.375,-0.13,0.525,0.05]),drawedges=False) cbar.set_label('Intensidad Backscatter $[10^{-9}m^{-1}sr^{-1}]$',size=17,fontproperties=prop_1, color=(0.45, 0.45, 0.45)) cbar.ax.set_xticklabels(['10','10$^{2}$', '10$^{3}$', '10$^{4}$'],fontsize=16, color = (0.45, 0.45, 0.45), fontproperties=prop_1) cbar.solids.set_edgecolor("face") cbar.outline.set_edgecolor('white') plt.savefig('test1.png', bbox_inches='tight',dpi=150) if Nombre_CLA is not None: os.system("scp test1.png lherrera@siata.gov.co:/var/www/laura/Pagina/Ceilometros/Ceilometro"+str(NombrePNG)+"_"+str(Nombre_CLA)+'_'+str(Ventana)+".png") else: os.system("scp test1.png lherrera@siata.gov.co:/var/www/laura/Pagina/Ceilometros/Ceilometro"+str(NombrePNG)+"_"+str(Ventana)+".png") os.system("rm test1.png") def CLA_Ceil(Backs_Ceil, freq_min): Backs = Backs_Ceil.resample(str(freq_min)+'Min', how='mean', label='right', closed='right') MLH_grad_filtro = [] MLH_var_filtro = [] for i in range(len(Backs)): try: Bz = mediamovil(Backs.values[i],9)[(9/2):] # Media movil 90m # Filtro Nubes Bz_filtro = cloudfilter_test(Bz, gradiente(Bz)) if np.max(Bz_filtro) > 1200: Bz_filtro = cloudfilter_test(Bz_filtro, gradiente(Bz_filtro)) if np.max(Bz_filtro) > 1200: Bz_filtro = cloudfilter_test(Bz_filtro, gradiente(Bz_filtro)) # Gradiente Bz_grad_filtro = gradiente(Bz_filtro) # Gradiente minimo Datos filtrados notnull = Bz_grad_filtro[0:200][np.isfinite(Bz_grad_filtro[:200])==True] if np.sum(Bz) > 450000: MLH_grad_filtro.append(np.NaN) else: MLH_grad_filtro.append((np.where((notnull) == min((notnull)))[0][0]*10)) # Varianza Bz_var_filtro = np.array(var_movil(Bz_filtro,20)) # Maxima Varianza Datos Filtrados if np.sum(Bz) > 450000: MLH_var_filtro.append(np.NaN) else: MLH_var_filtro.append((np.where(Bz_var_filtro[:200][np.isfinite(Bz_var_filtro[:200])==True] == max(Bz_var_filtro[:200][np.isfinite(Bz_var_filtro[:200])==True]))[0][0]*10)+100) except: MLH_var_filtro.append(np.NaN) MLH_grad_filtro.append(np.NaN) MLH_grad_filtro = np.array(MLH_grad_filtro).astype(float) MLH_grad_filtro[MLH_grad_filtro == 0]= np.NaN MLH_var_filtro = np.array(MLH_var_filtro).astype(float) MLH_var_filtro[MLH_var_filtro == 0]= np.NaN # Filtro Vecinos Filtro_Vecinos_var = filtro_vecinos(MLH_var_filtro) Filtro_Vecinos_grad = filtro_vecinos(MLH_grad_filtro) CapaLimite_Var_df = pd.DataFrame(Filtro_Vecinos_var,index = Backs.index) CapaLimite_grad_df = pd.DataFrame(Filtro_Vecinos_grad,index = Backs.index) return CapaLimite_Var_df, CapaLimite_grad_df def Estima_Richardson(Fecha_Inicio, Fecha_Fin, freq): # Se refiere a la frecuencia del RWP 5 o 60 Temperature_Zenit, RelativeHumidity_Zenit,VaporDensity_Zenit, Liquid_Zenit = Read_CSV(Fecha_Inicio[:8]) Altura, Fechas, Velocidad, Direccion, Zonal, Meridional, Omega = lee_datos_DiaVentana(Fecha_Inicio, 24, freq) Fechas = pd.date_range(Fecha_Inicio,Fecha_Fin,freq='1D') for i in Fechas[1:]: # Temperatura Radiometro Dia = datetime.strftime(i, '%Y%m%d') try: Temp_Zenit, Hum_Zenit, VaporDensity_Zenit, Liquid_Zenit = Read_CSV(Dia) Temperature_Zenit = Temperature_Zenit.append(Temp_Zenit) RelativeHumidity_Zenit = RelativeHumidity_Zenit.append(Hum_Zenit) # Vientos Radar Dia = datetime.strftime(i, '%Y%m%d %H:%M') Alt, Fech, Vel, Dir, Zon, Meri, Ome = lee_datos_DiaVentana(Dia, 24, freq) Velocidad = Velocidad.append(Vel) Zonal = Zonal.append(Zon) Meridional = Meridional.append(Meri) Direccion = Direccion.append(Dir) Omega = Omega.append(Ome) except: pass if freq == 5: Temperature_Zenit = Temperature_Zenit.resample('5Min', how='mean', label='right', closed='right') RelativeHumidity_Zenit = RelativeHumidity_Zenit.resample('5Min', how='mean', label='right', closed='right') Temperature_Zenit = Temperature_Zenit.reindex(pd.date_range(Fecha_Inicio, Fecha_Fin,freq='5Min')) RelativeHumidity_Zenit = RelativeHumidity_Zenit.reindex(pd.date_range(Fecha_Inicio, Fecha_Fin,freq='5Min')) else: Temperature_Zenit = Temperature_Zenit.resample('1H', how='mean', label='right', closed='right') RelativeHumidity_Zenit = RelativeHumidity_Zenit.resample('1H', how='mean', label='right', closed='right') Temperature_Zenit = Temperature_Zenit.reindex(pd.date_range(Fecha_Inicio, Fecha_Fin,freq='1H')) RelativeHumidity_Zenit = RelativeHumidity_Zenit.reindex(pd.date_range(Fecha_Inicio, Fecha_Fin,freq='1H')) Zonal["index"] = Zonal.index Zonal = Zonal.drop_duplicates(subset='index', keep='last') if freq == 5: Zonal = Zonal.reindex(pd.date_range(Fecha_Inicio, Fecha_Fin,freq='5Min')) else: Zonal = Zonal.reindex(pd.date_range(Fecha_Inicio, Fecha_Fin,freq='1H')) del Zonal['index'] Meridional["index"] = Meridional.index Meridional = Meridional.drop_duplicates(subset='index', keep='last') if freq == 5: Meridional = Meridional.reindex(pd.date_range(Fecha_Inicio, Fecha_Fin,freq='5Min')) else: Meridional = Meridional.reindex(pd.date_range(Fecha_Inicio, Fecha_Fin,freq='1H')) del Meridional['index'] # Valores interpolados para igualar los valores respecto a las alturas del RWP # Altura Referencia RWP # ======================================================================================================== Alt = Altura Altura_RWP = Alt[:,1] Altura_Radio = (np.array(Temperature_Zenit.keys().astype(np.float)) * 1000) + 60 Temperatura_Matriz = [] for i in range(len(Temperature_Zenit)): Temp_perfil = Temperature_Zenit.values[i] Temp_interp = [] for Alt_RWP in Altura_RWP: idd = 0 while not (Alt_RWP < Altura_Radio[idd + 1] and Alt_RWP > Altura_Radio[idd]): idd = idd + 1 Temp = Temp_perfil[idd] + (((Alt_RWP - Altura_Radio[idd])/(Altura_Radio[idd+1] - Altura_Radio[idd])) \ * (Temp_perfil[idd+1] - Temp_perfil[idd])) Temp_interp.append(Temp) Temperatura_Matriz.append(Temp_interp) Humedad_Matriz = [] for i in range(len(RelativeHumidity_Zenit)): Hum_perfil = RelativeHumidity_Zenit.values[i] Hum_interp = [] for Alt_RWP in Altura_RWP: idd = 0 while not (Alt_RWP < Altura_Radio[idd + 1] and Alt_RWP > Altura_Radio[idd]): idd = idd + 1 Hum = Hum_perfil[idd] + (((Alt_RWP - Altura_Radio[idd])/(Altura_Radio[idd+1] - Altura_Radio[idd])) \ * (Hum_perfil[idd+1] - Hum_perfil[idd])) Hum_interp.append(Hum) Humedad_Matriz.append(Hum_interp) # DataFrame con las Temperaturas interpoladas... Datos con los que se trabajará de ahora en adelante # Parámetros # ======================================================================================================== Temperature_Zenit = pd.DataFrame(Temperatura_Matriz, columns = Altura_RWP, index = Temperature_Zenit.index) RelativeHumidity_Zenit = pd.DataFrame(Humedad_Matriz, columns = Altura_RWP, index = RelativeHumidity_Zenit.index) Altura = Altura_RWP Press_est = 1013.25 * (1 - (2.25577e-5 * (Altura + 1547))) ** 5.25588 # Expresión cálculo Press atmosférica Theta = Temperature_Zenit * ((1000/np.array(Press_est))**0.286) # Estimación Temperatura Potencial #Dew Point Num = (243.04 * (np.log(RelativeHumidity_Zenit/100.)+((17.625 * Temperature_Zenit)/(243.04 + Temperature_Zenit)))) Den = 17.625 - np.log(RelativeHumidity_Zenit /100.)-(17.625 * Temperature_Zenit/(243.04 + Temperature_Zenit)) Dew_p = Num / Den # Mixing Ratio Sd = (7.5 * Dew_p) / (237.7 + Dew_p) a = 6.11 * (10 ** (Sd)) Mix_Rat = ((0.622 * a) / (Press_est - a)) # Theta_v Theta_v = Theta * (1 + (0.61 * Mix_Rat)) # Perfiles Verticales Individuales # ====================================================================================== if freq == 5: Fechas = pd.date_range(Fecha_Inicio,Fecha_Fin,freq='5Min') else: Fechas = pd.date_range(Fecha_Inicio,Fecha_Fin,freq='1H') MLH = [] for i in Fechas: Dia = datetime.strftime(i, '%Y%m%d') Hora = datetime.strftime(i, '%H:%M') # UTC # Radiometro Temp_perfil = Temperature_Zenit[Dia].between_time(Hora,Hora).values[0] Theta_perfil = Theta[Dia].between_time(Hora,Hora).values[0] Thetav_perfil = Theta_v[Dia].between_time(Hora,Hora).values[0] # Perfilador Zonal_perfil = Zonal[Dia].between_time(Hora,Hora).values[0] Meridional_perfil = Meridional[Dia].between_time(Hora,Hora).values[0] Ri_1 = Richardson_comp(Altura, Thetav_perfil, Zonal_perfil, Meridional_perfil) Rc = 1 try: # Interpolación para mejorar la resolución Alt_2 = Altura[np.where(Ri_1 > Rc)[0][0]] Alt_1 = Altura[np.where(Ri_1 > Rc)[0][0] - 1] R2 = Ri_1[Altura == Alt_2] R1 = Ri_1[Altura == Alt_1] m = (Alt_2 - Alt_1)/ (R2 - R1) MLH.append((((Rc - R1 ) * m ) + Alt_1)[0]) except: MLH.append(np.NaN) MLH = pd.DataFrame(MLH, Fechas) return MLH def Grafica_Perfil(Perfil, Ventana): if Ventana < 5: minute =mdates.MinuteLocator(interval = 25) elif Ventana == 6: minute =mdates.MinuteLocator(interval = 60) elif 7 < Ventana < 20: minute =mdates.MinuteLocator(interval = 90) elif 20 < Ventana < 48: minute =mdates.MinuteLocator(interval = 180) elif 48 <= Ventana: minute =mdates.MinuteLocator(interval = 500) plt.close('all') fig = plt.figure(figsize=[14,6]) ax = fig.add_subplot(111) ax = fig.add_subplot(111) ax.spines['bottom'].set_color((0.45, 0.45, 0.45)) ax.spines['top'].set_color((1,1,1)) ax.xaxis.label.set_color((0.45, 0.45, 0.45)) ax.tick_params(axis='x', colors=(0.45, 0.45, 0.45)) ax.spines['right'].set_color((1,1,1)) ax.spines['left'].set_color((0.45, 0.45, 0.45)) ax.yaxis.label.set_color((0.45, 0.45, 0.45)) ax.tick_params(axis='y', colors=(0.45, 0.45, 0.45)) ax.tick_params(right="off") ax.tick_params(left="on",width=1.5) ax.tick_params(top="off") Perfil_50 = Perfil[4] Perfil_50[Perfil_50.values > 1100] = np.NaN ax.plot(Perfil_50.index.to_pydatetime(), Perfil_50, lw='1.7', color=(0.97647059, 0.90196078, 0.22352941),label=u'50 metros') Perfil_300 = Perfil[49] Perfil_300[Perfil_300.values > 1100] = np.NaN ax.plot(Perfil_300.index.to_pydatetime(), Perfil_300, lw='1.7', color=(0.54509804, 0.73333333, 0.45490196),label=u'500 metros') Perfil_1000 = Perfil[99] Perfil_1000[Perfil_1000.values > 600] = np.NaN ax.plot(Perfil_1000.index.to_pydatetime(), Perfil_1000, lw='1.7', color=(0.16862745, 0.28235294, 0.41176471),label=u'1000 metros') ax.set_title(u'CEILÓMETRO ' + Nombre, fontsize=17,fontproperties=prop, color=(0.45, 0.45, 0.45)) ax.set_ylim([0,1100]) ax.set_ylabel(ur'Intensidad Backscatter [$10^{-9} m^{-1} sr^{-1}$]', fontsize=17, fontproperties=prop_1) plt.xticks(fontsize=17, fontproperties=prop_1) ax.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M\n\n%d-%b')) ax.xaxis.set_major_locator(minute) # Leyenda legend = ax.legend(bbox_to_anchor=(0., -0.55, 1, .3), loc=9,ncol=3, mode="expand", borderaxespad=0., prop={'size':15}, shadow=False) frame = legend.get_frame() frame.set_edgecolor('0.90') for text in legend.get_texts(): plt.setp(text, color = (0.45, 0.45, 0.45)) plt.savefig("perfil_Ceilo.png", bbox_inches='tight', dpi=150) os.system("scp perfil_Ceilo.png lherrera@siata.gov.co:/var/www/laura/Pagina/Ceilometros/Perfil"+str(Ceilometro)+'_'+str(Ventana)+".png") # ========================================================================================================== # Genera las Gráficas #for Nombre, Ceilometro in zip(['TorreSiata'],['siata']): for Nombre, Ceilometro in zip([],[]): try: print Nombre Backs_Ceil = lee_datos(Ceilometro, 72) CLA_Varianza, CLA_gradiente = CLA_Ceil(Backs_Ceil, '30') try: CLA_Richardson = Estima_Richardson(datetime.strftime(Backs_Ceil.index[0],'%Y%m%d %H:%M'),datetime.strftime(Backs_Ceil.index[-1],'%Y%m%d %H:%M'),5) except: traceback.print_exc() pass for Ventana in [3,12,24,48,72]: print Ventana try: Grafica_Ventanas(Backs_Ceil, Ventana, Nombre, CLA = CLA_gradiente, Nombre_CLA = 'Gradiente') Grafica_Ventanas(Backs_Ceil, Ventana, Nombre, CLA = CLA_Varianza, Nombre_CLA = 'Varianza') Grafica_Ventanas(Backs_Ceil, Ventana, Nombre) Grafica_Ventanas(Backs_Ceil, Ventana, Nombre, CLA = CLA_Richardson,Nombre_CLA = 'Rich') except: traceback.print_exc() # pass pass try: Backs_Perfil = Backs_Ceil[Backs_Ceil.index[-1] - timedelta(hours = Ventana):Backs_Ceil.index[-1]] Grafica_Perfil(Backs_Perfil, Ventana) except: traceback.print_exc() # pass pass except: traceback.print_exc() # pass pass #for Nombre, Ceilometro in zip(['AMVA','Itagui'],['amva','itagui']): for Nombre, Ceilometro in zip(['Itagui'],['itagui']): print Nombre try: Backs_Ceil = lee_datos(Ceilometro, 72) CLA_Varianza, CLA_gradiente = CLA_Ceil(Backs_Ceil, '30') for Ventana in [3,12,24,48,72]: print Ventana try: Grafica_Ventanas(Backs_Ceil, Ventana, Nombre, CLA = CLA_gradiente, Nombre_CLA = 'Gradiente') Grafica_Ventanas(Backs_Ceil, Ventana, Nombre, CLA = CLA_Varianza, Nombre_CLA = 'Varianza') Grafica_Ventanas(Backs_Ceil, Ventana, Nombre) Grafica_Ventanas(Backs_Ceil, Ventana, Nombre, CLA = CLA_Richardson,Nombre_CLA = 'Rich') except: traceback.print_exc() pass try: Backs_Perfil = Backs_Ceil[Backs_Ceil.index[-1] - timedelta(hours = Ventana):Backs_Ceil.index[-1]] Grafica_Perfil(Backs_Perfil, Ventana) except: traceback.print_exc() pass except: traceback.print_exc()