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:

library(rgdal)
library(rgeos)
library(tidyverse)
library(broom) # for tidy=fortify

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

poly1 <- SpatialPolygons(list(Polygons(list(Polygon(coords=matrix(c(0, 0, 2, 2, 0, 1, 1, 0),
                                                                  ncol=2, byrow=FALSE))), ID=c("a")),
                              Polygons(list(Polygon(coords=matrix(c(0, 0, 2, 2, 2, 3, 3, 2),
                                                                  ncol=2, byrow=FALSE))), ID=c("b"))))

poly2 <- SpatialPolygons(list(Polygons(list(Polygon(coords=matrix(c(0, 0, 2, 2, 1, 1, 0.5, 3, 3, 0, 0, 2),
                                                                  ncol=2, byrow=FALSE))), ID=c("c"))))

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:

plot(poly1, border="orange")
plot(poly2, border="blue", add=TRUE, lty=2, density=8, angle=30, col="blue")

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:

gI <- gIntersection(poly1, poly2, byid=TRUE, drop_lower_td=TRUE)
plot(gI, add=TRUE, border="red", lwd=3)

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:

gDisjoint(poly1, poly2)

## [1] FALSE

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.

gContains(poly1, poly2)

## [1] FALSE

Też fałsz.

gContainsProperly

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

gContainsProperly(poly1, poly2)

## [1] FALSE

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.

gCovers(poly1, poly2)

## [1] FALSE

gCoveredBy

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

gCoveredBy(poly1, poly2)

## [1] FALSE

gWithin

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

gWithin(poly1, poly2)

## [1] FALSE

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

wojewodztwa <- readOGR("../!mapy_shp/wojewodztwa.shp", "wojewodztwa")

wojewodztwa <- spTransform(wojewodztwa, CRS("+init=epsg:4326"))

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

ggplot(wojewodztwa) +
   geom_polygon(aes(long, lat, group=group, fill=group), color="gray", show.legend = FALSE) +
   coord_map() +
   theme_void()

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.

wojewodztwa_nazwy <- wojewodztwa@data %>% select(jpt_kod_je, jpt_nazwa_)
wojewodztwa_df <- tidy(wojewodztwa, region = "jpt_kod_je")
wojewodztwa_df <- left_join(wojewodztwa_df, wojewodztwa_nazwy, by=c("id"="jpt_kod_je"))

Z taką ramką ggplot też sobie poradzi:

ggplot(wojewodztwa_df) +
   geom_polygon(aes(long, lat, group=group, fill=jpt_nazwa_), color="gray") +
   coord_map() +
   theme_void()

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

srodki <- coordinates(wojewodztwa) %>%
   as_tibble() %>%
   set_names(c("long", "lat")) %>%
   mutate(nazwa = wojewodztwa@data$jpt_nazwa_,
          powierzchnia = wojewodztwa@data$jpt_powier) %>%
   select(nazwa, long, lat, powierzchnia)

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:

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

ggplot() +
   geom_polygon(data = wojewodztwa_df, aes(long, lat, group=group),
                fill = "gray80", color="gray30", show.legend = FALSE) +
   geom_point(data = srodki, aes(long, lat,
                                 size = powierzchnia, color = powierzchnia),
              show.legend = FALSE) +
   geom_text(data = srodki, aes(long, lat, label = nazwa), vjust = 1.7, size = 2.8) +
   coord_map() +
   theme_void() + 
   scale_color_gradient(low = "red", high = "yellow") + 
   scale_size(range=c(2,6))

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

punkt <- SpatialPoints(matrix(c(22.90027, 51.22072), ncol=2), CRS("+init=epsg:4326"))

i tylko jedno województwo:

mazowsze <- wojewodztwa[wojewodztwa@data$jpt_nazwa_ == "mazowieckie", ]

Czy środek lubelskiego wypada w mazowieckim?

gWithin(punkt, mazowsze)

## [1] FALSE

Nie leży. W zasadzie nie powinien ;-) A teraz zrobimy to dla całej listy środków województw:

# bierzemy wszystkie punkty
punkty <- SpatialPoints(srodki[,c(2:3)], CRS("+init=epsg:4326"))

# czy leżą w mazowieckim?
srodki$czy_mazowsze <- gWithin(punkty, mazowsze, byid = TRUE) %>% as.logical() 

# wynik:
srodki %>% select(nazwa, czy_mazowsze)

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:

gWithin(punkty, wojewodztwa, byid = TRUE)

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:

powiaty <- readOGR("../!mapy_shp/powiaty.shp", "powiaty")

powiaty <- spTransform(powiaty, CRS("+init=epsg:4326"))

powiaty_df <- tibble(powiat = powiaty@data$jpt_nazwa_)

powiaty_df$czy_mazowsze <- gWithin(powiaty, mazowsze, byid = TRUE) %>% as.logical()

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?

powiaty_df %>% count(czy_mazowsze)

Tak! Które to powiaty?

powiaty_df %>% filter(czy_mazowsze) %>% mutate(lp = row_number()) %>% select(lp, powiat)

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:

library(readxl)
miejscowosci <- read_excel("Wspolrzedne_miejscowosci.xlsx") %>%
   set_names(c("Miejscowosc", "Wojewodztwo", "Lat", "Long"))

ggplot() +
   # najpierw punkty, a na nie granice województw
   geom_point(data = miejscowosci, aes(Long, Lat, color = Wojewodztwo), alpha = 0.5, show.legend = FALSE, size = 0.9) +
   geom_polygon(data = wojewodztwa_df, aes(long, lat, group=group), fill = NA, color = "gray30") +
   facet_wrap(~Wojewodztwo) +
   coord_map() +
   theme_void()

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

# zwróć uwagę, że kolumny long i lat są odwrotnie!
miejscowosci_sp <- SpatialPoints(miejscowosci[,c(4,3)], CRS("+init=epsg:4326")) 

# macierz przypisania punktu do województwa - to trochę trwa
miejscowosci_mat <- gWithin(miejscowosci_sp, wojewodztwa, byid = TRUE) %>% matrix(ncol = 16, byrow = TRUE)

# dla każdego wiersza macierzy szukamy kolumny (województwa) z TRUE i wpisujemy do tabeli miejscowosci nazwę województwa
miejscowosci$wojewodztwo_mapa <- miejscowosci_mat %>%
   apply(1, function(x) min(which(x, arr.ind = TRUE), na.rm=TRUE)) %>%
   as.character(wojewodztwa@data$jpt_nazwa_)[.]

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

ggplot() +
   geom_point(data = miejscowosci, aes(Long, Lat, color = wojewodztwo_mapa), alpha = 0.5, show.legend = FALSE, size = 0.9) +
   geom_polygon(data = wojewodztwa_df, aes(long, lat, group=group), fill = NA, color = "gray30") +
   scale_color_manual(values = rainbow(17)) +
   facet_wrap(~wojewodztwo_mapa) +
   coord_map() +
   theme_void()

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

miejscowosci %>%
   count(Wojewodztwo, wojewodztwo_mapa) %>%
   ungroup() %>%
   ggplot() +
   geom_tile(aes(Wojewodztwo, wojewodztwo_mapa, fill = n), color = "gray80", show.legend = FALSE) +
   scale_fill_distiller(palette = "YlGnBu") +
   labs(x = "Województwo w danych", y = "Województwo odczytane ze współrzędnych") + 
   theme(axis.text.x = element_text(angle=90, hjust=1, vjust=0))

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.