Jakość powietrza w Polsce #3 – gdzie brakuje nam czujników?

Od jakiegoś czasu potrafimy się dobrać do danych z czujników powietrza w Polsce (sprawdź artykuły spod tagu powietrze). Potrafimy też określić, ile ich w ogóle jest i gdzie są rozmieszczone. Spróbujmy więc może określić, gdzie tych czujników jest najmniej. To znaczy, spróbujmy wyznaczyć największą czujnikową białą plamę na mapie Polski.

Nie jest to zagadnienie stricte związane z data science i uczeniem maszynowym. Jest to natomiast zagadnienie algorytmiczne. A jeżeli siedzimy w data science, to prędzej lub później i tak natrafimy na jakieś problemy, do których rozwiązania droga będzie prowadzić przez różne pomocnicze algorytmy.

Problem

Zastanówmy się, jak dokładnie wygląda nasz problem. Chcemy wyznaczyć obszar na płaszczyźnie, korzystając z punktów (na tej samej płaszczyźnie). Jak to będzie wyglądać w praktyce. Weźmy sobie garść punktów:

Punkty
Punkty

Na razie, szału nie ma. Mamy jakieś kropki. Musimy narysować jakiś kształt, który nie będzie ich zawierał. Czy naprawdę jednak potrzebujemy precyzyjny kształt? Może wystarczy nam okrąg? Aktualnie sądzę, że okrąg będzie wystarczającym rozwiązaniem. Największy możliwy okrąg przyklejony do jakiejś kropki. Na przykład:

Punkty i okrąg z jednym zaczepieniem
Punkty i okrąg z jednym zaczepieniem

No ok, jest to jakiś obszar faktycznie bez punktów. Problem jest tylko taki, że okrąg ten, możemy przykleić do każdego z najbardziej zewnętrznych punktów. Co gorsza, jeśli będziemy szukali największego takiego okręgu, to okaże się, że dla każdego punktu zewnętrznego możemy stworzyć okrąg o nieskończonym promieniu. Będziemy więc mieć co najmniej kilka takich okręgów o nieskończonych promieniach i każdy z nich będzie równie dobry (a może zły?). Jest to więc trochę lipa. Poszukajmy więc dalej:

Punkty i okrąg z dwoma zaczepieniami
Punkty i okrąg z dwoma zaczepieniami

Tutaj sytuacja wydaje się już być trochę lepsza, ale czy na pewno? No bo nie wciśniemy nieskończenie wielkiego okręgu tak, żeby dotykał te dwa punkty. Nie? A co jeśli będziemy mieć taki wielki okrąg, że będzie nam się wydawało, że te dwa punkty leżą na prostej linii? Coś w stylu niemożliwości zauważenia krzywizny ziemi na boisku piłkarskim (pomijając fakt, że ktoś mógłby takie boisko ekstra wyrównywać). Czyli, że okrąg, który ma dwa punkty zaczepienia, też niekoniecznie nam się tutaj przyda. Dokonajmy więc kolejnego kroku i poszukajmy okręgów z trzema punktami zaczepienia:

Punkty i okrąg z trzema zaczepieniami
Punkty i okrąg z trzema zaczepieniami

I tutaj chyba wszystko wygląda ok. Nawet jeśli nie „oszukamy” tak jak ja tutaj wybierając punkty „wewnątrz” naszej przestrzeni, to i tak, jeśli będziemy budować okręgi na „trójkach” to nie będą one nieskończenie wielkie.

I to jest mniej więcej sendo naszego problemu – znaleźć największy okrąg, który możemy opisać na trzech punktach. Dodatkowym ograniczeniem jest to, że nie może być żadnych innych punktów wewnątrz naszego okręgu. Musimy po prostu znaleźć największą białą i okrągłą kartkę, którą możemy położyć na punktach, tak żeby żadnego nie przykryć. Tak mniej więcej sformułowałbym ten problem prostymi słowami.

Jak więc go rozwiązać?

Diagram Woronoja

