Mapy, mapy raz jeszcze

W jakim województwie leży wskazany punkt?

Dzisiaj zajmiemy się obszarami, głównie w oparciu o pakiet rgeos. Wpis raczej techniczny.

Po załadowaniu poniższych pakietów:

przygotujemy sobie dwa testowe obszary, na których przykładzie zobaczymy jak działają wybrane funkcje z pakietu rgeos.

Szczerze powiedziawszy mechanika budowania polygonów jest koszmarna, ale korzystając z przykładów w sieci udało się je przygotować :) Zobaczmy jak wyglądają nasze obszary:

Pomarańczową obwódką zaznaczony jest obszar poly1, zakreskowany niebieski to poly2. Kształty są odpowiednio dobrane, tak aby obszary zachodziły na siebie

gIntersects

Funkcja gIntersects zwraca TRUE, jeśli podane jako parametry obszary mają co najmniej jeden punkt wspólny. Zobaczmy jak to zadziała:

Dodatkowe parametry pozwoliły na uzyskanie kształtu z częścią wspólną (zaznaczony na czerwono). Jak widać wyszło znakomicie!

gDisjoint

Kolejna funkcja – gDisjoint zwraca TRUE, jeśli obszary nie mają punktów wspólnych:

Oczywiście w naszym przypadku jest to fałsz.

gContains

gContains zwraca TRUE, jeśli żaden z elementów poly1 nie znajduje się poza poly2, a co najmniej jeden punkt poly2 mieści się w poly1.

Też fałsz.

gContainsProperly

gContainsProperly zwraca TRUE w tych samych warunkach co gContains z dodatkowym wymogiem, że poly2 nie przecina się z granicą poly1.

gCovers

gCovers zwraca TRUE, jeśli żaden punkt poly2 nie jest wewnątrz poly1. To trochę różni się od gContains, ponieważ nie wymaga punktu w poly1, co może stanowić problem, ponieważ granice nie są uważane za leżące “wewnątrz” obszaru.

gCoveredBy

gCoveredBy jest przeciwieństwem gCovers i jest równoważne zamianie kolejności parametrów poly1 i poly2.

gWithin

gWithin jest przeciwieństwem gContains i jest równoważne zamianie poly1 i poly2.

Ta funkcja przyda nam się najbardziej.

Na początek potrzebujemy plików z mapami – są ogólnodostępne i za darmo na stronie Centralnego Ośrodka Dokumentacji Geodezyjnej i Kartograficznej. Pobieramy paczkę PRG – jednostki administracyjne.

Zobaczmy co zawiera mapa z województwami. Wczytujemy ją do R i zmieniamy układ współrzędnych na taki znany z geografii (są zapisane w innym układzie):

Po tym przekształceniu możemy je narysować. ggplot2 poradzi sobie z danymi typu SpatialPolygons:

ale możemy z nich wyciągnąć interesujące nas informacje – współrzędne opisujące granice obszarów i nazwy obszarów (lub ich kody TERYT – przykład łączenia danych z mapą pokazywałem już dawniej – kod TERYT jest świetnym kluczem łączącym) i upakować je w ramkę danych.

Z taką ramką ggplot też sobie poradzi:

Wyznaczmy teraz środki obszarów (województw) – z pomocą przychodzi funkcja coordinates z pakietu sp (rgdal go ładuje, nie trzeba tego robić bezpośrednio).

Przy okazji wyciągnąłem inną informację zawartą w plikach shp – powierzchnię województwa (to jakaś dziwna jednostka; mazowieckie ma 35559.20km2, tutaj wartość ta jest przemnożona przez 100). Bo pliki shp to nie tylko obszary, ale też dane. Zobaczmy co uzyskaliśmy:

