Czy chcielibyście zobaczyć mapę Polski, na której zaznaczony jest każdy Polak? Wizualizację ponad 38 milionów osób rozrzuconych mniej lub bardziej losowo na obszarze ponad 300 tysięcy km^2? Jeżeli w tej chwili w waszych głowach zapala się pytanie 'po co?’, nie czytajcie dalej. Jeżeli zaś już widzicie taką mapę oczami wyobraźni, poniższy wpis bardzo Wam się spodoba.

Celem tego projektu było pokazanie w interesujący sposób informacji z Narodowego Spisu Powszechnego 2011. Jedną z inspiracji była interaktywna mapa przedstawiająca spis powszechny w Stanach Zjednoczonych o której przeczytać można tutaj. W ramach dzisiaj opisywanego projektu opracowano mapę Polski składającą się z ponad 38 milionów punktów, każdy punkt odpowiadający jednej osobie, rozmieszczenie punktów odpowiadające rzeczywistej gęstości zaludnienia.

Ten wpis jak i cały projekt został wykonany przez Pawła Wiechuckiego w ramach wolontariatu dla naszej fundacji.
Kiedyś podczas jakiejś prezentacji usłyszałem, że organizacje non-profit to ciekawe miejsca, ponieważ przyciągają osoby, które za darmo tworzą rzeczy bezcenne, dlatego że chcą by te rzeczy powstały. Projekt Pawła Wiechuckiego jest świetnym przykładem takiej aktywności. Nie było łatwo, trzeba było pokonać wiele trudności, których nie widać przy bardziej standardowych zastosowaniach [np. R nie potrafi wygenerować w formacie wektorowym kwadratu o boku mniejszym niż pół punktu tzn. 1/144 cala, nie można też wygenerować dowolnie dużego rysunku, a przynajmniej 100GB RAM to za mało, sama mapa to połowa sukcesu, bufor drukarki może być niewystarczający by taką mapę wydrukować itp], ale się udało.

Przy okazji dziękuję też firmie Iqor Polska za udostępnienie plotera drukującego na 42 calowej rolce, dzięki temu udało się wydrukować mapę na połączonych arkuszach o łącznej powierzchni ponad 6 metrów kwadratowych, powiesić na ścianie i zrobić zdjęcie młodej wskazującej kropkę przedstawiającą babcię.

Nieszablonowa mapa Polski

Paweł Wiechucki

Przepis na stworzenie wizualizacji danych o gęstości zaludnienia na bazie danych o spisie powszechnym 2011 jest z pozoru bardzo prosty. Potrzebne do tego są tylko dwa składniki. Pierwszy z nich to raport o ludności w gminach udostępniony przez GUS jako część wyników spisu ludności i mieszkań 2011. Drugi – mapa administracyjna Polski w podziale na gminy. Niestety, tak jak to zwykle bywa, również w tym przypadku nie mamy do czynienia z przepisem na wizualizację z torebki. Nie wystarczy wymieszać gotowych składników, żeby uzyskać efekt. Nawet jeśli pod ręką mamy statystyczny mikse-R i tak będziemy musieli się trochę ubrudzić, bo zanim użyjemy go do stworzenia naszej wizualizacji, dane musimy wcześniej odpowiednio wycisnąć, poszatkować i zmielić. Problem w tym, że wyciskarkę, noże i młynek musimy sobie w dużej mierze sami zbudować. Ale zacznijmy od początku.

Dane o ludności w gminach można pobrać ze strony internetowej GUS [link]. Dostępne są one w postaci pliku xls. Dane zebrane z poszczególnych województw, zostały umieszczone w osobnych arkuszach. Każdy arkusz ma sekcję nagłówka. Dane z gmin są pogrupowane pod rekordami zawierającymi zagregowane dane z odpowiadających im powiatów. Grupy porozdzielane są pustymi liniami.
Napisałem mały skrypt w Perlu, który pozbywa się wszystkich tych ozdobników, ułatwiających człowiekowi czytanie i tworzy przyjazny maszynie plik csv o kolumnach: Wojewodztwo, Powiat, Gmina, Symbol terytorialny, Ogolem, Miasto i Wies. Dodatkowo, znaki myślnika, wskazujące na brak danych w danej kategorii (Miasto, Wies), zastąpiłem zerami.

Dużo trudniej jest zdobyć i odpowiednio przygotować mapę administracyjną Polski w podziale na gminy w formacie Shapefile (.shp). Można się w tym celu odwołać do zasobów oferowanych przez projekt OpenStreetMap [link], jednak dane o granicach gmin udostępniane w jego ramach nie są jeszcze kompletne [link].

