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:

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:

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:

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:

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.

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):

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.