estatisticas.py 14.1 KB
Newer Older
1
# -*- coding:utf-8 -*-
2
# Programa de visualização e estatísticas de dados meteorológicos
3
# coletados pela estação Meteorolog do Centro de Tecnologia Acadêmica
4
# Dados no formato DataHora Temperatura Umidade_do_Ar Solo Pressão Luminosidade
5 6 7 8 9 10 11 12 13 14 15 16 17 18
# DataHora no formato  AAAAMMDDHHmmSS onde
#          AAAA = Ano
#          MM = Mês
#          DD = Dia
#          HH = Hora
#          mm = Minuto
#          ss = Segundo
#
# Mais informações:
#  http://cta.if.ufrgs.br/projects/estacao-meteorologica-modular/wiki/Meteorolog
#
# Autor: Rafael Pezzi - Centro de Tecnologia Acadêmica - IF/UFRGS
#
# Este programa é softare livre, disponível sob os termos da
19
# GNU Affero General Public License V3
20 21 22 23 24 25 26 27 28 29
# http://www.gnu.org/licenses/agpl-3.0.html

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import matplotlib.mlab as mlab
import random as random
import datetime as datetime
import sys , os

30 31 32
# Realizar médias de n_medidas aleatórias:
#n_medias = [1,5,10,50,100,1000,5000,10000]
#n_medias = [5,10,50,100,1000,5000]
33
n_medias = [5,10,50,100,1000]
34
#n_medias = [5,10] # Possibilidades reduzidas para debug
35

36
# Número de medias a realizar pelo método de Monte Carlo
37 38
n_MC=5000

39 40 41 42 43 44
# deg_of_freedom=0 para cálculo dos desvios padrão da população
# deg_of_freedom=1 para cálculo dos desvios padrão da amostra (com correção de Bessel)
deg_of_freedom=1
# Ref: http://docs.scipy.org/doc/numpy/reference/generated/numpy.std.html#numpy.std


45 46
print "Estatísticas Meteorológicas Meteorolog - Centro de Tecnologia Acadêmica"

47
# Testa se o número de argumentos está correto
48 49 50 51 52 53 54 55 56 57 58 59 60 61
if len(sys.argv) < 2:
    sys.stderr.write('Uso: python '+ sys.argv[0]+' arquivo_de_log\n' )
    sys.exit(1)
# e se o arquivo de log existe
if not os.path.exists(sys.argv[1]):
    sys.stderr.write('ERRO: Arquivo '+sys.argv[1]+' não encontrado!\n')
    sys.exit(1)

# Ativa/desativa seleção de intervalo
# Caso negativo (False), todos dados do arquivo de log serão computados
seleciona_intervalo = False
#seleciona_intervalo = True

# Defina aqui o intervalo de tempo a se analisado
62
#  TODO: Incluir seleção pela linha de comando de execução
63 64 65 66 67 68 69
sel_Inicio = datetime.datetime(2013,03,01,00,00,00)
sel_Fim = datetime.datetime(2013,03,31,23,59,59)

# Armazena os dados coletados em matrizes unidimensionais
data, temp, u_ar, u_solo, pressao, lum = np.loadtxt(sys.argv[1], unpack=True,
        converters={ 0: mdates.strpdate2num('%Y%m%d%H%M%S')})

70
## Seria útil incluir um algoritmo de filtragem para remover pontos expurios
71 72 73 74 75 76 77 78 79 80 81 82
##   observados em algumas partes por mil medidas de temperatura

# Exibe informações dos dados
print "Dados de "+str(mdates.num2date(data[0]).strftime('%Y-%m-%d %H:%M:%S'))+u" até "+str(mdates.num2date(data[len(data)-1]).strftime('%Y-%m-%d %H:%M:%S'))
print "Série com "+str(len(data))+" amostras\n"

# Seleciona dados em um intervalo de tempo
if (seleciona_intervalo == True):

    # Cria listas temporárias para seleção
    # Acho que pode ser otimizado, usando iteradores direto em arrays
    #  crítico para series de dados gigantes
83 84
    selec_data=[]; selec_temp=[];  selec_u_ar=[]
    selec_u_solo=[]; selec_pressao=[]; selec_lum=[]
