Jak szybko przygotować mapkę z (mniejszymi) miastami Polski?
Szybki temat, bo i potrzeba była szybka. Żona zapytała jak ma zrobić mapkę z lokalizacją oddziałów jej firmy w poszczególnych miastach Polski. I nie były to miasta wojewódzkie (to znaczy takie też), ale też powiatowe.
Szukanie mapy konturowej naszego kraju, a później szukanie gdzie na owej mapie zaznaczyć jakiś przysłowiowy Pcim Dolny to nie jest zadanie na godzinę dwudziestą trzecią. Zajęłoby (mi) to jakieś dwie albo trzy godziny. Od czego jest R? Rozwiązania na szybko przyszły mi do głowy dwa, oba na podstawie doświadczeń oczywiście. Lubię mapy i dane na mapach, to mam też trochę doświadczenia w tym temacie.
Wersja z Google Maps
Skorzystamy z biblioteki ggmap do znalezienia punktów na mapie, pobrania mapy oraz jej narysowania – ta biblioteka to po prostu interface R do GoogleMaps (ale nie tylko). Dodatkowe elementy narysujemy z pomocą ggplot2.
|
1 2 |
library(ggplot2) library(ggmap) |
Przygotujmy sobie dane. Wybrałem z Wikipedii trochę (losowo) miast powiatowych. Najpierw wtłoczymy je w ramkę danych i przygotujemy w niej odpowiednie pola na informacje dodatkowe:
|
1 2 3 4 5 6 7 8 9 10 11 |
miasta_df <- data.frame(nazwa = c("Legnica", "Grudziądz", "Włocławek", "Chełm", "Gorzów Wielkopolski", "Piotrków Trybunalski", "Tarnów", "Siedlce", "Krosno", "Tarnobrzeg", "Łomża", "Słupsk", "Bielsko-Biała", "Piekary Śląskie", "Żory", "Elbląg", "Kalisz", "Leszno", "Świnoujście"), stringsAsFactors = FALSE) # miejsce na długość i szerokość geograficzną miasta_df$long <- NA miasta_df$lat <- NA |
Teraz dla każdego z miast pobierzemy współrzędne geograficzne. Zapytamy o nie po prostu Google. Nazwę miasta konwertuję (iconv()) na UTF8 (system na którym używam R to Windows, więc znaki kodowane są w Latin2, z którym internety mają raczej problem) i dodaję “Polska”, coby doprecyzować moje zapytanie. W odpowiedzi dostajemy zmienną ze współrzędnymi – dopisujemy je do stosownego wiersza w całej tabeli.
|
1 2 3 4 5 |
for(i in 1:nrow(miasta_df)) { loc <- geocode(paste0(iconv(miasta_df[i, "nazwa"], to="UTF8"), ", Polska")) miasta_df[i, "long"] <- loc$lon miasta_df[i, "lat"] <- loc$lat } |
Mamy wszystko co trzeba:
| nazwa | long | lat |
|---|---|---|
| Legnica | 16.15532 | 51.20701 |
| Grudziądz | 18.75356 | 53.48375 |
| Włocławek | 19.06774 | 52.64833 |
| Chełm | 23.47120 | 51.14312 |
| Gorzów Wielkopolski | 15.23693 | 52.73253 |
| Piotrków Trybunalski | 19.70302 | 51.40517 |
| Tarnów | 20.98584 | 50.01210 |
| Siedlce | 22.29016 | 52.16760 |
| Krosno | 21.76605 | 49.68248 |
| Tarnobrzeg | 21.67907 | 50.57291 |
| Łomża | 22.05903 | 53.17812 |
| Słupsk | 17.02848 | 54.46415 |
| Bielsko-Biała | 19.05838 | 49.82238 |
| Piekary Śląskie | 18.92705 | 50.37891 |
| Żory | 18.70064 | 50.04472 |
| Elbląg | 19.40449 | 54.15606 |
| Kalisz | 18.08535 | 51.76728 |
| Leszno | 16.59375 | 51.84199 |
| Świnoujście | 14.24758 | 53.91003 |
więc pora to narysować. Pobierzmy najpierw tło – google’owę mapę naszego kraju. Wiecie, że środek Polski leży w okolicach Kutna? Dlatego poniżej pierwszym parametrem (środek obszaru) jest właśnie Kutno. Taka wiedza ze szkoły czasem się przydaje ;) Eksperymentalnie dobieramy (metodą prób i błędów) wielkość parametru zoom i w efekcie otrzymujemy coś takiego:
|
1 2 |
polska <- get_map("Kutno, Polska", zoom=6) ggmap(polska) |
Jest OK. Jak trochę przyciemnimy, dodamy kontrastowe punkty z miastami i ich nazwami będzie ok. Proszę bardzo:
|
1 2 3 4 5 6 7 8 |
ggmap(polska, darken = 0.6) + geom_point(data=miasta_df, aes(long, lat), size=3, color="red") + geom_text(data=miasta_df, aes(long, lat, label=nazwa), size=3, color="white", hjust=-0.1) + theme_void() |
Wersja bez Google Maps
Może jednak z różnych względów nam się to nie podobać (bo na mapie z GMaps są napisy, bo za dużo kolorów, bo cokolwiek). Wykorzystamy więc dane wektorowe (pliki Shapefile, więcej mięsa po angielsku znajdziecie).
Najpierw jednak trzeba mieć te pliki. Są ogólnodostępne i za darmo na stronie Centralnego Ośrodka Dokumentacji Geodezyjnej i Kartograficznej i przydadzą się wielokrotnie. Interesuje nas paczka “PRG – jednostki administracyjne”, która zawiera dane o obrębach gmin, powiatów oraz województw. Archiwum trzeba ściągnąć i rozpakować. Najpierw jednak potrzebujemy w samym R biblioteki rgdal ze wszystkimi dodatkami. Zainstalujemy ją z CRANa.
|
1 2 |
install.packages("rgdal") library(rgdal) |
W kodzie musimy odpowiednio podać ścieżkę do katalogu z plikami dla województw:
|
1 2 3 4 |
wojewodztwa <- readOGR("ŚCIEŻKA_DO_MAP", "wojewodztwa") wojewodztwa <- spTransform(wojewodztwa, CRS("+init=epsg:4326")) wojewodztwa <- fortify(wojewodztwa, region = "jpt_kod_je") |
W wyniku tego działania dostajemy ramkę (dłuuugą) z punktami wyznaczającymi granice województw. Możemy ją narysować i od razu nałożyć na nią nasze punkty z miastami:
|
1 2 3 4 5 6 7 8 9 10 11 12 |
ggplot() + geom_polygon(data=wojewodztwa, aes(long, lat, group=group), fill="white", color="gray") + geom_point(data=miasta_df, aes(long, lat), size=3, color="red") + geom_text(data=miasta_df, aes(long, lat, label=nazwa), size=3, color="red", hjust=-0.2) + coord_map() + theme_void() |
Gotowe! Jakieś 15-20 minut roboty, tylko kilkanaście linii kodu.
Gotowy obrazek eksportujemy do pliku (można przez ggsave() od razu albo z poziomu RStudio) i przekazujemy żonie do prezentacji. A nie minęła jeszcze północ…
Potrzeba matką wynalazku ;)
Chyba tracę chęci:
> wojewodztwa <- readOGR("D:/R/PRG", "wojwodztwa")
Error in ogrInfo(dsn = dsn, layer = layer, encoding = encoding, use_iconv = use_iconv, : Cannot open layer
Googlam i chyba rozbija się tu o wersję R lub biblioteki. Nie znalazłem jednak odpowiedzi, która by mnie utwierdziła w tym przekonaniu. Czuję się jak 17 lat temu gdy bawiłem się lnuxami i była zabawa by tak zgrać wersje by zależności między wszystkimi potrzebnymi pakietami zatrybiły. Chyba za stary jestem na takie zabawy.
Przeczytaj dokładnie nazwę warstwy (drugi argument readOGR w Twoim kodzie, który wklejasz)… literka po literce.
A jeśli to nie to, to update pakietów może pomoże?
Łukasz, nie wiem jakie masz doświadczenie z plikami shp ale może natknąłeś się na taki problem. Gdy zamiast:
wojewodztwa <- fortify(wojewodztwa, region = "jpt_kod_je")
dam:
wojewodztwa <- fortify(wojewodztwa, region = "jpt_nazwa_")
to mapa źle się rysuje, jakby dziecko pokreśliło je kredką. Jest to dziwne bo QGIS zauważyłem, że w tym polu jest nazwa województwa a ona powinna odpowiadać jednoznacznie polu jpt_kod_je. Zależało mi by w dataframe mieć nazwę województwa gdyż z bazy danych pobieram sobie jakieś tam dane dla województw i dzięki temu łatwo i wygodnie mogę matchować te dane ze sobą.
Dodam, że robię to z kolorami bo koloruję mapę z użyciem heatmapy
Kod bez samego zapytania do bazy wklejam poniżej. Z tego do czego sam doszedłem, to, że wszystko rozwala się po tym jak merguję wartości kolorów z data frame z danymi mapy:
districts <- readOGR("D:/R/PRG", "wojewodztwa")
districts <- spTransform(districts, CRS("+init=epsg:4326"))
districts <- fortify(districts, region = "jpt_nazwa_")
distcol <- data.frame(
id = sqlquery["WOJEWODZTWO"],
value = heat.colors(16)
)
distplot <- merge (distcol,districts, by.x = "WOJEWODZTWO", by.y = "id")
ggplot() +
geom_polygon(data=distplot,
aes(long, lat, group=group, fill = I(value)),
color="#000000", show.legend = FALSE) +
coord_map() +
theme_void()
Sprawdź kolejność punktów w data frame przed i po połączeniu, jest jedna kolumna która odpowiada za to, widać po danych jak wartości w niej rosną. Być może merge coś z nią robi – dlatego wolę używać left_join.
Ggplot rysuje linie i łączy je po kolei w polygony. Jeśli linie (punkty końcowe linii) nie będą po kolei to wyjdą takie maziaje.
A dane łączyć lepiej jakimś uniwersalnym kluczem niż nazwą – ja polecam kod TERYT. Mogą się trafić dwa rejony o tej samej nazwie (chyba są dwa powiaty, a na pewno tak jest z gminami).
Do robienia map i analiz przestrzennych polecam rozwiązania GISowe np. QGIS -darmowy opensource i bardzo przyjemny w nauce a dający wiele możliwości wizualizacyjnych i analitycznych (przestrzennie) :)
Nigdy nic nie robiłem w R, ale pracuję z bazami danych i wizualizacjami. Moje pytanie: jak w R narysować na mapie Polski strukturę składającą się np z kilku powiatów i do całości przypiąć wykres? Załóżmy, że posługujemy się shapefile i GUS TERYT.