Zmieńmy na chwilę temat (pozornie). Kojarzysz może taki twór jak diagram Woronoja? Jest to diagram, który powstał jako podział płaszczyzny na komórki. Komórka taka reprezentuje obszar, w którym wszystkie możliwe punkty będą miały najbliżej do wyznaczonego przez nas punktu. Dla każdego wskazanego punku możliwe będzie stworzenie odpowiedniej komórki.

Diagram Woronoja
CC BY-SA 3.0, Link

Zagadnienie to najłatwiej sobie zobrazować wyobrażając sobie problem urzędu pocztowego. Mamy sobie jakieś miasto, a w tym mieście np. 15 urzędów pocztowych. Jako że urzędy te już tam stoją, nie możemy za bardzo dostosować ich położenia np. względem gęstości zaludnienia. Innymi słowy, nie możemy ich poprzesuwać. To, co możemy jednak zrobić, to upewnić się, że każdy adres jest przyporządkowany do najbliższego urzędu pocztowego. W tym celu wyliczamy sobie odległość każdego adresu do każdego urzędu i odpowiednio je przyporządkowujemy (trochę podobnie jak w k-średnich). Jeżeli oznaczymy sobie te obszary różnymi kolorami, to otrzymamy coś podobnego do diagramu Woronoja. Dlaczego podobnego? Bo nasz diagram będzie miał białe plamy, a w diagramie Woronoja wyznaczamy również granicę pomiędzy obszarami. Najlepiej będzie, jak spojrzysz na obrazek po prawej stronie – jest to pełnoprawny diagram Woronoja.

Po co może być nam taki diagram? Nie szukamy przecież komórek wokół czujników powietrza? No właśnie. Jeżeli przyjrzysz się powyższemu obrazkowi w powiększeniu, zauważysz, że obojętnie jakich by te komórki nie miały kształtów, zawsze jakaś krawędź będzie zaczynać się i kończyć w wierzchołku. A czym jest taki wierzchołek? Wierzchołek ten jest punktem, który ma taką samą odległość do centrum trzech przylegających do niego komórek. Przypomnijmy sobie, co jest w tym centrum? Jest tam czujnik. Wierzchołek taki reprezentuje więc miejsce, które jest równo odległe od trzech najbliższych czujników. Jest więc środkiem okręgu. BUM, czegoś takiego właśnie szukaliśmy.

Największy okrąg

Wystarczy więc namierzyć największy okrąg na powstałym diagramie i mamy już problem z głowy. A jak to ugryźć? Implementacja, z której skorzystałem (scipy.spatial.Voronoi) daje mi tak utworzone wierzchołki. Sprawdzam więc sobie odległość każdego wierzchołka do wszystkich czujników (overkill, ale jeszcze nie wpadłem jak sprawdzić odległość to tych trzech czujników, z których powstał). Następnie dla każdego wierzchołka zapisuję sobie najmniejszą zaobserwowaną odległość. Czyli okrąg, do którego nie wpada żaden czujnik. Jak już przejdę przez wszystkie wierzchołki, to wybieram największą odległość. Problem rozwiązany.

Wizualizacja

Czymże byłaby ta analiza bez chociażby jednej wizualizacji. Okazuje się, że mogę przedstawić ich tutaj aż siedem. Po jednej dla każdego rodzaju czujników. Zobaczmy więc, co uzyskałem (pod każdą lokalizacją jest teoretycznie bardziej ludzki adres środka okręgu):

