Jakość powietrza w Polsce #1 – (re)Start

W 2017 roku zacząłem się interesować tematem jakości powietrza. Pchnął mnie do tego fragment jakiegoś materiału w jakimś programie informacyjnym. Nie pamiętam, o co w nim dokładnie chodziło, pamiętam jednak, że wydawało mi się, że dość mocno w nim manipulują, jeśli chodzi o interpretację danych. Było słabo, a prowadząca mówiła, że jest dobrze. Poczułem się odpowiedzialny za zbadanie tematu na własną rękę.

Wziąłem więc ten temat na blogowy tapet i napisałem nawet trochę postów. Leżą one na moim starym porzuconym już blogu tutaj. Posty te zdecydowanie nie zamykają tematu, ale stanowią chyba fajną bazę do kontynuacji moich rozważań.

Post, który czytasz, jest pierwszym z serii postów atakujących ten temat. Postanowiłem dalej dłubać w tej sprawie, bo … trudno mi wyobrazić sobie coś ważniejszego niż powietrze, którym oddychamy. Nie chodzi mi tutaj o sianie paniki, lobbowanie polityczne czy obwinianie kogokolwiek. Bardziej interesuje mnie, co będę w stanie powiedzieć na ten temat, bazując na danych. Nie mam tutaj żadnego grafiku ani planu. Do tematu będę wracał, jeśli wpadnę na jakiś ciekawy pomysł, pojawią się jakieś nowe ciekawe dane albo znajdę jakiś błąd we wcześniejszym kodzie i będę chciał go poprawić. Mam nadzieję, że prezentowane przeze mnie sposoby pracy z danymi zaciekawią Cię na tyle drogi Czytelniku, że sam zaczniesz drążyć temat. Zacznijmy więc.

2013

Jakiś czas temu natrafiłem na artykuł na Wikipedii pod tytułem Lista miast w Polsce o najbardziej zanieczyszczonym powietrzu. Motywacją (zapewne) do powstania tego artykułu była publikacja średniej rocznej stanów polutantów PM2.5 i PM10. Wartości te były mierzone w 113 miastach, a na końcu tabela została ułożona według wartości PM2.5. Zalecenie WHO odnośnie do maksymalnego poziomu rocznej średniej PM2.5 wynosiło wtedy 10 μg/m³, a w tabeli wylądowały miejscowości, które przekroczyły ten poziom co najmniej dwukrotnie, czyli osiągnęły poziom 20 μg/m³. Dla mieszkańców tych miejsc w 2013 jest to chyba dość istotna informacja. Ja wtedy mieszkałem we Wrocławiu, wiec czuję się „objęty” tą informacją.

Dzisiaj?

Niestety z nieznanych mi powodów, ranking taki nie istnieje dla kolejnych lat. A mogłoby to być bardzo ciekawe przedsięwzięcie. Można byłoby sprawdzić które miasta się ogarnęły i wychodzą na prostą, a które „jadą po bandzie”. Można by też użyć takiej tabelki, żeby sobie naocznie określić czy miejsce, do którego się przeprowadzamy, nas nie zabije.

Faktem jest, że mamy dostęp do pełnych danych z pomiarów oferowanych przez GIOŚ od 2000 roku. Więc przy odrobinie wprawy możemy sobie sami skonstruować podobną tabelkę dla poszczególnych lat. Nie uda nam się jej jednak w pełni odtworzyć, bo nie wszystkie miasta mają stacje GIOŚ, które mają czujniki PM2.5 i PM10. No i właśnie tym problemem chciałbym się zająć w pierwszym poście z tego cyklu.

Idea

Idea, którą chciałbym zrealizować jak najlepiej, jest następująca:

• namierzenie najbliższego dla danego miasta czujnika danego polutanta
• pobranie aktualnych danych z czujnika
• wyliczenie średniej
• ustalenie rankingu miast (uwzględniającego te same miasta co w rankingu z 2013 roku)

Do zrealizowania tej idei użyję Pythona i API dostarczonego przez GIOŚ.

Pierwsze kroki