Województwo Długość geograficzna Szerokość geograficzna Powierzchnia
opolskie 17.89988 50.64711 941272
świętokrzyskie 20.76909 50.76339 1171136
kujawsko-pomorskie 18.48822 53.07270 1797058
mazowieckie 21.09645 52.34576 3555920
pomorskie 17.98619 54.15424 1831001
śląskie 18.99410 50.33108 1233406
warmińsko-mazurskie 20.82493 53.85721 2417419
zachodniopomorskie 15.54329 53.58476 2289315
dolnośląskie 16.41069 51.08950 1994777
wielkopolskie 17.24310 52.33078 2982774
łódzkie 19.41760 51.60487 1821720
podlaskie 22.92931 53.26452 2018598
małopolskie 20.26933 49.85895 1518007
lubuskie 15.34275 52.19617 1398751
podkarpackie 22.16912 49.95367 1784523
lubelskie 22.90027 51.22072 2512291

Przygotujmy więc mapę województw z zaznaczonymi ich środkami, nazwami i informacją o powierzchni (wielkość i jasność kropki):

Przechodząc do postawionego na początku pytania (w jakim województwie leży dany punkt) weźmy środek lubelskiego:

i tylko jedno województwo:

Czy środek lubelskiego wypada w mazowieckim?

Nie leży. W zasadzie nie powinien ;-)

A teraz zrobimy to dla całej listy środków województw:

Województwo (środek) Czy w mazowieckim?
opolskie FALSE
świętokrzyskie FALSE
kujawsko-pomorskie FALSE
mazowieckie TRUE
pomorskie FALSE
śląskie FALSE
warmińsko-mazurskie FALSE
zachodniopomorskie FALSE
dolnośląskie FALSE
wielkopolskie FALSE
łódzkie FALSE
podlaskie FALSE
małopolskie FALSE
lubuskie FALSE
podkarpackie FALSE
lubelskie FALSE

Tylko jeden punkt powinien leżeć w mazowieckim i tak właśnie jest. Oczywiście jest to punkt środkowy tego województwa.

Możemy też sprawdzić hurtowo wszystkie punkty:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
0 T f f f f f f f f f f f f f f f
1 f T f f f f f f f f f f f f f f
2 f f T f f f f f f f f f f f f f
3 f f f T f f f f f f f f f f f f
4 f f f f T f f f f f f f f f f f
5 f f f f f T f f f f f f f f f f
6 f f f f f f T f f f f f f f f f
7 f f f f f f f T f f f f f f f f
8 f f f f f f f f T f f f f f f f
9 f f f f f f f f f T f f f f f f
10 f f f f f f f f f f T f f f f f
11 f f f f f f f f f f f T f f f f
12 f f f f f f f f f f f f T f f f
13 f f f f f f f f f f f f f T f f
14 f f f f f f f f f f f f f f T f
15 f f f f f f f f f f f f f f f T

Oczywiście po przekątnej mamy TRUE, bo środki województw ułożone są w tej samej kolejności co województwa.

Tyle o punktach, przejdźmy do obszarów. Na przykład poszukując województwa, w którym leży powiat. Potrzebujemy mapy (precyzyjniej: regionów) z powiatami. Jest w paczce ściągniętej z CODiG. Wczytujemy dane, dostosowujemy układ współrzędnych, tworzymy tabelkę z nazwami i dodajemy do niej informację czy powiat leży w województwie mazowieckim:

Ile powiatów mamy? Google mówi:

Województwo mazowieckie składa się z 37 powiatów i 5 miast na prawach powiatu. Powiaty dzielą się na 314 gmin – 35 miejskich, 50 miejsko-wiejskich i 229 wiejskich

Powiatów łącznie zatem jest 37+5=42, tak?

czy_mazowsze n
FALSE 338
TRUE 42

Tak! Które to powiaty?