PM2.5 - S11 (koncepcja), Sieranie, gmina Świeszyno, powiat koszaliński, zachodniopomorskie, 76-022, RP
PM2.5 – S11 (koncepcja), Sieranie, gmina Świeszyno, powiat koszaliński, zachodniopomorskie, 76-022, RP
PM10 - Dąbek, gmina Stupsk, powiat mławski, mazowieckie, 06-561, RP
PM10 – Dąbek, gmina Stupsk, powiat mławski, mazowieckie, 06-561, RP
C6H6 - Psie Głowy, gmina Czaplinek, powiat drawski, zachodniopomorskie, 78-550, RP
C6H6 – Psie Głowy, gmina Czaplinek, powiat drawski, zachodniopomorskie, 78-550, RP
NO2 - 33, Czarkówka Duża, gmina Perlejewo, powiat siemiatycki, podlaskie, RP
NO2 – 33, Czarkówka Duża, gmina Perlejewo, powiat siemiatycki, podlaskie, RP
SO2 - 33, Czarkówka Duża, gmina Perlejewo, powiat siemiatycki, podlaskie, RP
SO2 – 33, Czarkówka Duża, gmina Perlejewo, powiat siemiatycki, podlaskie, RP
O3 - Psie Głowy, gmina Czaplinek, powiat drawski, zachodniopomorskie, 78-550, RP
O3 – Psie Głowy, gmina Czaplinek, powiat drawski, zachodniopomorskie, 78-550, RP
CO - Świelubie, gmina Dygowo, powiat kołobrzeski, zachodniopomorskie, 78-113, RP
CO – Świelubie, gmina Dygowo, powiat kołobrzeski, zachodniopomorskie, 78-113, RP

Jak widać, kilka regionów się powtarza – głównie województwo zachodniopomorskie. Nie wiem, z czego to wynika, ale na miejscu mieszkańców Kołobrzegu i Koszalina bym się tym mocniej zainteresował.

Problem

