Jak wszyscy wiedzą, z mediów społecznościowych można się o każdym dowiedzieć wszystkiego. Powszechny dostęp do wielu z tych serwisów został w ciągu ostatniego roku mocno ograniczony (np. Facebook) albo też obwarowany zdobyciem serii kodów dostępu etc. Na potrzeby dzisiejszych zajęć zajmiemy się Twitterem; dostęp do niego jest możliwy dzięki bibliotece rtweet, niestety, prawdopodobnie nie uda się jej uruchomić na kompuetrach w Pracowni... Twitter wymaga aż czterech różnych kodów, aby można było korzystać z API - te poniżej zostały stworzone specjalnie na potrzeby zajęć TEXT i po około tygodniu zostaną zresetowane.
# PRZYKŁAD 10.1
library(rtweet)
create_token(
app = "TEXTclass",
consumer_key = "mjBZhJWhlcYYkjTNQX1zC9z1t",
consumer_secret = "J8KjngzwuimbpanNzGuWJOhx6flNID5D9wyhguCpxI5HMzjWTt",
access_token = "414634084-T9H7IJR7nX0k2wYuJE8SWydgyaPduV4L0CweL3M1",
access_secret = "EKQBhBnAf6Ew9Q9adeY1bn7YJEykKdsYpUbTUuyknp4ZW")## <Token>
## <oauth_endpoint>
## request: https://api.twitter.com/oauth/request_token
## authorize: https://api.twitter.com/oauth/authenticate
## access: https://api.twitter.com/oauth/access_token
## <oauth_app> TEXTclass
## key: mjBZhJWhlcYYkjTNQX1zC9z1t
## secret: <hidden>
## <credentials> oauth_token, oauth_token_secret
## ---Po tym uwierzytelnaniu można już popbrać dane - poniżej zbierzemy 10000 tweetów w konta CNN
# PRZYKŁAD 10.2
tw <- get_timeline("cnn", n = 10000)
tw## # A tibble: 3,209 x 90
## user_id status_id created_at screen_name text source
## <chr> <chr> <dttm> <chr> <chr> <chr>
## 1 759251 12163519... 2020-01-12 13:31:09 CNN Form... Socia...
## 2 759251 12163443... 2020-01-12 13:01:06 CNN Two ... Socia...
## 3 759251 12163405... 2020-01-12 12:46:08 CNN A se... Socia...
## 4 759251 12163368... 2020-01-12 12:31:10 CNN For ... Socia...
## 5 759251 12163329... 2020-01-12 12:15:46 CNN A ne... Socia...
## 6 759251 12163311... 2020-01-12 12:08:43 CNN Mill... Socia...
## 7 759251 12163292... 2020-01-12 12:01:05 CNN A nu... Socia...
## 8 759251 12163277... 2020-01-12 11:55:00 CNN "Ste... Socia...
## 9 759251 12163254... 2020-01-12 11:46:05 CNN The ... Socia...
## 10 759251 12163232... 2020-01-12 11:37:16 CNN Volu... Socia...
## # ... with 3,199 more rows, and 84 more variables: display_text_width <dbl>,
## # reply_to_status_id <chr>, reply_to_user_id <chr>,
## # reply_to_screen_name <chr>, is_quote <lgl>, is_retweet <lgl>,
## # favorite_count <int>, retweet_count <int>, quote_count <int>,
## # reply_count <int>, hashtags <list>, symbols <list>, urls_url <list>,
## # urls_t.co <list>, urls_expanded_url <list>, media_url <list>,
## # media_t.co <list>, media_expanded_url <list>, media_type <list>,
## # ext_media_url <list>, ext_media_t.co <list>, ext_media_expanded_url <list>,
## # ext_media_type <chr>, mentions_user_id <list>, mentions_screen_name <list>,
## # lang <chr>, quoted_status_id <chr>, quoted_text <chr>,
## # quoted_created_at <dttm>, quoted_source <chr>, quoted_favorite_count <int>,
## # quoted_retweet_count <int>, quoted_user_id <chr>, quoted_screen_name <chr>,
## # quoted_name <chr>, quoted_followers_count <int>,
## # quoted_friends_count <int>, quoted_statuses_count <int>,
## # quoted_location <chr>, quoted_description <chr>, quoted_verified <lgl>,
## # retweet_status_id <chr>, retweet_text <chr>, retweet_created_at <dttm>,
## # retweet_source <chr>, retweet_favorite_count <int>,
## # retweet_retweet_count <int>, retweet_user_id <chr>,
## # retweet_screen_name <chr>, retweet_name <chr>,
## # retweet_followers_count <int>, retweet_friends_count <int>,
## # retweet_statuses_count <int>, retweet_location <chr>,
## # retweet_description <chr>, retweet_verified <lgl>, place_url <chr>,
## # place_name <chr>, place_full_name <chr>, place_type <chr>, country <chr>,
## # country_code <chr>, geo_coords <list>, coords_coords <list>,
## # bbox_coords <list>, status_url <chr>, name <chr>, location <chr>,
## # description <chr>, url <chr>, protected <lgl>, followers_count <int>,
## # friends_count <int>, listed_count <int>, statuses_count <int>,
## # favourites_count <int>, account_created_at <dttm>, verified <lgl>,
## # profile_url <chr>, profile_expanded_url <chr>, account_lang <lgl>,
## # profile_banner_url <chr>, profile_background_url <chr>,
## # profile_image_url <chr>Pakiet rtweet ma wbudowaną funkcję ts_plot(), która korzysta z pakietu ggplot2 i do której, podobnie jak do ggplota, można dodawać kolejne warstwy. Argumentem funkcji ts_plot jest rozdzielczość z jaką chcemy wykreślić aktywnośc (czyli liczbę tweetów) w danych.
# PRZYKŁAD 10.3
ts_plot(tw, "5 hour", color = "gray") + geom_point(color = "blue")
Biblioteka oferuje całą masę różnych opcji (polecam prezentację oraz dokumentację). My zrobimy tylko jeden prosty wykres, ale wykorzystamy do niego wcześniej zdobyte metody - poliyczymy walencję w kolejnych tweetach, uzywając klasyfikatora słownikowego, opartego na zbiorze Warriner et al. z Laboratorium 6.
# PRZYKŁAD 10.4
library(tidytext)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionemo <- as_tibble(read.csv("http://www.fizyka.pw.edu.pl/~julas/TEXT/lab/Ratings_Warriner_et_al.csv", stringsAsFactors = F))
emo <- emo %>%
select(word = Word, valence = V.Mean.Sum, arousal = A.Mean.Sum)
tw.data <- tw %>% transmute(id = row_number(), text = text) %>% unnest_tokens(word, text)
tw.sent <- tw.data %>%
inner_join(emo) %>%
group_by(id) %>%
summarise(valence = mean(valence))## Joining, by = "word"tw$valence[tw.sent$id] <- tw.sent$valence## Warning: Unknown or uninitialised column: 'valence'.ggplot(tw, aes(x = created_at, y = valence)) +
geom_line(alpha=0.25, color = "red") +
geom_smooth(method = "loess", span = 0.01)
Google Maps
Google jest prawdziwym potentatem i w zasadzie monopolistą jeśli chodzi o szereg usług - możemy albo z nich nie korzystać albo też zgadzać się na różne praktyki. Nie da się jednak ukryć, że poniższe usługi, do których dedykowane są konkretne API są dość przydatne:
- Distance Matrix,
- Maps Static,
- Directions,
- Geocoding,
- Maps
Nas interesuje w tym momencie szczególnie ta ostatnia pozycja. R oferuje kilka pakietów do obsługi Google Maps - my skorzystamy z ggmap. Jak sama nazwa wskazuje, pakiet korzysta z ggplot2 i jest z nim kompatybilny pod względem idei użycia warstw. Na począwek wygenerujemy prostą mapę scentrowaną na środku Warszawy. Do obsługi map potrzebny jest nam klucz API - poniższy jest wygenerowany na użytek zajęć i zostanie zresetowany po dwóch tygodniach.
# PRZYKŁAD 10.6
library(ggmap)## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.## Please cite ggmap if you use it! See citation("ggmap") for details.register_google("AIzaSyCzA85Bwcm8Ott7ow28VpyXVRSlJomUdN0")
map <- get_googlemap(center = c(lon = 21.00, lat = 52.25), zoom = 11, maptype = "road")## Source : https://maps.googleapis.com/maps/api/staticmap?center=52.25,21&zoom=11&size=640x640&scale=2&maptype=roadmap&key=xxxggmap(map)
Na mapę naniesiemy teraz pozycje przystanków sieci transportu miejskiego - zbiór wejściowy zosta pobrany ze serwera ZTM.
# PRZYKŁAD 10.7
stops <- readLines("http://www.if.pw.edu.pl/~julas/TEXT/lab/stops.txt")
head(stops)## [1] " 1001 KIJOWSKA, -- WARSZAWA"
## [2] " *PR 9"
## [3] " 100101 2 Ul./Pl.: TARGOWA, Kier.: AL.ZIELENIECKA, Y= 52.248670 X= 21.044260 "
## [4] " L 6 - stały: 125 135 138 166 509 517 "
## [5] " L 1 - na żądanie: N21 "
## [6] " 100102 2 Ul./Pl.: TARGOWA, Kier.: ZĄBKOWSKA, Y= 52.249020 X= 21.044540 "Na oko zbiór ma dość nieuporządkowaną postać, na szczęście okazuje się, że poszczególne dane, które chcielibyśmy na nasze potrzeby otrzymać (numer przystanku oraz współrzędne) łatwo poddają się ekstrakcji za pomocą prostych wyrażeń regularnych
# PRZYKŁAD 10.8
get.coor <- function(line) {
p <- regexpr("[0-9]{6}", line)
y <- regexpr("Y= 52.[0-9]{5,6}", line)
x <- regexpr("X= 2[0-1].[0-9]{5,6}", line)
p <- substr(line, p[1], p[1] + attr(p, "match.length") - 1)
y <- substr(line, y[1] + 3, y[1] + attr(y, "match.length") - 1)
x <- substr(line, x[1] + 3, x[1] + attr(x, "match.length") - 1)
data.frame(stop_id = as.numeric(p), long = as.numeric(y), lat = as.numeric(x))
}
stops <- get.coor(stops[grep("Y= ", stops)])
head(stops)## stop_id long lat
## 1 100101 52.24867 21.04426
## 2 100102 52.24902 21.04454
## 3 100103 52.24900 21.04398
## 4 100104 52.24990 21.04173
## 5 100105 52.25035 21.04386
## 6 100106 52.25001 21.04371Korzystając z powyższej ramki danych, możemy teraz nanieść poszczególne punkty na mapę, dodając je jako kolejną warstwę ggplota.
# PRZYKŁAD 10.9
ggmap(map) + geom_point(data = stops, aes(x = lat, y = long), size=0.5, shape = 19)## Warning: Removed 1286 rows containing missing values (geom_point).
Na bazie takich danych można się pokusić o np. naniesienie warstwy z interpolowaną gestością przystanków (warstwa stat_density2d). Przy okazji dokonamy również drobnej zmiany podstawowej mapy, ustawiając język na polski oraz kolor na czarno-biały.
# PRZYKŁAD 10.10
map <- get_googlemap(center = c(lon = 21.00, lat = 52.25), zoom = 11, maptype = "road", color = "bw", language = "PL")## Source : https://maps.googleapis.com/maps/api/staticmap?center=52.25,21&zoom=11&size=640x640&scale=2&maptype=roadmap&language=PL&key=xxxggmap(map) +
geom_point(data = stops, aes(x = lat, y = long), size=0.5, shape = 19) +
stat_density2d(data = stops, aes(x = lat, y = long, fill = stat(level)), geom = "polygon", alpha = 0.05, bins = 30) + scale_fill_gradient(low = "blue", high = "red")## Warning: Removed 1286 rows containing non-finite values (stat_density2d).## Warning: Removed 1286 rows containing missing values (geom_point).
Tramwaje w Warszawie
Pozycje przystanków są na pewno interesujące, ale zwykle więcej uwagi poświęca się danym dynamicznym, zmieniającym się w czasie. Takie możliwości oferuje np. serwis Otwarte Dane pod egidą Urzędu Miasta Warszawy. Rejestrując tam konto otrzymujemy klucz API do różnych danych - nas interesuje serwis udostępniający w czasie rzeczywistym pozycje tramajów w Warszawie.
Do opracowania wyników zapytań potrzebujemy biblioteki jsonlite, która obsługuje format JSON zwracany przez API
# PRZYKŁAD 10.11
library(jsonlite)##
## Attaching package: 'jsonlite'## The following object is masked from 'package:rtweet':
##
## flattenum.waw.api <- "9c9a80e9-68d1-4f74-8d49-c241f3f5649f"
um.waw.url <- "https://api.um.warszawa.pl/api/action/busestrams_get/?resource_id=f2e5503e-927d-4ad3-9500-4ab9e55deb59&apikey="
url.api <- paste(um.waw.url, um.waw.api, "&type=2", sep = "")
trams <- fromJSON(url.api)
head(trams$result)## Lines Lon VehicleNumber Time Lat Brigade
## 1 6 21.00723 1222 2020-01-12 14:56:54 52.26009 6
## 2 6 20.92970 1243 2020-01-12 14:56:53 52.29192 4
## 3 22 21.04993 1261+1262 2020-01-10 19:35:58 52.23782 9
## 4 6 20.95281 1267 2020-01-12 14:56:52 52.28741 2
## 5 28 20.94414 1277+1278 2020-01-12 14:56:54 52.27151 3
## 6 28 20.94396 1321+1322 2020-01-12 14:56:50 52.27134 2Teraz swobodnie wykorzystujemy naszą mapę Google, aby nanieść na nią pozycje tramwajów.
# PRZYKŁAD 10.12
map <- get_googlemap(center = c(lon = 21.00, lat = 52.25), zoom = 12, maptype = "road", color="bw", language = "PL")## Source : https://maps.googleapis.com/maps/api/staticmap?center=52.25,21&zoom=12&size=640x640&scale=2&maptype=roadmap&language=PL&key=xxxggmap(map) +
geom_point(data = trams$result, aes(x = Lon, y = Lat))## Warning: Removed 25 rows containing missing values (geom_point).
Mapę można trochę "uatrakcyjnić", nakładając na punkty numer linii oraz dodatkowo kolorując według tych numerów.
# PRZYKŁAD 10.13
trams$result$Lines <- as.numeric(gsub(" ", "", trams$result$Lines))
ggmap(map) +
geom_point(data = trams$result, aes(x = Lon, y = Lat, color = as.factor(Lines)), size = 4) +
geom_text(data = trams$result, aes(x = Lon, y = Lat, label = Lines), size = 3) +
scale_color_discrete(guide = FALSE)## Warning: Removed 25 rows containing missing values (geom_point).## Warning: Removed 25 rows containing missing values (geom_text).
Zanieczyszczenie powietrza
Jak wiadomo, zanieczyszczenie powietrza (czyli popularny smog) stanowi coraz większy problem w Polsce - jest to wypadkowa natężonego ruchu pojazdów oraz ogrzewania domów za pomocą kiespkiej jakości paliw, przy czym ten drugi czynnik jest raczej decydujący. W rezultacie takiej sytuacji, pojawiają się coraz to nowe przedsiębiorstwa, oferujące dostęp do metod moniotorowania poziomu zanieczyszczenia. Jedną z takich form jest Airly, która udostępnia mapę oraz apki, umożliwiające ocene stanu powietrza w danej chwili, jak również obejrzenie 24-godzinnej historii takich pomiarów oraz prognoze na kolejną dobę. Airly dostarcza również API: po zarejestrowaniu się, otrzymujemy klucz dający umożliwiający wykonanie maks. 50 zapytań na minutę oraz do 1000 na dobę.
Naszym pierwszym krokiem będzie zdobycie danych 45 czujników z obszaru Warszawy.
# PRZYKŁAD 10.14
airly.api = "0oRd9wGrqtnhN91H1UaAhPwKjRhGrmxy"
airly.url1 <- "https://airapi.airly.eu/v2/installations/nearest?lat=52.25&lng=21.00&maxDistanceKM=30&maxResults=45&apikey="
query <- paste(airly.url1, airly.api, sep = "")
czujniki <- fromJSON(query)
head(czujniki)## id location.latitude location.longitude address.country address.city
## 1 8941 52.24435 21.00421 Poland Warszawa
## 2 2191 52.25574 20.99337 Poland Warszawa
## 3 10238 52.25662 20.99424 Poland Warszawa
## 4 1074 52.24233 20.99543 Poland Warszawa
## 5 8882 52.24433 20.98954 Poland Warszawa
## 6 2350 52.25445 20.98454 Poland Warszawa
## address.street address.number address.displayAddress1
## 1 Antoniego Corazziego 4 Warszawa
## 2 Inflancka 4b Warszawa
## 3 Inflancka 4 Warszawa
## 4 aleja "Solidarności" 113 Warszawa
## 5 Nowolipki 16 Warszawa
## 6 aleja Jana Pawła II 80 Warszawa
## address.displayAddress2 elevation airly sponsor.id sponsor.name
## 1 Antoniego Corazziego 111.70 TRUE 511 eAzymut.pl
## 2 Inflancka 4A 101.07 TRUE 22 Aviva
## 3 Inflancka 99.72 TRUE 582 KPMG
## 4 aleja Solidarności 113 114.64 TRUE 49 Anonimowy
## 5 Nowolipki 116.42 TRUE 519 ZGN Wola
## 6 Jana Pawła II 80 108.31 TRUE 22 Aviva
## sponsor.description
## 1 Airly Sensor's sponsor
## 2 Airly Sensor's sponsor
## 3 Airly Sensor's sponsor
## 4 Airly Sensor's sponsor
## 5 Airly Sensor's sponsor
## 6 Airly Sensor's sponsor
## sponsor.logo
## 1 https://cdn.airly.eu/logo/eAzymutpl_1558097164668_807603640.jpg
## 2 https://cdn.airly.eu/logo/Aviva_1538146740542_399306786.jpg
## 3 https://cdn.airly.eu/logo/KPMG_1571916912072_123401498.jpg
## 4 https://cdn.airly.eu/logo/Ambasador_Airly.jpg
## 5 https://cdn.airly.eu/logo/ZGNWola_1559737411624_2068990889.jpg
## 6 https://cdn.airly.eu/logo/Aviva_1538146740542_399306786.jpg
## sponsor.link
## 1 https://www.eazymut.pl/
## 2 https://wiemczymoddycham.pl/
## 3 <NA>
## 4 <NA>
## 5 <NA>
## 6 https://wiemczymoddycham.pl/Jak widać, każdy z czujników ma unikalny numer (id) oraz podaną lokalizajcę. Będziemy teraz chcieli odpytać każdy z czujników - naszym celem jest stan powietrza w danej chwili (reprezentowany przez tzn indeks ACQ), ale przy okazji otrzymamy też całą masę innych danych, które wykorzystamy później. W tym momencie stworzymy sobie ramkę danych, zawierającą wartości indeksu oraz pozycje czujników, a następnie korzystając z ggmap naniesiemy te dane na mapę.
# PRZYKŁAD 10.15
get.air.mes <- function(id, api.key) {
query <- sprintf("https://airapi.airly.eu/v2/measurements/installation?installationId=%d&apikey=%s", id, api.key)
print(id)
#Sys.sleep(1)
fromJSON(query)
}
pomiary <- lapply(czujniki$id, get.air.mes, api.key = airly.api)## [1] 8941
## [1] 2191
## [1] 10238
## [1] 1074
## [1] 8882
## [1] 2350
## [1] 3190
## [1] 8883
## [1] 8824
## [1] 9928
## [1] 9932
## [1] 337
## [1] 658
## [1] 2724
## [1] 9933
## [1] 8881
## [1] 6532
## [1] 7741
## [1] 7755
## [1] 6864
## [1] 9944
## [1] 9064
## [1] 6961
## [1] 26
## [1] 3410
## [1] 8880
## [1] 3507
## [1] 8879
## [1] 6959
## [1] 2679
## [1] 7602
## [1] 2330
## [1] 9112
## [1] 9939
## [1] 8931
## [1] 817
## [1] 8174
## [1] 6630
## [1] 2720
## [1] 3328
## [1] 27
## [1] 7603
## [1] 7702
## [1] 2143
## [1] 9931indexes <- sapply(1:length(pomiary), function(i) pomiary[[i]]$current$indexes$value)
data <- cbind(czujniki$location, indexes)
head(data)## latitude longitude indexes
## 1 52.24435 21.00421 50.67
## 2 52.25574 20.99337 46.08
## 3 52.25662 20.99424 39.50
## 4 52.24233 20.99543 49.52
## 5 52.24433 20.98954 44.26
## 6 52.25445 20.98454 44.74ggmap(map) + geom_point(data = data, aes(x = longitude, y = latitude, color = indexes), size = 4) +
scale_color_gradient2("AQI", low = "green", high = "violet", mid = "red", midpoint = 85)
Jak wspomiano, nie są to jedyne dane, które przechowujemy w zmiennej pomiary. Dla każdego czujnika mamy dostęp do hisorii oraz do prognozy: wykorzystajmy tę pierwszą informację do policzenia korelacji pomiędzy pomiarami poszczególnych czujników.
# PRZYKŁAD 10.16
library(corrplot)## corrplot 0.84 loadedhistoria <- sapply(1:length(pomiary), function(i) sapply(pomiary[[i]]$history$indexes, function(x) x$value))
C <- cor(historia)
corrplot(C, type = "upper", na.label = "-", diag = F)
Jak widać, korelacje są raczej wysokie - tylko jeden czujnik (nr 11) jest wyraźnie inny. Dysponując takimi danymi można się pokusić o sprawdzenie jak korelacja zachowuje się w funkcji odległości. Potrzebujemy do tego funkcji wyznaczającej odegłości pomiędzy punktami określonami przez współrzędne geograficzne (wykorzystamy bibliotekę geosphere). Poza tym użyjemy bibliotek fields i Hmisc do binowania oraz prezentowania słupków niepewności.
# PRZYKŁAD 10.17
library(geosphere, quietly = T, warn.conflicts = F)
library(fields, quietly = T, warn.conflicts = F)## Spam version 2.5-1 (2019-12-12) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.##
## Attaching package: 'spam'## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve## See https://github.com/NCAR/Fields for
## an extensive vignette, other supplements and source codelibrary(Hmisc, quietly = T, warn.conflicts = F)
D <- distm(data[c(2, 1)], fun = distHaversine)/1000
d <- stats.bin(D[upper.tri(D)], C[upper.tri(C)])
x <- d$centers
y <- d$stats[2,]
yerr <- d$stats[3] / sqrt(d$stats[1,])
par(mfrow = c(1,2))
plot(D[upper.tri(D)], C[upper.tri(C)], pch = 19, cex = 0.5, col = "gray", tcl = 0.25, xlab = "odelglosc [km]", ylab = "korelacja")
errbar(x, y, y + yerr, y - yerr, add = T, col = "blue", errbar.col = "blue")
errbar(x, y, y + yerr, y - yerr, col = "blue", errbar.col = "blue", tcl = 0.25, xlab = "odelglosc [km]", ylab = "korelacja")
par(mfrow = c(1,1))