Nie będę tutaj omawiał każdego kawałka kodu. Przede wszystkim dlatego, że jest to kod napisany przeze mnie na kolanie na potrzebę sprawdzenia, czy i jak można to zrobić ;). Będę zmieniał go jeszcze wielokrotnie, więc nie przywiązuję się do niego za bardzo. Ma też już na tym etapie kilka słabych fragmentów, więc nie miałoby to za dużo sensu.

Ok, to od czego więc możemy zacząć? Przydałaby nam się lista miast, którymi jesteśmy zainteresowani. Możemy je sobie skopiować ze wspomnianego artykułu na Wikipedii. Skoro mamy już miasta, to musimy namierzyć najbliższy czujnik dla każdego polutanta. No i tutaj pojawia się największe wyzwanie.

API GIOŚ

API GIOŚ udostępnia nam cztery zapytania, z których możemy swobodnie korzystać:

• Stacje pomiarowe
• Stanowiska pomiarowe
• Dane pomiarowe
• Indeks jakości powietrza

To, co nas najbardziej interesuje to Dane pomiarowe. Mamy tam informacje o polutancie, godzinie i średniej wartości pomiaru dla tej godziny. Bardzo fajna sprawa. Jedynym problemem jest to, że dane te mamy dla trzech ostatnich dni kalendarzowych. Czyli jeżeli zadajemy zapytanie np. 24 listopada, to dane będziemy mieć od 22 listopada od godziny 1-szej w nocy. Nie jest to idealna sytuacja, ale da się z tym żyć. Gdy chcemy uzyskać te dane, to musimy wykonać zapytanie zawierające sensorId.

Skąd wziąć sensorId? Uzyskamy je z odpowiedzi z zapytania o Stanowisko pomiarowe. W takiej odpowiedzi mamy sensorId dla każdego polutanta mierzonego w danej stacji. Czyli jeśli namierzymy stację, to otrzymamy informację o wszystkich sensorach w niej zamontowanych. Dlaczego napisałem „namierzymy”? Bo żeby uzyskać odpowiedź, musimy w zapytaniu podać stationId. Którego, na razie nie znamy.

Jak więc poznać stationId? Ha, tego dowiemy się z zapytania o Stacje pomiarowe. Zapytanie to da nam listę wszystkich stacji pomiarowych, z których dane udostępnia GIOŚ. Mamy tam współrzędne geograficzne stacji, jej stationId, i informacje o położeniu administracyjnym stacji. W ten sposób zbieramy sobie po kolei wszystkie kawałki naszej układanki i możemy przystępować do namierzania najbliższego czujnika.

Najbliższy czujnik

Żeby wyznaczyć najbliższy czujnik dla interesującego nas punktu, wypadałoby znać współrzędne tego punktu i współrzędne wszystkich czujników. Najlepiej byłoby znać te współrzędne w trzech wymiarach (długość i szerokość geograficzna oraz wysokość nad poziomem morza). Niestety dane GIOŚ nie zawierają wysokości nad poziomem morza, musimy więc ograniczyć się do dwóch wymiarów.

Do tej pory mieliśmy jednak tylko listę miast, jeśli chcemy więc poznać, gdzie jest najbliższy czujnik dla miasta Wrocław, musimy przekuć tę nazwę na współrzędne geograficzne. Użyłem do tego modułu GeoPy i funkcji geocode z domyślnymi parametrami.

Stworzyłem więc pętlę, która przechodzi po miastach i przekuwa ich nazwy na współrzędne geograficzne. Super, mamy więc już współrzędne wszystkich punktów, dla których danych będziemy poszukiwać. Teraz pora na czujniki.

Jeśli czytałeś uważnie opis API, zauważyłeś pewnie, że czujniki same w sobie nie mają współrzędnych. Trzeba je w sprytny sposób spiąć z listą stacji. Zrobiłem to w następujący sposób:

1. Pobrałem listę wszystkich stacji pomiarowych i przerobiłem ją na ramkę danych. Mam w niej kolumny ze współrzędnymi stacji.
2. Przechodzę po każdej stacji i pobieram listę czujników, które są tam dostępne. Na tej liście dokonuję operacji one hot encode i doklejam ją do ramki danych z punktu 1.
3. Przechodzę przez listę miast. Dla każdego miasta wyliczam odległość do każdej stacji.
4. W powyższej pętli przechodzę przez listę polutantów i wybieram najbliższą stację która, posiada taki czujnik.

W ten sposób każde miasto dostaje dane dla każdego z siedmiu polutantów z najbliższych możliwych czujników.

Dane

Wszystko, co udało mi się uzyskać w powyższym procesie, zapisuję na dysku w surowych plikach JSON. Oprócz tego, od razu wyliczam też średnią z zebranych pomiarów i wrzucam w ramkę danych w odpowiednią kolumnę. Najtrudniejsze za nami.

Gdy już mamy średnie z najbliższych sensorów dla każdego miasta możemy przystąpić do posortowania tych wartości. Sortowanie, na razie robię poprzez sortowanie kolejnych kolumn od wartości największej do najmniejszej. Biorę więc ramkę danych i sortuję ją według wartości w kolumnie PM2.5 (chyba najważniejszy polutant). Następnie wszystkie miejsca, gdzie jest remis, sortuję według PM10. Wszystkie remisy, które wciąż są obecne, sortuję kolejno po C6H6, NO2, SO2, O3 i CO. Gdy już mam posortowaną ramkę danych, wyliczam różnicę w położeniu względem tabelki z 2013. Dodatkowo jeszcze koloruję poszczególne komórki w zależności, do jakiej z 5 kategorii zdefiniowanych przez GIOŚ wpadają. Przykładowy wynik poniżej.

Pełny kod tworzący tabelę taką jak powyżej, znajduje się tutaj.

Dyskusja

1. Sposób pobierania danych uzależniłem od ilości miast na liście. Jeśli więc będzie tam mało miast, to ostatecznie odpytam małą liczbę sensorów o dane. Idea jest taka, że każdy punkt powinien otrzymać dane z 7 sensorów. Dla trzech miast będzie to więc maksymalnie 21 zapytań do API Dane pomiarowe. Maksymalnie, bo niektóre miasta mogą mieć te same najbliższe czujniki. Problem zaczyna się przy dużej ilości miast. Aktualnie GIOŚ oferuje dostęp do 648 czujników, co przy optymalnym podziale (tyle samo czujników dla każdego polutanta) daje 92 unikatowe stacje pomiarowe. Jeśli więc liczba miast będzie większa niż 92, to na pewno niektóre czujniki będą się powtarzać. W rzeczywistości powtórzenie czujników (czyli powtórzenie zapytań do API) nastąpi wcześniej. Z racji, że na liście mam 113 miast, mam sytuację, w której odpytuję niektóre czujniki wielokrotnie. Marnotrawstwo. Optymalnie byłoby zapisywać danie, po każdym zapytaniu i przy kolejnym sprawdzać, czy już ich nie mamy. Jeśli je mamy to nie łączyć się do API. Można by też się pokusić o pobranie wszystkich danych z czujników i później odpytywać lokalny zasób.
2. Niektóre komórki mają wartości nan. Wynika to z faktu, że niektóre czujniki mogą być uszkodzone, ale dalej są wystawione w API. Aktualnie nie sprawdzam, czy czujnik dostarcza w ogóle jakieś wartości, stąd też niekiedy w tabeli lądują wartości nan.
3. Model. Aktualny model szacowania zanieczyszczenia w danym punkcie jest trywialny: weź wartość z najbliższego geograficznie sensora. Jest to ogromna naiwność i prawdopodobnie najsłabszy element całego skryptu. Należałoby tutaj skonstruować lepszy model przestrzenny. Prawdopodobnie dla każdego polutanta. Najwięcej czujników mamy dla NO2 – 135. Nie wiem, ile z nich jest aktualnie sprawne, ale tutaj można by się pokusić o jakieś modelowanie. Najsłabiej wypada PM2.5 – 46 (sic!). Tutaj za dużo nie ugramy. Pozostaje jeszcze w ogóle pytanie, czy da się coś sensownego w ogóle.
4. Ranking. Aktualnie buduję ranking, sortując ramkę danych, tak jak opisałem powyżej. A może nie ma to sensu? No bo możemy sobie wyobrazić miasto, które ma niski poziom PM2.5, a kosmiczny NO2. W rankingu na najgorsze miejsce będzie więc nisko, a ogólnie może tam się nie dać oddychać. Zastanawiam się więc, czy nie przeskalować każdego z polutantów tak, żeby wartość graniczna dla klasy „bardzo zły” dawała 1. Wtedy każde miejsce miałoby serię wartości w odniesieniu do 1 dla każdego polutanta. I wtedy można by ranking skonstruować na dwa sposoby. Posortować taką listę malejąco i rozbić na kolumny. Wtedy w pierwszej kolumnie byłyby najgorsze czynniki (w tej samej skali), w drugiej byłby by drugie najgorsze i tak dalej. Wtedy nie wybieralibyśmy arbitralnie najgorszego czynnika (np. PM2.5) ale po prostu brany byłby ten, który ma największą wartość. Co jednak z sytuacją, gdy mamy miejsce, gdzie pierwszy czujnik będzie miał wartość 2 (2 razy wartość graniczna „bardzo zły”) a pozostałe czujniki będą miały wartości 0.1. Będziemy mieć też drugie miejsce, które wszystkie czujniki będzie miało na poziomie około 1. Pierwsze miejsce wygra w rankingu. Natomiast w drugim miejscu może znaleźć się całkiem nisko. A czy powinno? Może w pierwszy miejscu mamy awarię filtrów w fabryce, a drugie miejsce znajduje się niedaleko miejsca jakiejś katastrofy naturalnej? Zastanawiam się więc, czy po prostu nie sumować tych wartości. Wtedy miejsca, które wszystkie czynniki mają na wysokim poziomie, będą mogły w rankingu wylądować wyżej niż miejsca, które mają jeden czynnik na bardzo wysokim poziomie. Ale to już zdecydowanie kwestia do dyskusji.