Na części z powyższych wizualizacji możemy zauważyć, że narysowane okręgi tak średnio stykają się z czujnikami. A przecież takie było założenie rozwiązania naszego problemu. Aktualnie mam dwie hipotezy odnośnie do przyczyn tej nieścisłości:

  1. Popełniłem gdzieś błąd w kodzie. Cóż, kod, do którego link znajdziesz na końcu, jest baaardzo roboczy. Może gdzieś tam się pomyliłem i widzimy właśnie tego efekty.
  2. Któryś z użytych przeze mnie modułów źle obsługuje współrzędne geograficzne. W tym miejscu mam trzy podejrzenia: scipy.spatial.Voronoi, geopy.distance.vincenty i folium.vector_layers.Circle. Przy czym drugi i trzeci moduł są stworzone z myślą o współrzędnych geograficznych, o ile więc nie są skopane, powinny dobrze ogarnąć temat. Pierwszy z kolei raczej oczekuje tak samo wyskalowanych liniowych wartości, podejrzewam więc, że tu jest pies pogrzebany ;(.

Wierzchołki poza krajem

Okazuje się, że obiekt Voronoi wyznacza też wierzchołki, które znajdują się poza obszarem punktów. Może więc się zdarzyć, że największy okrąg wyjdzie nam np. w Czechach. Pierwszym pomysłem na rozwiązanie tego problemu było ograniczenie listy wierzchołków to tych mieszczących się w kwadracie wyznaczonym przez maksymalne położenie czujników w kierunkach północ, południe, wschód, zachód. Działało to fajnie przez chwilę, bo i tak Czechy się pojawiły w przypadku południowo-zachodniej części kraju. Później miałem pomysł na stworzenie najbardziej wypukłego kształtu na bazie czujników (w miarę prosta sprawa), ale wtedy też mogłoby się okazać, że wyrzucam dobre punkty (np. okolice Koszalina przy polutancie PM2.5). Teraz więc po prostu robię reverse geocode i sprawdzam, czy dostaję tam [„country”] == „RP”. Tak jest prościej.

Wnioski

Przy pomocy diagramów Woronoja możemy całkiem sensownie sobie powydzielać obszary wokół jakiś interesujących nas punktów. Tutaj na przykład możemy znaleźć zastosowania w naukach przyrodniczych, medycynie, inżynierii i informatyce. Jeśli już chcemy używać tych diagramów, to powinniśmy pomyśleć o lepszej reprezentacji współrzędnych punktów geograficznych (może coś z tego albo tego wątku?), albo o bardziej dostosowanym algorytmie. Jeśli tego nie zrobimy, to nieliniowość reprezentacji długość-szerokość geograficzna mogą nam popsuć szyki. No, a jeśli poszukujemy największych okręgów niezawierających żadnych punktów, to ten sposób będzie chyba najlepszym rozwiązaniem.

Pełny kod użyty w artykule znajduje się tutaj.

6 myśli na temat “Jakość powietrza w Polsce #3 – gdzie brakuje nam czujników?

    1. Dzięki!
      Diagramy Woronoja to fajna „zabawka” (tak o nich pomyślałem kilka lat temu gdy na nie trafiłem), która pojawia się w różnych kontekstach. Warto mieć ją w swoim arsenale ;).

  1. Z tego co piszesz wynika, że Diagramy Woronoja dają ciekawe możliwości zastosowania. Wydaje mi się jednak, że nie w przypadku rozmieszczenia czujników zanieczyszczenia. O ich rozmieszczeniu decydują zupełnie inne czynniki np. obecność źródła emisji.

    1. Hej!
      Tak, tutaj zdecydowanie otwiera się pole do dyskusji. Jeśli się przyjrzymy tym „białym plamom” to niektóre wypadają na największym poligonie w europie. Mogło by się więc wydawać, że jest to marnotrawstwo, żeby montować stację czujników w środku takiego miejsca. I pewnie tak jest.

      Jeśli stać nas na 50-100 czujników i mamy zaplecze naukowe które wskazuje nam arbitralnie gdzie jest największa lipa to faktycznie od tych miejsc powinniśmy zacząć. Jest to eksploatacja naszej wiedzy domenowej. Ale jeśli nie dokonamy eksploracji w pozornie bezsensowne miejsca, to możemy przeoczyć coś istotnego – wpływ manewrów wojskowych na jakość powietrza (w czołgach chyba nikt się nie szczypie z filtrami spalania, zmasowany ostrzał altymetrii rakietowej też może być ciekawy do zbadania), co się dzieje przy płonących wysypiskach śmieci?

      Fakt, nie używał bym diagramów do rozplanowania położenia wszystkich przyszłych stacji pomiarowych – wtedy eksploracji było by za dużo, a „znane” miejsca mogły by zostać zignorowane. Ale jeden czujnik rocznie w środku tak wyliczonej białej plamy byłby chyba do zrobienia? Takie moje rozmyślania 😉

  2. Ciekawy artykuł. (Pozytywnie 😉 ) nie zgodzę się jednak ze zdaniem: „Nie jest to zagadnienie stricte związane z data science i uczeniem maszynowym”, bo jest – tylko poza głównym nurtem. Wykorzystując podział obszaru na komórki, a następnie wyznaczając średnicę tych komórek, koncepcyjnie wypracowałeś metrykę, która pozwoliła Tobie na klasyfikację obszarów i odpowiedź na pytanie postawione na początku artykułu (na których obszarach gęstość czujników jest zbyt mała?). Taka metryka jest bardzo przydatna w estymacji błędów interpolacji, ale mówi również o potrzebie nakładów inwestycyjnych na czujniki punktowe ALBO o tym, że czujniki punktowe, koniec końców, nie nadają się do pomiarów wielkoskalowych zanieczyszczeń.

    Właściwie problem już rozwiązałeś, ale rzucę Ci jeszcze jedną metryką, która pozwoli określić, czy czujniki są rozmieszczone poprawnie: Morans I. Możesz sprawdzić, czy nie ma zjawiska zbytniego grupowania czujników (co jest nieporządane w tym wypadku).

    Ogólnie obserwuję Twoje wpisy od pewnego czasu i polecam zainteresować się geostatystyką, bo dużo algorytmów i rozwiązań jest tam dostępnych a badasz zjawiska przestrzenne i coś może być dla Ciebie przydatne.

    1. Dziękuję za cenny komentarz!

      Metryki Moran’s I faktycznie nie znałem, dopisuję ją sobie do backlogu tematów do zbadania. O geostatystyce w ogóle nie myślałem jeszcze. Ale zdaje się to mieć sens, szczególnie że jeszcze planuję drążyć temat jakości powietrza. Mam wrażenie, że badam na razie wierzchołek góry lodowej 🙂

Dodaj komentarz

Twój adres email nie zostanie opublikowany. Pola, których wypełnienie jest wymagane, są oznaczone symbolem *