85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106
    for j in enumerate(data):
        if mdates.date2num(sel_Inicio) <= j[1] <= mdates.date2num(sel_Fim):
            selec_data.append(data[j[0]])
            selec_temp.append(temp[j[0]])
            selec_u_ar.append(u_ar[j[0]])
            selec_u_solo.append(u_solo[j[0]])
            selec_pressao.append(pressao[j[0]])
            selec_lum.append(lum[j[0]])

    # Copia dados selecionados para matrizes originais
    data=np.asarray(selec_data)
    temp=np.asarray(selec_temp)
    u_ar=np.asarray(selec_u_ar)
    u_solo=np.asarray(selec_u_solo)
    pressao=np.asarray(selec_pressao)
    lum=np.asarray(selec_lum)

    # Exibe informações sobre o periodo selecionado
    print "Selecionado período:"
    print str(sel_Inicio)+" até "+str(sel_Fim)
    print "Série com "+str(len(data))+" amostras\n"

107 108
print "\n############## Séries Temporais e Histogramas ##################"

109 110 111 112
## Análise da temperatura
T_min=temp.min()
T_max=temp.max()
T_med=temp.mean()
113
T_std=temp.std(ddof=deg_of_freedom)
114 115 116 117 118 119
print "Temperatura: \nmin: "+str(T_min)+"°C\tmax: "+str(T_max)+"°C\tmédia: "+str("%.2f" % T_med)+"°C\tDesvio: "+str("%.1f" % T_std)+"°C"

## Análise da umidade relativa do Ar
U_ar_min=u_ar.min()
U_ar_max=u_ar.max()
U_ar_med=u_ar.mean()
120
U_ar_std=u_ar.std(ddof=deg_of_freedom)
121 122 123 124 125 126
print "\nUmidade Relativa do Ar: \nmin: "+str(U_ar_min)+"%\tmax: "+str(U_ar_max)+"%\tmédia: "+str("%.1f" % U_ar_med)+"%\tDesvio: "+str("%.1f" % U_ar_std)+"%"

# Análise da pressão atomsférica
P_min=pressao.min()
P_max=pressao.max()
P_med=pressao.mean()
127
P_std=pressao.std(ddof=deg_of_freedom)
128 129 130 131 132 133 134

print "\nPressão: \nmin: "+str(int(P_min))+" Pa\tmax: "+str(int(P_max))+" Pa\tmédia: "+ str("%.0f" % P_med)+" Pa\tDesvio: "+str("%.2f" % P_std)+" Pa"

# Análise da luminosidade
L_min=lum.min()
L_max=lum.max()
L_med=lum.mean()
135
L_std=lum.std(ddof=deg_of_freedom)
136 137 138 139 140 141 142 143 144 145 146 147

print "\nLuminosidade: \nmin: "+str(int(L_min))+" u.a.\tmax: "+str(int(L_max))+" u.a.\tmédia: "+ str("%.0f" % L_med)+" u.a.\tDesvio: "+str("%.2f" % L_std)+" u.a."


# Apresentação de histrograma:
# Temperatura
plt.figure() # Cria nova figura
plt.ylabel(u"Número de eventos")
plt.xlabel(u"Temperatura (°C)")
plt.hist(temp,bins=25)
#Inclui título com média e desvio padrão
#plt.text(plt.xlim()[0]+0.5,0.9*plt.ylim()[1],u"Média: "+str("%.2f" % temp.mean())+u"°C\nDesvio Padrão: "+str("%.2f" % temp.std())+u"°C") # Este coloca dentro da figura
148
plt.title(u"Temperatura média: "+str("%.2f" % T_med)+u"°C  |  Desvio Padrão: "+str("%.2f" % temp.std(ddof=deg_of_freedom))+u"°C")
149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167
plt.grid(True)

# Apresentação de histrograma:
# Umidade Rel. do Ar
plt.figure()
plt.ylabel(u"Número de eventos")
plt.xlabel("Umidade Relativa (%)")
plt.hist(u_ar,bins=range(26,100,2))
#plt.text(plt.xlim()[0]+1,0.9*plt.ylim()[1],u"Média: "+str("%.2f" % u_ar.mean())+u" %\nDesvio Padrão: "+str("%.2f" % u_ar.std())+" %")
plt.title(u"UR média: "+str("%.2f" % u_ar.mean())+u" %  |  Desvio Padrão: "+str("%.2f" % u_ar.std())+" %")
plt.grid(True)

