Po przeszukaniu połowy internetu znalazłem tylko dwa posty w języku polskim na temat tworzenia map z wykorzystaniem języka R(z czego wynika że jeszcze gdzieś są jakieś dwa inne).
Nie jest to popularny temat, dlatego też zdecydowałem się parę słów napisać o tym tak żeby ktoś nie musiał przechodzić przez tą samą drogę co ja.
Cała moja wiedza w tym temacie pochodzi z kilku bardzo dobrych stron do których linki dam na końcu. Tutaj zaprezentuje streszczenie tego i pokażę działający przykład.
Najpierw należy się zaopatrzyć w mapy, można je pobrać bezpłatnie z tego linka http://www.gugik.gov.pl/pzgik/
Jest tam do wyboru kilka formatów, nas będzie interesowało rozszerzenie shp i zestaw map: "PRG – jednostki administracyjne". Jest to paczka z podziałami na Gminy, Powiaty, Województwa, Polska.
Zaimportujmy je teraz do RStudio. Służy do tego polecenie readOGR znajdziemy je w bibliotece rgdal. Są oczywiście inne metody ale ta działa*.
*Gdy chciałem wczytać po raz pierwszy te dane to natknąłem się na znany problem ale nie rozwiązany do tej pory: "dlaczego u jednego działa a u drugiego nie?", nie będę wnikał w szczegóły ale chodziło o wersje R.
Niektóre artykuły z przed 5 lat i więcej mogą zawierać nieaktualne składnie poleceń, trzeba uważać na tego typu rzeczy bo czasem się zdarzają.
Składnia polecenia readOGR w szczegółach nie jest aż tak istotna, najważniejsze są dla nas dwie zmienne:
dns - ścieżka do pliku shp,
layer - nazwa warstwy, najczęściej jest to nazwa pliku bez rozszerzenia shp, dla nas warstwami będą Gminy, Powiaty, Województwa, Polska bo taki mamy podział terytorialny, w mapach dla innych krajów będzie inaczej.
Jeśli chcemy zaimportować np. województwa to uruchamiamy kod:
wojewodztwa<-readOGR(dsn = "ścieżka do pliku : Województwa.shp",layer= "Województwa")
Tak utworzona zmienna jest wektorem przestrzennym(można tu mówić o danych geoprzestrzennych). Możemy sobie ją podglądnąć i zobaczymy, że ma ona 5 slotów, z czego najważniejsze są dla nas : data, polygons, proj4string.
W data mamy wszystkie informacje o danej warstwie włącznie z nazwami(JPT_NAZWA) i kodami terytorialnymi(JPT_KOD_JE lub JPT_KJ_I_1). Niestety nazwy nie mają polskich znaków, chyba że ja coś robię nie tak.
Natomiast kody terytorialne są unikatowe dla każdej jednostki administracyjnej, w przypadku województw są to numery różniące się o 2 : '02','04'... dostosowane do kolejności alfabetycznej województw rosnąco.
Możemy podglądnąć sobie wygląd tej mapy stosując polecenie:
plot(wojewodztwa)
Funkcja plot daje nam kilka możliwości modyfikacji takiej mapy. Pokolorujmy na niebiesko wszystkie województwa których kody(JPT_KJ_I_1) są większe niż '20'.
Dobrze jest sprawdzić na początek jaki typ zmiennej ma JPT_KJ_I_1. W tym celu wykonajmy polecenie:
sapply(X=wojewodztwa@data,FUN=
gdzie x to nasz wektor a parametr FUN określa co chcemy wyświetlić, w tym przypadku "strukturę wektora". Zwracam też uwagę, że do slotów odwołujemy się poprzez "@" a nie "$".
Możemy zobaczyć że JPT_KJ_I_1 ma typ "factor" co nas nie cieszy, bo jeśli będziemy chcieli wybrać te województwa to musimy użyć ">" co nie ma sensu dla factor. Dlatego też przekonwertujemy go na typ numeric.
wojewodztwa$JPT_KJ_I_1 <-as.numeric(wojewodztwa$JPT_
Składnia jest bardzo prosta. W takich operacjach należy uważać na potencjalne obcięte 0 na początkach, natomiast w naszym przypadku nie należy się tym przejmować.
W jaki sposób teraz pokolorować te województwa?
Plan jest prosty. Wybierzemy tylko te województwa, które nas interesują i pokolorujemy je na niebiesko:
woj<-wojewodztwa$JPT_KJ_I_1>20
powstaje zmienna logiczna true/false, true przy tych województwach które chcemy pokolorować
teraz użycie polecenia
plot(wojewodztwa[woj,])
narysuje nam tylko te wybrane województwa a pokolorowanie ich to polecenie
plot(wojewodztwa[woj,], col="blue")
Chcielibyśmy jednak mieć całą mapę i te wybrane województwa w tym celu nałożymy na siebie dwie mapy lub bardziej poprawnie dwie warstwy:
plot(data=wojewodztwa) # pierwsza warstwa
plot(data=wojewodztwa[woj,], col="blue", add=TRUE) # druga warstwa
opcja "add" pozwala na narysowanie dwóch wykresów na jednym równocześnie. Uff...
Coś jeszcze z tego plot można by wycisnąć ale bardziej interesujące polecenie to ggplot2. Daje ono dużo większe możliwości no i same wykresy czy też mapy wyglądają zdecydowanie lepiej.
Znajdziemy je w pakiecie tidyverse, który polecam zainstalować w całości.
Jest jeden problem z ggplot2, przyjmuje on typ danych: data frame a nasze dane to dane geoprzestrzenne, należy zatem przekonwertować je do odpowiedniego formatu.
Użyjemy do tego funkcji tidy() z pakietu broom(opcjonalnie można też skorzystać z funkcji fortify() jakaś jest różnica ale dla nas bez znaczenia).
wojewodztwa_df <- tidy(wojewodztwa)
Na wyniku otrzymujemy zwykły wektor z siedmioma kolumnami: long, lat, order, hole, piece, gruop, id.
A mówiąc ściślej 440tyś obserwacji dla 7 zmiennych, skąd tego tyle skoro mamy 16 województw? W dalszej części znajdziemy odpowiedź.
Narysujmy w takim razie mapkę przy użyciu ggplot. Teraz się zaczyna...
Może o samym ggplot dwa zdania. Pod tym linkiem jest świetne wprowadzenie w tą tematykę. Jest to pierwszy rozdział z książki "Jezyk R Kompletny zestaw narzędzi dla analityków danych"
http://www.super-booki.pl/pdf/
Szablon dla tej funkcji jest następujący:
ggplot(data =<DANE>) + <funkcja geometryczna>(mapping=aes(<
Mamy tutaj również do czynienia z dwoma warstwami, pierwsza z nich rysuje pusty wykres a druga ( zamiast add=true mamy znak "+") umieszcza na nim w zależności od funkcji geometrycznej
punkty/figury/wielokąty, wykresy skrzynkowe i wiele innych... Mapping definiuje sposób mapowania zmiennych ze zbioru danych(data) na właściwości wizualne. Argument mapping zawsze
występuje w parze z funkcją aes(), natomiast sama funkcja aes() przyjmuje podstawowe argumenty x i y określające, które zmienne
należy zmapować na osie x i y. W naszym przypadku lat, long czyli długość i szerokość.
Funkcją geometryczną, którą wykorzystamy jest geom_polygon(), rysuje ona wielokąty. Z tąd można wnioskować, że każde województwo(czy inna warstwa mapy) jest wielokątem.
Każda z 440tyś obserwacji reprezentuje jeden punkt wielokąta i ma ona współrzędne (lat,long) i teraz bardzo prosto to działa: funkcja aes() otrzymuje współrzędne i je rysuje.
W takim razie uruchamiamy polecenie:
ggplot(data=wojewodztwa_df) + geom_polygon(aes(x=lat,y=long)
coś poszło nie tak... Kilka rzeczy. Po pierwsze dziwna sprawa ale mapy od gis-u mają zamienione współrzędne... dlatego należy zamienić (x=long, y =lat)
ggplot(data=wojewodztwa_df) + geom_polygon(aes(x=long,y=lat)
wygląda lepiej ale dalej coś nie gra. Musimy przekazać funkcji informację jak grupować obserwacje w wielokąty. jest to domyślna wartość: group=group.
Ten zapis powinien się znaleźć w aes() jako, że tyczy się wizualizacji zmiennych.
ggplot(data=wojewodztwa_df) + geom_polygon(aes(x=long,y=lat, group=group))
Tym razem wygląda to już całkiem dobrze. z tym że mamy tylko mapę Polski bez podziału na województwa. Musimy dodać kontury województw czyli ich kolor(color) i grubość lini(size).
ggplot(data=wojewodztwa_df) + geom_polygon(aes(x=long,y=lat, group=group),color="black",
Domyślne wypełnienie wielokątów jest ustawione na kolor "grey20" dlatego też zmieńmy je na jaśniejszy kolor.
ggplot(data=wojewodztwa_df) + geom_polygon(aes(x=long,y=lat, group=group),color="black",
Jeszcze gdyby ktoś chciał można ustawić przeźroczystość(alpha), jakieś modyfikacje z liniami i tym podobne ale już się w to nie bawmy.
Generalnie wygląda to już przyzwoicie. A teraz największa zagadka, czyli dlaczego ta mapka jest tak rozciągnięta? Jesteśmy przyzwyczajeni do nieco innego wyglądu.
Chodzi tutaj o Geodezyjny System Odniesienia. Różne mapy mogą mieć różny. Te pobrane przez nas mają GRS80 lub inaczej WSG84 lub inaczej EPSG:4326.
Można to sprawdzić na dwa sposoby. Otwierając plik wojewodztwa.prj lub z użyciem funkcji
st_crs(wojewodztwo) która znajduje się w pakiecie sf.
WSG84 - Jest to jednolity dla całego świata układ współrzędnych. Dodatkowo jest powszechnie stosowany w urządzeniach do nawigacji. Charakteryzuje go m.in. równoleżnikowe rozciągnięcie danych. My natomiast potrzebujemy EPSG:2180. Jest to jednolity dla obszaru całej Polski układ stosowany dla map w skali poniżej 1:10 000.
Gdyby kogoś interesowało to tutaj jest dużo wiedzy na ten temat, całkiem nieźle opisanej: https://gis-support.pl/co-to-
Zmiana polega na dołożeniu informacji, o tym z jakiego układu współrzędnych będziemy korzystać.
EPSG2180 <- CRS('+init=epsg:2180') # zachowujemy informację
wojewodztwa <- spTransform(wojewodztwa, EPSG2180) # transformujemy nasze dane
Zmiany dokonujemy na danych geoprzestrzennych dlatego też aby móc narysować mapkę najpierw ponownie trzeba zmienić nasze zaktualizowane dane na data frame.
wojewodztwa_df <- tidy(wojewodztwa)
ggplot(data=wojewodztwa_df) + geom_polygon(aes(x=long,y=lat, group=group),color="black",
Wygląda prawie dobrze, jest jednak pewien problem z wyświetlaniem i zachowaniem skali na obu osiach, funkcja coord_equal()
jest tak zdefiniowana, aby dać równe "skale" na obu osiach.
ggplot(data=wojewodztwa_df) + geom_polygon(aes(x=long,y=lat, group=group),color="black",
Jest to bardzo podstawowa mapa. Domyślny punkt wyjścia do dalszych prac nad nią. Poniżej zebrałem cały kod który należy wykonać by uzyskać powyższą mapkę:
wojewodztwa<-readOGR(dsn = "ścieżka do pliku : Województwa.shp",layer= "Województwa")
EPSG2180 <- CRS('+init=epsg:2180') # zachowujemy informację
wojewodztwa <- spTransform(wojewodztwa, EPSG2180) # transformujemy nasze dane
wojewodztwa_df <- tidy(wojewodztwa)
ggplot(data=wojewodztwa_df) + geom_polygon(aes(x=long,y=lat, group=group),color="black",
Jak na razie bez większych problemów natomiast teraz się zaczyna trudniejsza część. Chcemy zrobić to samo co przy użyciu funkcji plot().
Zapewne można nałożyć warstwy jak poprzednio ale chciałbym to zrobić przy użyciu jednego polecenia i taką możliwość daje nam ggplot.
Najpierw musimy, jak poprzednio, dołożyć kolumnę "grupującą" województwa >20 i <=20.
Jak przyglądamy się wektorowi wojewodztwa_df i wektorowi przestrzennemu wojewodztwa widzimy, że różnią się one zdecydowanie.
Funkcja tidy() oprócz zmiany danych na data frame obcięła też zakres danych do minimum.
Zatem należy dodać interesującą nas kolumnę z wojewodztwa@data, w tym przypadku JPT_KOD_JE (kody terytorialne).
Problem polega na tym że nie mamy klucza połączenia między wojewodztwa_df i województwa nie mają one wspólnej kolumny.
Na szczęście możemy wykorzystać to w jaki sposób powstał wektor wojewodztwa_df.
Przyjrzyjmy się kolumnie wojewodztwa_df$id jest w niej przechowywana informacja o tym do której warstwy są przypisane obserwacje. Jeśli województw mamy 16 to id ma zakres 0-15.
Gdyby udało nam się przypisać dany nr id do kodu województwa to otrzymamy to co chcemy.
Podczas tworzenia wojewodztwa_df jest zachowywana kolejność czyli bierzemy pierwsze województwo z wojewodztwa@data$JPT_KOD_JE i nadajemy my id=0 dorzucamy obserwacja i bierzemy kolejny nr itd.
Zatem ponumerujmy, nie zmieniając układu wojewodztwa@data$JPT_KOD_JE kolejnym nr 0-15 i zamieńmy na typ character gdyż wojewodztwa_df$id ma typ charakter.
temp_df <- data.frame(wojewodztwa@data$
temp_df$id <- as.character(seq(0,nrow(temp_
połączmy się teraz temp_df z wojewodztwa_df
wojewodztwa_df <- left_join(wojewodztwa_df, temp_df, by="id") # nie tylko się łączymy ale dokładamy wszystko z temp_df poza kluczem połączenia
Na podglądzie widać, że mamy teraz 2 dodatkowe zmienne w tablicy wojewodztwa_df. Pozostaje nam podzielenie województw na grupy i zapisanie informacji do dowolnej kolumny, powiedzmy id2.
wojewodztwa_df$id2 <-substr(wojewodztwa_df$
Aby narysować mapkę w takim grupowaniu należy dołożyć do podstawowego polecenia ryzującego mapę informacje, że chcemy ją pokolorować/wypełnić kolorem wielokąty na które wskazuje podział id2.
aes() należy dodać fill="id2" wtedy jest to informacja dla funkcji geom_polygon, że dla tych samych grup ma dać te same kolory( na razie domyślne).
ggplot(data=wojewodztwa_df) + geom_polygon(aes(x=long,y=lat,
Zwracam uwagę, że fragment: fill="lightgrey" jest już niepotrzebny skoro kolorujemy domyślnie według grupowania id2.