Rzetelnym źródłem danych, na których zdobyciu mi zależało, okazał się być Geoportal [link]. Udostępnia on zbiór usług zgodnych ze standardami OGC, wśród których znajduje się usługa Polskiego Rejestru Granic. Wyszczególnione usługi są usługami WMS (Web Map Service), które w odpowiedzi na zadane parametry zwracają rastrowy obraz mapy. Niestety w tej formie dane o granicach gmin byłyby dla nas kompletnie bezużyteczne. Jednak tajemnicą poliszynela jest, że te same dane udostępniane są przez analogiczne usługi WFS (Web Feature Service) w formacie GML (Geography Markup Language). Dzięki temu, budując odpowiednie żądania, możemy pobrać z serwera informacje o granicach administracyjnych, które następnie możemy wykorzystać do stworzenia naszej mapy.

W tym miejscu należy zaznaczyć, że zgodnie z regulaminem Geoportalu zasoby, które udostępnia nieodpłatnie, można wykorzystywać jedynie w celach „niezarobkowych i niekomercyjnych”. Nie można też ich rozpowszechniać w żadnej postaci. Dlatego nie ma przeciwwskazań, żeby skorzystać z nich w celach hobbystycznych, jednak już np. ich publikacja na blogu lub wykorzystanie w prezentacji bez odpowiedniej zgody było by naruszeniem regulaminu. W szczególności, danymi z Polskiego Rejestru Granic dysponuje Centralny Ośrodek Dokumentacji Geodezyjnej i Kartograficznej (CODGiK) [link].

Nawiasem mówiąc, w trakcie kiedy pracowałem nad stworzeniem wizualizacji, pojawiła się informacja, że polska społeczność zgromadzona wokół OpenStreetMap walczy o wolny dostęp do danych z Geoportalu [link].

Jeśli w przeglądarce wpiszemy następujący adres
http://sdi.geoportal.gov.pl/wfs_prg/wfservice.aspx?request=getcapabilities&version=1.1.0&service=wfs
zobaczymy opis usługi WFS Państwowego Rejestru Granic. Możemy z niego wyczytać, że udostępnia ona informacje o województwach, powiatach i gminach, a także jakie są kody EPSG układów współrzędnych w których udostępniane są dane. Jeśli zaimportujemy dane w domyślnym układzie współrzędnych opisywanym kodem EPSG:4326, mapy które uzyskamy będą zdeformowane – zamiast kształtu Polski, który znamy i który prawie można wpisać w kwadrat, otrzymamy Polskę w postaci szerokiego i niskiego placka. Żeby tego uniknąć albo musimy pobrać dane w jednym z dwóch alternatywnych układów współrzędnych (EPSG:2180, EPSG:3785) albo musimy zmienić układ domyślny na inny, korzystając z dodatkowego narzędzia np. Quantum GIS.

Żeby uzyskać dane w formacie GML w domyślnym układzie współrzędnych należy w przeglądarce wpisać:
dla województw: http://sdi.geoportal.gov.pl/wfs_prg/wfservice.aspx?request=getfeature&version=1.1.0&service=wfs&typename=WOJEWODZTWA
dla powiatów: http://sdi.geoportal.gov.pl/wfs_prg/wfservice.aspx?request=getfeature&version=1.1.0&service=wfs&typename=POWIATY
dla gmin: http://sdi.geoportal.gov.pl/wfs_prg/wfservice.aspx?request=getfeature&version=1.1.0&service=wfs&typename=GMINY

Zmianę układu współrzędnych uzyskujemy przez dopisanie &srsName=EPSG:2180 lub &srsName=EPSG:3785 na końcu zapytania.

Wyniki, które pojawią się nam w przeglądarce, należy zapisać na dysku jako pliki o rozszerzeniu GML. Teraz możemy przekonwertować je do formatu Shapefile, który posłuży nam do stworzenia wizualizacji danych w R. Mnie osobiście pomogły w tym wskazówki, które znalazłem na blogu Geoinformacja Pawła Gajewskiego [link]. Sama konwersja nie nastręcza problemów, trzeba jednak wiedzieć, że dane pobierane z Geoportalu mają zamienione miejscami współrzędne x i y. W swoim wpisie Pan Gajewski tłumaczy, jak sobie z tym poradzić.

Zanim będziemy mogli pokazać dane na mapie, do rozwiązania pozostał nam jeszcze jeden problem: jak powiązać rekordy danych z GUS z obszarami w pliku Shapefile, które odpowiadają konkretnym gminom. Jeśli jeszcze raz podejrzymy zawartość pliku GML dla gmin zobaczymy, że każda gmina opisana jest przez nazwę i identyfikator. Identyfikatory niestety nie odpowiadają symbolom terytorialnym z bazy TERYT, które można znaleźć w bazie danych z GUS. Pozostaje nam próba powiązania odpowiednich rekordów po nazwach.

Dane z GUS opatrzone są przedrostkami „M.”, „M-W.”, „G.” lub „m.”. Jeden rekord często zawiera informację o ludności i miasta i wsi pod wspólną nazwą np.