# Apresentação de histrograma:
# Pressão
plt.figure()
plt.ylabel(u"Número de eventos")
plt.xlabel(u"Pressão atmosférica (Pa)")
plt.hist(pressao,bins=25)
#plt.text(1.001*plt.xlim()[0],0.9*plt.ylim()[1],u"Média: "+str("%.2f" % pressao.mean())+u" Pa\nDesvio Padrão: "+str("%.2f" % pressao.std())+" Pa")
168
plt.title(u"Pressão média: "+str("%.2f" % pressao.mean())+u" Pa  |  Desvio Padrão: "+str("%.2f" % pressao.std(ddof=deg_of_freedom))+" Pa")
169 170 171 172 173 174 175 176 177
plt.grid(True)

# Apresentação de histrograma:
# Luminosidade
plt.figure()
plt.ylabel(u"Número de eventos")
plt.xlabel(u"Luminosidade (u.a.)")
plt.hist(lum,bins=25)
#plt.text(1*plt.xlim()[0]+60,0.9*plt.ylim()[1],u"Média: "+str("%.2f" % lum.mean())+u"\nDesvio Padrão: "+str("%.2f" % lum.std())+" ")
178
plt.title(u"Luminosidade média: "+str("%.2f" % lum.mean())+u"  |  Desvio Padrão: "+str("%.2f" % lum.std(ddof=deg_of_freedom))+" ")
179 180 181 182 183 184 185 186 187
plt.grid(True)

## Apresentação das series temporais de medidas
# Temperatura
plt.figure()
# Gira o texto do eixo temporal
plt.xticks(rotation=70)
plt.plot_date(x=data, y=temp, fmt="r.")

188
# Define os limites temporais do gráfico
189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236
#   útil para evitar escalas exageradas devido a dadas incorretas - ocorreu na estação sem relógio
plt.xlim(data[0], data[len(data)-1])
plt.ylabel(u"Temperatura (°C)")
plt.grid(True)
# Altera as margens para a margem inferior comportar a data
plt.subplots_adjust(left=0.1, right=0.9, top=0.9, bottom=0.2)

# Umidade Relativa do Ar
plt.figure()
plt.xticks(rotation=70)
plt.plot_date(x=data, y=u_ar, fmt="r.")
plt.xlim(data[0], data[len(data)-1])
plt.ylabel("Umidade relativa do ar (%)")
plt.grid(True)
plt.subplots_adjust(left=0.1, right=0.9, top=0.9, bottom=0.2)

'''
# Condutividade do solo
plt.figure()
plt.xticks(rotation=70)
plt.plot_date(x=data, y=u_solo, fmt="r.")
plt.xlim(data[0], data[len(data)-1])
plt.ylabel("Umidade do solo (u.a.)")
plt.grid(True)
plt.subplots_adjust(left=0.1, right=0.9, top=0.9, bottom=0.2)
'''

# Pressão
plt.figure()
plt.xticks(rotation=70)
plt.plot_date(x=data, y=pressao, fmt="r.")
plt.xlim(data[0], data[len(data)-1])
plt.ylabel(u"Pressão atmosférica (Pa)")
plt.grid(True)
plt.subplots_adjust(left=0.13, right=0.9, top=0.9, bottom=0.2)

# Luminosidade
plt.figure()
plt.xticks(rotation=70)
plt.plot_date(x=data, y=lum, fmt="r.")
plt.xlim(data[0], data[len(data)-1])
plt.ylabel("Luminosidade (u.a.)")
plt.grid(True)
plt.subplots_adjust(left=0.1, right=0.9, top=0.9, bottom=0.2)


#'''
# Análise Monte Carlo para as médias
237
# Para fins de ilustração didática do teorema do limite central e desvio da média
238 239 240
print "\n############## Simulação Monte Carlo ##################"
print n_MC, "médias calculadas para", n_medias, "medidas"

241 242