TODO

• Stworzyć plik ze słownikiem „lokacja”: współrzędne. Dokonywać konwersji przy pomocy API, tylko jeśli pojawi się nowa lokacja.
• Sprawdzić, czy mamy dane z danej godziny, jeśli tak, to wtedy pobrać wszystkie.
• Określić jakie wartości, jakiej substancji są bardziej groźne. Skalowanie typu (górna granica „Bardzo dobry”) / (dolna granica „Bardzo zły”) daje różne wyniki. Np. O3 i SO2 – wartości < 71 i < 51 wpadają w przedział „Bardzo dobry” jednakże pierwsza to ok. 29% a druga to ok. 10% przedziału „Bardzo źle”. Jak więc je posortować?
• Lepszy model przestrzenny? Aktualny to: dane z najbliższego czujnika.
• Stworzenie największych możliwych okręgów na mapie Polski bez czujników danego rodzaju. Środek takiego okręgu będzie aż się prosił w GIOŚ o czujnik ;).
• Przejść z plików JSON na bazę danych. SQLite?
• Zebrać to wszystko i stworzyć moduł.
• Prognozowanie.
• Ranking wszystkich 930 miast.

Podsumowanie

Jak widać, dane z czujników powietrza mogą dostarczyć dość dużo pola do analiz. Mamy paczki z danymi archiwalnymi i dane aktualne prosto z API. Oprócz modelowania przestrzennego można też modelować prognozy szeregów czasowych dla poszczególnych czujników. Fakt, czujników tych jest nieco mało (~600 sprawnych), ale przypuszczam, że wraz z rozwojem technologii będą coraz tańsze (więc GIOŚ będzie mógł ich rozmieszczać coraz więcej). A i świadomość ludzi jest coraz większa, więc politycznie też pojawia się coraz większa presja na poprawę jakości powietrza. A nie da się czegoś poprawić, jeśli się tego nie mierzy.