Wojewodztwo,Powiat,Gminy,Symbol terytorialny,Ogolem,Miasto,Wies
Wojewodztwo Dolnośląskie,Powiat górowski,M-W.Góra,0204013,20963,12575,8388

Dane mapowe z Geoportalu nie są opisane kodami terytowymi a tereny wiejskie i miejskie w sytuacjach, jak powyżej mają swoje własne obszary. Nazwy gmin są opisane inaczej niż przez GUS, mianowicie albo wcale nie mają opisu albo widnieje przy nich dopisek „miasto”, „gm. miejska” lub „obszar wiejski”. Nazwy dużych miast poprzedzone są przedrostkiem „M.”.

Jeśli występuje jeden rekord mapowy odpowiadający z nazwy jednemu rekordowi z GUS, wtedy nie ma wątpliwości, jak je powiązać. Często jednak jest kilka rekordów z GUS, które można powiązać z jednym rekordem mapowym. Skłoniło mnie to do wykorzystania informacji o zawieraniu się gmin w województwach. Dla danych z GUS ta informacja jest dana, jeśli chodzi o dane z Geoportalu takie mapowanie trzeba sobie stworzyć. Wykorzystałem do tego skrypt R, który dla każdego obszaru gminy wyliczał centroidę a następnie sprawdzał jej zawieranie w obszarach województw. Tu pojawia się dodatkowo problem identycznych nazw gmin w tym samym województwie. Można by go rozwiązać w analogiczny sposób na poziomie powiatu, ale niestety Geoportal nie podaje nazw powiatów, kiedy korzysta się z ogólnodostępnych usług, więc informacja ta byłaby bezużyteczna.

Inny pomocniczy skrypt wykorzystałem do pobrania z pliku Shapefile dla gmin powiązania pomiędzy ich identyfikatorami a nazwą. Było to łatwiejsze rozwiązanie niż parsowanie w tym celu pliku GML.

Uzyskane w wyniku działania opisanych powyżej skryptów pliki csv z mapowaniami oraz plik ze sparsowanymi danymi z GUS zasiliły skrypt Perl, który posłużył mi do znalezienia ile osób zamieszkuje poszczególne obszary wskazane na mapie gmin. Tyle właśnie punktów należało wygenerować losowo dla każdego z obszarów przy generowaniu wizualizacji. Dla kilkunastu gmin nie udało mi się znaleźć jednoznacznej odpowiedzi na pytanie, które rekordy powiązać. Tu nie zostało mi nic innego, tylko rysowanie konturów gmin na mapie w R a potem szukanie odpowiadającej gminy w Google Maps. Tam też znajdowałem informację o powiecie, w którym leży dana gmina. Znając nazwę gminy i powiatu uzupełniałem informację o liczebności korzystając z danych z GUS. Teraz możemy wreszcie przejść do tworzenia wizualizacji.

Do tego celu można skorzystać z różnych narzędzi, jednak wszystkie przygotowania, które do tej pory opisywałem, miały na celu takie wstępne przetworzenie danych, żeby łatwo można było z nich skorzystać do stworzenia mapy w R. Pomijając problem spłaszczonej Polski, o którym pisałem na początku, napotkałem jeszcze dwa związane bezpośrednio z tworzeniem grafiki.

Po pierwsze, okazało się, że funkcja runifpoint, z której korzystałem, żeby wygenerować losowe punkty rozłożone równomiernie na obszarze każdej z gmin, ma problemy kiedy ilość punktów przekracza milion (wielkość wyznaczona empirycznie). Z tego powodu, jeśli liczba punktów była większa, generowałem po 500 tysięcy punktów dopóki pozostałe nie mogły być wygenerowane bezproblemowo. Bez zastosowania takiej procedury, nie można by wygenerować punktów np. dla Warszawy.

Po drugie, bez zastosowania półprzezroczystości i zmniejszenia rozmiaru rysowanego punktu, mapa była wielką, czarną plamą.

Ale nawet po rozwiązaniu powyższych problemów generowany plik pdf był na tyle duży, że nie sposób było go otworzyć. Dopiero konwersja do formatu postscript pozwoliła cieszyć się widokiem przygotowanej wizualizacji.

Źródła

Skrypty w językach R, perl i python, wykorzystywane do opracowania map pobrać można z tego linku ProjektNSP2011.zip. Z uwagi na licencje, pliki z serwisu geoportal należy pobrać oddzielnie/samodzielnie.
Trwają rozmowy dotyczące odpłatności za licencje na zaprezentowanie map, na razie więc w powyższym pliku nie ma wygenerowanych map. Można te mapy legalnie na użytek własny wygenerować z załączonych skryptów. Po uregulowaniu kwestii licencyjnych opracowane mapy również zostaną dodane do powyższego pliku.