# Cria variaveis para calculo das médias
243 244 245 246 247 248
temp_std_n=[]
umid_std_n=[]
pressao_std_n=[]
lum_std_n=[]

for tentativas in n_medias:
249 250 251 252 253 254 255 256
 temp_med = 0.
 temp_med_vec=[]
 umid_med = 0.
 umid_med_vec=[]
 pressao_med = 0.
 pressao_med_vec=[]
 lum_med = 0.
 lum_med_vec=[]
257
 for j in range(n_MC):
258 259 260 261 262 263 264 265 266 267 268 269 270 271 272
   for n in range(tentativas):
       amostra=int(random.random()*len(data))
       lum_med+=lum[amostra]
       temp_med+=temp[amostra]
       umid_med+=u_ar[amostra]
       pressao_med+=pressao[amostra]

#   print temp_med/tentativas
#   print umid_med/tentativas
#   print pressao_med/tentativas
#   print ("\n")
   temp_med_vec.append(temp_med/tentativas)
   umid_med_vec.append(umid_med/tentativas)
   pressao_med_vec.append(pressao_med/tentativas)
   lum_med_vec.append(lum_med/tentativas)
273

274
   temp_med = 0
275
   umid_med =0
276 277
   lum_med = 0
   pressao_med = 0
278 279 280 281 282 283 284 285 286 287


 # Convertentndo listas para matrizes para melhorar desempenho (pode - deve - ser feito anteriormente no programa)
 temp_med_vec=np.asarray(temp_med_vec)
 umid_med_vec=np.asarray(umid_med_vec)
 pressao_med_vec=np.asarray(pressao_med_vec)
 lum_med_vec=np.asarray(lum_med_vec)



288
 # Armazenando as médias dos resultados do cálculo Monte Carlo
289 290 291 292 293 294 295 296 297
 temp_med =  temp_med_vec.mean()
 umid_med = umid_med_vec.mean()
 lum_med = lum_med_vec.mean()
 pressao_med = pressao_med_vec.mean()
 lum_med_vec=np.asarray(lum_med_vec)


 print " "
 print str(tentativas)+" medidas:"
298
 temp_med_std=temp_med_vec.std(ddof=deg_of_freedom)
299 300 301 302
 temp_std_n.append(temp_med_std)
 print "Desvio Temp Média para MC com "+str(tentativas)+" medidas: "+str(temp_med_std)

# umid_med_vec.append(umid_med/tentativas)
303
 umid_med_std=umid_med_vec.std(ddof=deg_of_freedom)
304 305 306 307
 umid_std_n.append(umid_med_std)
 print "Desvio Umidade Média para MC com "+str(tentativas)+" medidas: "+str(umid_med_std)

# pressao_med_vec.append(pressao_med/tentativas)
308
 pressao_med_std=pressao_med_vec.std(ddof=deg_of_freedom)
309 310 311
 pressao_std_n.append(pressao_med_std)
 print "Desvio Pressão Média para MC com "+str(tentativas)+" medidas: "+str(pressao_med_std)

312
 lum_med_std=lum_med_vec.std(ddof=deg_of_freedom)
313 314 315 316
 lum_std_n.append(lum_med_std)
 print "Desvio lumin. Média para MC com "+str(tentativas)+" medidas: "+str(lum_med_std)


317
 plt.figure()
318 319 320 321
 plt.title(u"Temp. MC  "+str(tentativas)+" medidas: sigma = "+str("%0.3f" % temp_med_std)+u" | Média = "+str(temp_med))
 n, bins, patches =  plt.hist(temp_med_vec, bins=50,normed=1)
 y = mlab.normpdf( bins,temp_med_vec.mean() , temp_med_std)
 l = plt.plot(bins, y, 'k--', linewidth=1.5)
322 323 324
 plt.xlim([T_min,T_max])
# plt.text(1.001*plt.xlim()[0],0.9*plt.ylim()[1],u"Média: "+str("%.2f" % temp_med_vec.mean())+u" °C\nDesvio Padrão: "+str("%.2f" % temp_med_vec.std())+" °C")
 plt.figure()
