Przez ostatnie dni media zdominował temat tragicznej jakości powietrza w polskich miastach. Nie jest to problem wirtualny, bo wystarczy wyjść na dwór by zobaczyć i poczuć marną jakość powietrza. Do wytworzenia smogu potrzebujemy dwóch czynników. Pierwszy to produkcja szkodliwych pyłów i zanieczyszczeń, która zimą jest mocno nasilona ze względu na sezon grzewczy, szczególnie w indywidualnych instalacjach C.O. z piecem na paliwo stało (drewno, węgiel), często opalanym odpadami. Elektrociepłownie, mimo, że opalane węglem wydalają proporcjonalnie mniej zanieczyszczeń na kW energii ze względu na zastosowane w instalacji filtry. Drugim, niezbędnym do powstania smogu warunkiem jest wystąpienie warstwy inwersyjnej w atmosferze, która blokuje wymianę powierza przy ziemi, z tym powyżej granicy inwersji. Jeżeli wysokość warstwy inwersyjnej wynosi 100-200 metrów, to znajdujemy się w zamkniętym “słoiku”, do którego wydalamy ogromne ilości zanieczyszczeń. Tyle skróconej teorii, przejdźmy do praktyki.
Wykorzystując bardzo duże stężenie pyłów w Warszawie postanowiłem sprawdzić jak problem wygląda z góry. Udało mi się wykonać dwa bardzo wymowne ujęcia:
Zdjęcie po lewej przedstawia kominy elektrociepłowni Siekierki i warszawskie drapacze chmur w tle. Zdjęcie po prawej to widok na centrum Warszawy z okolicy mostu Gdańskiego. Niestety zdjęcia nie pokazują chmur a właśnie zanieczyszczenia, którymi byliśmy zmuszeni oddychać.
Równolegle uruchomiłem pomiar stężenia pyłów w powietrzu na zewnątrz domu w Legionowie. Zakupiłem laserowy czujnik i połączyłem go z RaspberryPi. Czujnik wisi sobie bezpośrednio za oknem i komunikuje się z jeżyną przez interfejs szeregowy:
Każdy pomiar daje trzy wyniki odpowiadające pyłom o różnych frakcjach, w tym “najpopularniejszy” PM2.5. Program odbierający pomiary napisany w Pythonie:
import serial
from time import gmtime, strftime
port = serial.Serial("/dev/serial0", baudrate=9600, timeout=1.5)
data = port.read(32);
f = open("/home/pi/" + strftime("%Y-%m-%d", gmtime()) + ".txt", "a")
if ord(data[0]) == 66 and ord(data[1])==77:
suma = 0
for a in range(30):
suma += ord(data[a])
if suma == ord(data[30])*256+ord(data[31]):
PM01 = str(ord(data[4])*256+ord(data[5]))
PM25 = str(ord(data[6])*256+ord(data[7]))
PM10 = str(ord(data[8])*256+ord(data[9]))
str = strftime("%Y-%m-%d %H:%M:%S", gmtime()) + "\t" + PM01 + "\t" + PM25 + "\t" + PM10;
f.write(str + "\n")
print(str)
port.close()
f.close()
Powyższy kod uruchamiany jest co minutę. Dane zapisywane są w zwykłym pliku tekstowym:
2017-01-09 20:43:02 185 470 729
2017-01-09 20:44:02 181 440 682
2017-01-09 20:45:02 175 460 701
2017-01-09 20:46:02 172 445 725
2017-01-09 20:47:02 175 448 680
2017-01-09 20:48:02 177 449 725
2017-01-09 20:49:02 170 450 754
2017-01-09 20:50:02 176 458 746
2017-01-09 20:51:02 165 441 733
2017-01-09 20:52:02 174 434 707
2017-01-09 20:53:02 170 451 715
2017-01-09 20:54:02 170 436 702
Następnie raz na 15 minut uruchamiany jest kod rysujący wykres. Dla poprawy czytelności rysowane są średnie kroczące z pomiarów, a nie same pomiary:
import numpy
import ftplib
from ftplib import FTP
import sys
from matplotlib.dates import strpdate2num
from datetime import datetime, timedelta
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
date = sys.argv[1];
X = numpy.loadtxt('/root/AIR/OLSZANKOWA50A/'+date+'.txt', delimiter='\t', dtype=object, converters={0: lambda x: datetime.strptime(x.decode("utf-8"), "%Y-%m-%d %H:%M:%S"), 1: numpy.float, 2: numpy.float, 3: numpy.float})
window = 900;
step = 60;
ts = min(X[:,0])
tss = datetime(ts.year, ts.month, ts.day)
t = tss;
te = t + timedelta(days=1)
T = []
PM01 = []
PM25 = []
PM10 = []
while t < te:
t = t + timedelta(seconds=step)
idx = numpy.logical_and(X[:,0] > (t - timedelta(seconds=window)), X[:,0] < (t + timedelta(seconds=window)))
#print(t, sum(idx))
T.append(t)
if sum(idx) > 0:
PM01.append(numpy.mean(X[idx,1]))
PM25.append(numpy.mean(X[idx,2]))
PM10.append(numpy.mean(X[idx,3]))
else:
PM01.append(numpy.nan)
PM25.append(numpy.nan)
PM10.append(numpy.nan)
x = []
label = []
for a in range(0,24,1):
x.append(tss+timedelta(hours=a))
if a%2 ==0:
label.append((tss+timedelta(hours=a)).strftime('%H:%M:%S'))
else:
label.append('')
plt.figure(figsize=(20, 10))
plt.plot(T,PM01, 'b', lw=2, label='PM 1.00')
plt.plot(T,PM25, 'r', lw=4, label='PM 2.50')
plt.plot(T,PM10, 'g', lw=2, label='PM 10.0')
plt.xlabel('Time [UTC]')
plt.ylabel('ug/m3')
plt.legend(loc='upper left', shadow=True)
plt.xlim([ts, te])
plt.ylim([0, (numpy.ceil(max(PM10)/50)+1)*50])
plt.grid()
plt.title(date)
plt.xticks(x,label, rotation=0)
plt.savefig('/root/AIR/OLSZANKOWA50A/'+date+'.png',dpi=300)
File2Send = '/root/AIR/OLSZANKOWA50A/'+date+'.png'
Output_Directory = "/public_html/parking-legionowo.pl/AIR/"
ftp = FTP('***************')
ftp.login('***************', '***************')
file = open(File2Send, "rb")
ftp.cwd(Output_Directory)
ftp.storbinary('STOR ' +date+'.png', file)
print("STORing File now...")
ftp.quit()
file.close()
print("File transfered")
Po przygotowaniu wykresu jest on zapisywany do pliku png i wysyłany na serwer, gdzie można na bieżąco podglądać pomiary.
Tak wygląda normalny dzień “bez smogu” (te chwilowe wzrosty to prawdopodobnie momenty, gdy sąsiedzi dokładają węgla do starego pieca):
A tak wygląda dzień “smogowy” (akurat wtedy robiłem zdjęcia pokazane na początku):
Aktualne pomiary można zobaczyć na stronie sponsora: http://parking-legionowo.pl/AIR/