lp powiat
1 powiat makowski
2 powiat Radom
3 powiat łosicki
4 powiat żyrardowski
5 powiat nowodworski
6 powiat Ostrołęka
7 powiat Warszawa
8 powiat ostrołęcki
9 powiat ostrowski
10 powiat otwocki
11 powiat piaseczyński
12 powiat płocki
13 powiat płoński
14 powiat pruszkowski
15 powiat przasnyski
16 powiat przysuski
17 powiat miński
18 powiat mławski
19 powiat pułtuski
20 powiat Siedlce
21 powiat radomski
22 powiat siedlecki
23 powiat sierpecki
24 powiat sokołowski
25 powiat sochaczewski
26 powiat szydłowiecki
27 powiat warszawski zachodni
28 powiat węgrowski
29 powiat Płock
30 powiat białobrzeski
31 powiat ciechanowski
32 powiat garwoliński
33 powiat gostyniński
34 powiat grodziski
35 powiat grójecki
36 powiat kozienicki
37 powiat zwoleński
38 powiat legionowski
39 powiat lipski
40 powiat wołomiński
41 powiat wyszkowski
42 powiat żuromiński

Do czego możemy to wykorzystać?

Kiedyś znalazłem jakąś listę polskich miejscowości z ich współrzędnymi geograficznymi. Nazwa, województwo, szerokość i długość geograficzna. Nie jestem pewien przypisania miejscowości do województwa. Na obrazku wygląda to tak:

Klikając w obrazek możesz go powiększyć

Specjalnie rozdzieliłem mapę na województwa, aby było widać że przypisania są błędne. Mamy 44555 miejscowości, całkiem sporo – być może są to wszystkie. Nie pamiętam źródła tych danych, ale jak widać nie są idealne. Co się dzieje w kujawsko-pomorskim (zaanektowało łódzkie), małopolskim (połowa lubuskiego została do niego przypisana), mazowieckim (nieco się rozrosło), mamy województwo “Polska” i tak dalej. No słabo.

Spróbujmy przypisać punkty do poprawnych województw opisanych mapą:

Mając takie dane możemy sprawdzić czy udało się poprawnie przypisać województwa:

Klikając w obrazek możesz go powiększyć

Nie wszystko udało się przypisać (niektóre punkty są poza mapą Polski – widać to już było na pierwszej mapce), ale teraz wygląda to już sensowniej. Nie było też województwa łódzkiego.

Zobaczmy macierz różnic (w formie graficznej):

Gdyby wszystko było od razu poprawnie przypisane mielibyśmy znaczniki tylko na przekątnej i tylko żółte. Już z mapy widzieliśmy, że tak nie było.

Mam też drugą listę – z kodami pocztowymi (ulica, miejscowość, gmina, kod pocztowy). Tylko co z tym zrobić? Zwykły join dwóch tabel nie pomoże. Pierwszy z brzegu (pierwsza z alfabetu) Adamów występuje na liście miejscowości 27 razy (rekord do Nowa Wieś – 116 razy) i bądź tu mądry który jest gdzie. A dodatkowo z powyższych mapek widać, że przypisania do województwa są niepoprawne.

Przy małych miejscowościach nie ma nazwy ulicy, więc przypisanie kodu pocztowego musi nastąpić na podstawie na przykład gminy (bo nazwy miejscowości nie są unikalne).

Bierzemy zatem listę miejscowości, mapę gmin (ta sama paczka z mapami), przypisujemy (analogicznie do przypisania miejscowości do województw wyżej) miejscowości do gmin i już mamy bardziej kompletne dane. Na koniec wystarczy złączyć odpowiednio dane z tabel miejscowość-gmina oraz gmina-kod pocztowy i otrzymamy pełny wykaz kodów pocztowych poszczególnych miejsc. Powinno pasować, nie sprawdzałem.

2 komentarze do wpisu „Mapy, mapy raz jeszcze

  1. Pytanie początkującego:
    Jakie znaczenie ma symbol %>% który używasz w kodzie?
    Nie ma chyba takiej definicji w R.
    To jest definicja z jakiegoś pakietu?

    • To tzw pipe. Pochodzi z pakietu tidyverse, w uproszczeniu przekazuje wynik funkcji która jest przed nim do kolejnej, jako jej pierwszy parametr.

Dodaj komentarz