325 326 327 328 329
# plt.title(u"Desvio Umidade Média para MC com "+str(tentativas)+" medidas: "+str(umid_med_std))
 plt.title(u"Umid. MC  "+str(tentativas)+" medidas: sigma = "+str("%0.3f" % umid_med_std)+u" | Média = "+str(umid_med))
 n, bins, patches = plt.hist(umid_med_vec, bins=50,normed=1)
 y = mlab.normpdf( bins,umid_med_vec.mean() , umid_med_std)
 l = plt.plot(bins, y, 'k--', linewidth=1.5)
330 331 332
 plt.xlim([U_ar_min,U_ar_max])
# plt.text(1.001*plt.xlim()[0],0.9*plt.ylim()[1],u"Média: "+str("%.2f" % umid_med_vec.mean())+u" %\nDesvio Padrão: "+str("%.2f" % umid_med_vec.std())+" %")
 plt.figure()
333 334 335 336 337 338
# plt.title(u"Desvio Pressão Média para MC com "+str(tentativas)+" medidas: "+str(pressao_med_std))
 plt.title(u"Pressão MC  "+str(tentativas)+" medidas: sigma = "+str("%0.3f" % pressao_med_std)+u" | Média = "+str(pressao_med))

 n, bins, patches = plt.hist(pressao_med_vec, bins=50,normed=1)
 y = mlab.normpdf( bins,pressao_med_vec.mean() , pressao_med_std)
 l = plt.plot(bins, y, 'k--', linewidth=1.5)
339 340 341
 plt.xlim([P_min,P_max])
# plt.text(1.001*plt.xlim()[0],0.9*plt.ylim()[1],u"Média: "+str("%.2f" % pressao_med_vec.mean())+u" Pa\nDesvio Padrão: "+str("%.2f" % pressao_med_vec.std())+" Pa")
 plt.figure()
342 343 344 345 346 347
# plt.title(u"Desvio lumin. Média para MC com "+str(tentativas)+" medidas: "+str(lum_med_std))
 plt.title(u"Lum. MC  "+str(tentativas)+" medidas: sigma = "+str("%0.3f" % lum_med_std)+u" | Média = "+str(lum_med))

 n, bins, patches = plt.hist(lum_med_vec, bins=50, normed=1)
 y = mlab.normpdf( bins,lum_med_vec.mean() , lum_med_std)
 l = plt.plot(bins, y, 'k--', linewidth=1.5)
348 349 350 351 352 353 354 355 356 357 358 359
 plt.xlim([L_min,L_max])


# Todo: Ajuste Gaussiano dos histogramas
# x = np.arange(L_min,L_max,0.1)
# print x
# plt.plot(x)
# plt.plot(x,1000*mlab.normpdf(x,L_med,L_std))

#plt.show()
#'''

360 361 362
#print n_medias
#print temp_std_n

363
## Desvio padrão das médias calculadas pelo método de Monte Carlo
364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390

plt.figure()
plt.title(u"Desvio padrão para Temperatura média de N medidas")
plt.scatter(n_medias, temp_std_n)
for i,txt in enumerate(temp_std_n):
    plt.annotate("%.3f" % txt, (n_medias[i]*1.01,temp_std_n[i]*1.01))

#for i, txt in enumerate(n):
#    ax.annotate(txt, (z[i],y[i]))

plt.figure()
plt.title(u"Desvio padrão para Pressão média de N medidas")
plt.scatter(n_medias, pressao_std_n)
for i,txt in enumerate(pressao_std_n):
    plt.annotate("%.3f" % txt, (n_medias[i]*1.01,pressao_std_n[i]*1.01))

plt.figure()
plt.title(u"Desvio padrão para Umidade Relativa média de N medidas")
plt.scatter(n_medias, umid_std_n)
for i,txt in enumerate(umid_std_n):
    plt.annotate("%.3f" % txt, (n_medias[i]*1.01,umid_std_n[i]*1.01))

plt.figure()
plt.title(u"Desvio padrão para Luminosidade média de N medidas")
plt.scatter(n_medias, lum_std_n)
for i,txt in enumerate(lum_std_n):
    plt.annotate("%.3f" % txt, (n_medias[i]*1.01,lum_std_n[i]*1.01))
391 392


393 394
#Apresenta os gráficos gerados
plt.show()