Geolokalizacja w R na przykładzie instalacji termicznego przekształcania osadów ściekowych (ITPOŚ)

Wpis z serii technicznych, spowodowany konieczność jaka zaistniała. Przygotowując się do prezentacji na konferencję „Innowacje w Inżynierii Środowiska” zaistniała bowiem potrzeba przygotowania slajdu z lokalizacją instalacji termicznego przekształcania osadów ściekowych (spalarni osadów ściekowych) na mapie Polski i wyróżnienie tych instalacji, które odpowiedziały na przesłaną uprzednio ankietę.

W Polsce funkcjonuje aktualnie 11 takich obiektów, które są zlokalizowane następująco:

• Tama Pomorzańska 8, 70-030 Szczecin
• Długa 28, 84-230 Dębogórze Wybudowanie
• Benzynowa 26, 83-011 Gdańsk
• Leśna 10-173, 11-041 Olsztyn
• gen. Tadeusza Bora-Komorowskiego 72, 85-950 Bydgoszcz
• Zjazd 23, 18-400 Łomża
• Czajki 4/6, 03-054 Warszawa
• Sanitariuszek 66, 91-867 Łódź
• Sportowa, 66-016 Łężyca
• Przemysłowa 53, 26-052 Sitkówka-Nowiny
• Kosiarzy 3, 33-332 Kraków

Myślę, że warto z zastosowanymi w nich technologiami również się zaponać, więc pewnie w jednym z kolejnych wpisów będę chciał to zrobić. Dziś jednak zademonstruję, jak w prosty sposób uzyskać współrzędne geograficzne, które zidentyfikują nam położenie instalacji termicznego przekształcenia osadów ściekowych (ITPOŚ).

Na początku utworzę strukturę danych (ramka danych), która zawierać będzie następujące pola: adres lokalizacji ITPOŚ, długość i szerokość geograficzną oraz informację o tym, czy ITPOŚ wypełnił ankietę. Zatem do pracy.

miejsca_df <- data.frame(adres=c("Tama Pomorzańska 8, 70-030 Szczecin","Długa 28, 84-230 Dębogórze Wybudowanie","Benzynowa 26, 83-011 Gdańsk","Leśna 10-173, 11-041 Olsztyn","gen Tadeusza Bora-Komorowskiego 72, 85-950 Bydgoszcz","Zjazd 23, 18-400 Łomża","Czajki 4/6, 03-054 Warszawa","Sanitariuszek 66, 91-867 Łódź","Sportowa, 66-016 Łężyca","Przemysłowa 53, 26-052 Sitkówka-Nowiny","Kosiarzy 3, 33-332 Kraków"), stringsAsFactors = FALSE)
miejsca_df$long <- NA
miejsca_df$lat <- NA
miejsca_df$ankieta <- c("tak","tak","tak","nie","tak","nie","tak","tak","tak","tak","tak")

Teraz dla każdego miejsca pobierzemy współrzędne geograficzne.

Jeszcze nie tak dawno można było wykorzystać jedną z wielu dostępnych aplikacji jak np. http://pl.mygeoposition.com/. Korzystały one z darmowej usługi Google Maps. W czerwcu ubiegłego roku Google wprowadziło zmiany. O tym, jak wygląda obecny model i co zrobić by z niego korzystać za darmo w załączonym artykule. Napomknę tylko, że nadal można korzystać z darmowego dostępu do usługi Google Maps, gdy wykonywane działania nie przekroczą 300 USD miesięcznie. Jak wyglądają ceny korzystania z usługi można znaleźć tutaj.

Ponieważ pracuję w systemie Windows dokonam konwersji kodowania znaków z Latin2 do UTF8 za pomocą funkcji iconv() i dodam na końcu szukanej lokalizacji tekst Polska.

for(i in 1:nrow(miejsca_df)) {
  loc <- geocode(paste0(iconv(miejsca_df[i, "adres"], to="UTF8"), ", Polska"))
  miejsca_df[i, "long"] <- loc$lon
  miejsca_df[i, "lat"] <- loc$lat
}

Jak widać w celu geolokalizacji skorzystałem z funkcji geocode() z biblioteki ggmap, stanowiącej łącznik pomiędzy programem R a Google Maps. W rezultacie otrzymujemy:

Teraz wystarczy to wyrysować na mapie. Skorzystam ze schematu przedstawionego już uprzednio we wpisie http://ekordo.pl/mimo-ze-wytwarzamy-coraz-wiecej-odpadow-komunalnych-to-nadal-jestesmy-w-koncowce-europy/. Pobieram dane wektorowe jednostek administracyjnych Polski w postaci pliku shapefile, a następnie korzystając z biblioteki rgdal po uzupełnieniu o projekcję otrzymujemy ramkę danych stanowiącą granice województw.

library(rgdal)
wojewodztwa <- readOGR("Dane/ITPOS/Shp", "Wojewodztwa")
wojewodztwa <- spTransform(wojewodztwa, CRS("+init=epsg:4326"))
wojewodztwa <- fortify(wojewodztwa, region = "JPT_KOD_JE")

Stąd teraz możemy przystąpić do wyrysowania odpowiedniej mapki z lokalizacją ITPOS z wyróżnieniem tego czy dana ITPOŚ udzieliła odpowiedzi na zadane pytania czy też nie. W tym celu skorzystam z biblioteki ggplot2. Lokalizacja zostanie przedstawiona w postaci punktów, gdzie zmienna colour została zmapowana na odpowiedź z wysłanej ankiety zapisanej w zmienna ankieta

ggplot() +
  geom_polygon(data=wojewodztwa,
               aes(long, lat, group=group),
               fill="white", color="gray") +
  geom_point(data=miejsca_df,
             aes(long, lat, color=ankieta), size = 5) +
  geom_text(data=miejsca_df,
            aes(long, lat, label=adres),
            size=4, color="black", hjust=-0.2) +
  theme_void()

W rezultacie powyższych działań został przygotowany slajd do prezentacji. Organizatorom bardzo dziękuję za zaproszenie.

Dodatkowe materiały

• Informacje o konferencji Innowacje w inżynierii środowiska połączonej z jubileuszem 25-lecia pracy naukowej prof. Hanny Obarskiej-Pempkowiak
• Zmiany w korzystaniu z Google Maps API
• Stosowanie geom_point w rysunkach typu scatterplot