Smogowy armageddon

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/

|

Komentarze (3) »

  1. Qavtan:

    12 Jan 2017 @ 12:29

    Jak wiarygodne są takie wyniki w kontekście porównania ze stacjami pomiarowymi np. wojewódzkich inspektoratów ochrony środowiska?

    Gratuluję zacięcia ZróbToSam :]

  2. Wiktor:

    13 Jan 2017 @ 07:56

    Hej,

    Nie podałeś jaki masz czujnik. W zależności od czujnika ramka jest inna. W czynnikach jakie ja mam PMS3003 i PMS7003 ramki są bardzo podobne. W bajtach 2-3 jest podana długość ramki, a dane są w tym samym miejscu. Możesz przerobić ten kod na bardziej uniwersalny.

    Pozdrawiam

  3. jacek:

    20 Sep 2017 @ 20:04

    chciałem coś podobnego zrobić mam czujnik PMS7003, generalnie chciałem to odczytywać komendą:
    sudo od –endian=big -x -N10 < /dev/ttyAMA0
    z opisy wynika że początek transmisji powinien się zaczynać od 424d , ale jak kolejne razy odpalam komendę daje coś takiego:
    3500 8200 0a54 008c 000a
    8100 9f00 0a78 0094 000a
    2900 6400 2900 6400 2900

    jeśli odczyt robie przez pythona:
    import serial
    from time import gmtime, strftime
    port = serial.Serial("/dev/serial0", baudrate=9600, timeout=1.5)
    data = port.read(32);
    PM01 = str(ord(data[5])*256+ord(data[6]))
    PM25 = str(ord(data[7])*256+ord(data[8]))
    PM10 = str(ord(data[9])*256+ord(data[10]))
    print(PM01)
    print(PM25)
    print(PM10)

    to też podaje błędne wartości, ma ktoś jakiś pomysł? dlaczego transmisja każdego pakietu danych nie zaczyna sie od 424d?

RSS komentarzy · adres TrackBack

Odpowiedz