Financial Times zamieścił wykres wględnego tempa wzrostu (rate of growth) czyli procentu liczonego jako liczba-nowych / liczba-ogółem-z-okresu-poprzedniego x 100%. Na wykresie wględnego tempa wzrostu zachorowań na COVID19 wszystkim spada: Every day the Covid-19 virus is infecting an increasing number of people, but the rate of growth in cases in some of the worst-hit countries is starting to slow. Powyższe Czerscy przetłumaczyli jako m.in. trend dotyczy niemal wszystkich krajów rozwiniętych. [he, he... Rozwiniętych pod względem liczby chorych, pewnie chcieli uściślić, ale się nie zmieściło]
Spróbowałem narysować taki wykres samodzielnie:
library("dplyr") library("ggplot2") library("ggpubr") ## surl <- "https://www.ecdc.europa.eu/en/publications-data/download-todays-data-geographic-distribution-covid-19-cases-worldwide" today <- Sys.Date() tt<- format(today, "%d/%m/%Y") #d <- read.csv("covid19_C.csv", sep = ';', header=T, na.string="NA", stringsAsFactors=FALSE); d <- read.csv("covid19_C.csv", sep = ';', header=T, na.string="NA", colClasses = c('factor', 'factor', 'factor', 'character', 'character', 'numeric', 'numeric')); d$newc <- as.numeric(d$newc) d$newd <- as.numeric(d$newd)
Zwykłe read_csv
skutkowało tym, że newc/newd
nie były liczbami całkowitymi, tylko
czynnikami. Z kolei dodanie colClasses kończyło się błędem. W końcu stanęło na tym, że czytam
dane w kolumnach newc/newd zadeklarowanych jako napisy a potem konwertuję na liczby. Czy to jest
prawidłowa strategia to ja nie wiem...
Kolejny problem: kolumny newc/newd zawierają NA, wykorzystywana
później funkcja cumsum
z pakietu dplyr
,
obliczająca szereg kumulowany nie działa poprawnie jeżeli szereg
zawiera NA
. Zamieniam od razu NA
na zero.
Alternatywnie można korzystać
z funkcji replace_na
(pakiet dplyr):
# change NA to 0 d[is.na(d)] = 0 # Alternatywnie replace_na #d %>% replace_na(list(newc = 0, newd=0)) %>% # mutate( cc = cumsum(newc), dd=cumsum(newd))
Ograniczam się tylko do danych dla wybranych krajów, nie starszych niż 16 luty 2020:
d <- d %>% filter(as.Date(date, format="%Y-%m-%d") > "2020-02-15") %>% as.data.frame str(d) last.obs <- last(d$date) c1 <- c('IT', 'DE', 'ES', 'UK', 'FR') d1 <- d %>% filter (id %in% c1) %>% as.data.frame str(d1)
Obliczam wartości skumulowane (d
zawiera już skumulowane wartości,
ale obliczone Perlem tak nawiasem mówiąc):
t1 <- d1 %>% group_by(id) %>% summarise(cc = sum(newc, na.rm=T), dd=sum(newd, na.rm=T)) t1c <- d %>% group_by(id) %>% mutate(cum_cc = cumsum(newc), cum_dd = cumsum(newd)) %>% filter (id %in% c1) %>% filter(as.Date(date, format="%Y-%m-%d") > "2020-02-15") %>% as.data.frame str(t1c)
Wykres wartości skumulowanych:
pc1c <- ggplot(t1c, aes(x= as.Date(date, format="%Y-%m-%d"), y=cum_cc)) + geom_line(aes(group = id, color = id), size=.8) + xlab(label="") + theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) + ggtitle(sprintf("COVID19: total confirmed cases (%s)", last.obs), subtitle=sprintf("%s", surl)) ggsave(plot=pc1c, "Covid19_1c.png", width=15)
Kolumny cum_lcc/cum_ldd
zawierają wartości z kolumny cum_cc/cum_dd
ale opóźnione o jeden okres (funkcja lag
):
## t1c <- t1c %>% group_by(id) %>% mutate(cum_lcc = lag(cum_cc)) %>% as.data.frame t1c <- t1c %>% group_by(id) %>% mutate(cum_ldd = lag(cum_dd)) %>% as.data.frame t1c$gr_cc <- t1c$newc / (t1c$cum_lcc + 0.01) * 100 t1c$gr_dd <- t1c$newd / (t1c$cum_ldd + 0.01) * 100 ## Początkowo wartości mogą być ogromne zatem ## zamień na NA jeżeli gr_cc/dd > 90 t1c$gr_cc[ (t1c$gr_cc > 90) ] <- NA t1c$gr_dd[ (t1c$gr_dd > 90) ] <- NA
Wykres tempa wzrostu:
pc1c_gr <- ggplot(t1c, aes(x= as.Date(date, format="%Y-%m-%d"), y=gr_cc, colour = id, group=id )) + ##geom_line(aes(group = id, color = id), size=.8) + geom_smooth(method = "loess", se=FALSE) + xlab(label="") + theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) + ggtitle(sprintf("COVID19: confirmed cases growth rate (smoothed)"), subtitle=sprintf("%s", surl)) ggsave(plot=pc1c_gr, "Covid19_1g.png", width=15)
To samo co wyżej tylko dla PL/CZ/SK/HU:
c2 <- c('PL', 'CZ', 'SK', 'HU') t2c <- d %>% group_by(id) %>% mutate(cum_cc = cumsum(newc), cum_dd = cumsum(newd)) %>% filter (id %in% c2) %>% filter(as.Date(date, format="%Y-%m-%d") > "2020-02-15") %>% as.data.frame ##str(t2c) t2c.PL <- t2c %>% filter (id == "PL") %>% as.data.frame t2c.PL head(t2c.PL, n=200) pc2c <- ggplot(t2c, aes(x= as.Date(date, format="%Y-%m-%d"), y=cum_cc)) + geom_line(aes(group = id, color = id), size=.8) + xlab(label="") + theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) + ggtitle(sprintf("COVID19: total confirmed cases (%s)", last.obs), subtitle=sprintf("Total: %s\n%s", lab1c, surl)) ggsave(plot=pc2c, "Covid19_2c.png", width=15) t2c <- t2c %>% group_by(id) %>% mutate(cum_lcc = lag(cum_cc)) %>% as.data.frame t2c <- t2c %>% group_by(id) %>% mutate(cum_ldd = lag(cum_dd)) %>% as.data.frame t2c$gr_cc <- t2c$newc / (t2c$cum_lcc + 0.01) * 100 t2c$gr_dd <- t2c$newd / (t2c$cum_ldd + 0.01) * 100 ## zamień na NA jeżeli gr_cc/dd > 90 t2c$gr_cc[ (t2c$gr_cc > 90) ] <- NA t2c$gr_dd[ (t2c$gr_dd > 90) ] <- NA t2c.PL <- t2c %>% filter (id == "PL") %>% as.data.frame t2c.PL pc2c_gr <- ggplot(t2c, aes(x= as.Date(date, format="%Y-%m-%d"), y=gr_cc, colour = id, group=id )) + ##geom_line(aes(group = id, color = id), size=.8) + geom_smooth(method = "loess", se=FALSE) + xlab(label="") + theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) + ggtitle(sprintf("COVID19: confirmed cases growth rate (smoothed)"), subtitle=sprintf("%s", surl)) ggsave(plot=pc2c_gr, "Covid19_2g.png", width=15)
Koniec
Trzeba coś robić w czasie kwarantanny
## https://b-rodrigues.github.io/modern_R/ ## https://gist.github.com/imartinezl/2dc230f33604d5fb729fa139535cd0b3 library("eurostat") library("dplyr") library("ggplot2") library("ggpubr") ## options(scipen=1000000) dformat <- "%Y-%m-%d" eu28 <- c("AT", "BE", "BG", "HR", "CY", "CZ", "DK", "EE", "FI", "FR", "DE", "EL", "HU", "IE", "IT", "LT", "LU", "LV", "MT", "NL", "PL", "PT", "RO", "SK", "SI", "ES", "SE") eu6 <- c("DE", "FR", "IT", "ES", "PL") ### Demo_mor/ Mortality monthly ### ### ### dm <- get_eurostat(id="demo_mmonth", time_format = "num"); dm$date <- sprintf ("%s-%s-01", dm$time, substr(dm$month, 2, 3)) str(dm) ## There are 12 moths + TOTAL + UNKN dm_month <- levels(dm$month) dm_month ## Only new data dm28 <- dm %>% filter (geo %in% eu28 & as.Date(date) > "1999-12-31") str(dm28) levels(dm28$geo) ## Limit to DE/FR/IT/ES/PL: dm6 <- dm28 %>% filter (geo %in% eu6) str(dm6) levels(dm6$geo) pd1 <- ggplot(dm6, aes(x= as.Date(date, format="%Y-%m-%d"), y=values)) + geom_line(aes(group = geo, color = geo), size=.4) + xlab(label="") + ##scale_x_date(date_breaks = "3 months", date_labels = "%y%m") + scale_x_date(date_breaks = "6 months", date_labels = "%m\n%y", position="bottom") + theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) + ggtitle("Deaths", subtitle="https://ec.europa.eu/eurostat/data/database (demo_mmonth)") ## Newer data dm6 <- dm6 %>% filter (as.Date(date) > "2009-12-31") pd2 <- ggplot(dm6, aes(x= as.Date(date, format="%Y-%m-%d"), y=values)) + geom_line(aes(group = geo, color = geo), size=.4) + xlab(label="") + scale_x_date(date_breaks = "3 months", date_labels = "%m\n%y", position="bottom") + theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) + ggtitle("Deaths", subtitle="https://ec.europa.eu/eurostat/data/database (demo_mmonth)") ggsave(plot=pd1, file="mort_eu_L.png", width=12) ggsave(plot=pd2, file="mort_eu_S.png", width=12)
## Live births (demo_fmonth) ### ### ### dm <- get_eurostat(id="demo_fmonth", time_format = "num"); dm$date <- sprintf ("%s-%s-01", dm$time, substr(dm$month, 2, 3)) str(dm) ## There are 12 moths + TOTAL + UNKN dm_month <- levels(dm$month) dm_month dm28 <- dm %>% filter (geo %in% eu28 & as.Date(date) > "1999-12-31") str(dm28) levels(dm28$geo) dm6 <- dm28 %>% filter (geo %in% eu6) str(dm6) levels(dm6$geo) pd1 <- ggplot(dm6, aes(x= as.Date(date, format="%Y-%m-%d"), y=values)) + geom_line(aes(group = geo, color = geo), size=.4) + xlab(label="") + ##scale_x_date(date_breaks = "3 months", date_labels = "%y%m") + scale_x_date(date_breaks = "6 months", date_labels = "%m\n%y", position="bottom") + theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) + ggtitle("Births", subtitle="https://ec.europa.eu/eurostat/data/database (demo_fmonth)") ## dm6 <- dm6 %>% filter (as.Date(date) > "2009-12-31") pd2 <- ggplot(dm6, aes(x= as.Date(date, format="%Y-%m-%d"), y=values)) + geom_line(aes(group = geo, color = geo), size=.4) + xlab(label="") + scale_x_date(date_breaks = "3 months", date_labels = "%m\n%y", position="bottom") + theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) + ggtitle("Births", subtitle="https://ec.europa.eu/eurostat/data/database (demo_fmonth)") ggsave(plot=pd1, file="birt_eu_L.png", width=12) ggsave(plot=pd2, file="birt_eu_S.png", width=12)
## Population (only) yearly ### ### ## ## Population change - Demographic balance and crude rates at national level (demo_gind) dp <- get_eurostat(id="demo_gind", time_format = "num"); dp$date <- sprintf ("%s-01-01", dp$time) str(dp) dp_indic_dic <- get_eurostat_dic("indic_de") dp_indic_dic dp28 <- dp %>% filter (geo %in% eu28 & time > 1999 & indic_de == "JAN") str(dp28) dp6 <- dp28 %>% filter (geo %in% eu6) pdp1 <- ggplot(dp6, aes(x= as.Date(date, format="%Y-%m-%d"), y=values)) + geom_line(aes(group = geo, color = geo), size=.4) + xlab(label="") + ##scale_x_date(date_breaks = "3 months", date_labels = "%y%m") + ##scale_x_date(date_breaks = "6 months", date_labels = "%m\n%y", position="bottom") + theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) + ggtitle("Population", subtitle="https://ec.europa.eu/eurostat/data/database (demo_fmonth)") ggsave(plot=pdp1, file="pdp1", width=12)
Na stronie http://policja.pl/pol/form/1,Informacja-dzienna.html
udostępniane są dzienne dane dotyczące liczby interwencji,
zatrzymanych na gorącym uczynku, zatrzymanych poszukiwanych, pijanych
kierujących, wypadków, zabitych w wypadkach, rannych w wypadkach.
Ściągam wszystkie dane:
#!/bin/bash rm pp.html for ((i=0; i <= 274; i++)) do if [ ! -f ${i}.html ] ; then curl -o ${i}.html "http://policja.pl/pol/form/1,Informacja-dzienna.html?page=${i}" ; grep 'data-label' ${i}.html >> pp.html sleep 6 else grep 'data-label' ${i}.html >> pp.html echo done fi done
zamieniam prostymi skryptami na plik CSV, który ma następującą strukturę:
data;interwencje;zng;zp;znk;wypadki;zabici;ranni 2008-12-01;NA;873;344;447;135;1;1
okazuje się że liczba interwencji jest podawana od roku 2018, wcześniej nie była. Nic to wstawiamy NA.
Na przyszłość dane będą aktualizowane w ten sposób, że codziennie
(przez odpowiedni wpis w pliku crontab)
będzie pobierany plik http://policja.pl/pol/form/1,Informacja-dzienna.html
:
#!/usr/bin/perl use LWP::Simple; $PP="http://policja.pl/pol/form/1,Informacja-dzienna.html"; $PPBase="pp.csv"; $content = get("$PP"); $content =~ s/\r//g; # dla pewności usuń @content = split (/\n/, $content); foreach (@content) { chomp(); unless ($_ =~ m/data-label=/ ) { next } if ($_ =~ m/Data statystyki/ ) { $d = clean($_); } elsif ($_ =~ m/Interwencje/ ) { $i = clean($_); } elsif ($_ =~ m/Zatrzymani na g/ ) { $zg = clean($_); } elsif ($_ =~ m/Zatrzymani p/ ) { $zp = clean($_); } elsif ($_ =~ m/Zatrzymani n/ ) { $zn = clean($_); } elsif ($_ =~ m/Wypadki d/ ) { $w = clean($_); } elsif ($_ =~ m/Zabici/ ) { $z = clean($_); } elsif ($_ =~ m/Ranni/ ) { $r = clean($_); $l = "$d;$i;$zg;$zp;$zn;$w;$z;$r"; $last_line = "$l"; $last_date = "$d"; ## pierwszy wpis powinien zawierać dane dotyczące kolejnego dnia ## więc po pobraniu pierwszego można zakończyć last; } } ### read the database open (PP, "<$PPBase") || die "cannot open $PPBase/r!\n" ; while (<PP>) { chomp(); $line = $_; @tmp = split /;/, $line; } close(PP); ### append the database (if new record) open (PP, ">>$PPBase") || die "cannot open $PPBase/w!\n" ; unless ("$tmp[0]" eq "$last_date") { print PP "$last_line\n" } else {print STDERR "nic nowego nie widzę!\n"} close(PP); sub clean { my $s = shift; $s =~ s/<[^<>]*>//g; $s =~ s/[ \t]//g; return ($s); }
Zaktualizowana baza jest wysyłana na githuba. Tutaj jest:
https://github.com/hrpunio/Nafisa/tree/master/PP
Agregacja do danych tygodniowych okazała się nietrywialna
Niektóra lata zaczynają się od tygodnia numer 0 a inne od 1. Okazuje się,
że tak ma być (https://en.wikipedia.org/wiki/ISO_week_date#First_week
):
If 1 January is on a Monday, Tuesday, Wednesday or Thursday, it is in W01. If it is on a Friday, it is part of W53 of the previous year. If it is on a Saturday, it is part of the last week of the previous year which is numbered W52 in a common year and W53 in a leap year. If it is on a Sunday, it is part of W52 of the previous year.
Nie bawię się w subtelności tylko tygodnie o numerze zero dodaję do tygodnia z poprzedniego roku.
Sprawdzam czy jest OK i się okazuje że niektóre tygodnie mają 8 dni. W plikach html są błędy:
Błędne daty 2019-10-30 winno być 2019-09-30; podobnie błędne 2019-03-28 (winno być 2019-02-28), 2018-11-01 (2018-12-01), 2018-12-01 (2017-12-01), 2016-04-30 (2016-03-30), 2009-08-31 (2009-07-31). Powtórzone daty: 2016-03-10, 2010-07-25, 2010-01-10 (zdublowane/różne/arbitralnie usuwamy drugi) Ponadto brak danych z następujących dni: 2015-12-04--2015-12-07, 2015-04-17--2015-04-20, 2014-10-02--2014-10-05, 2014-01-23 i jeszcze paru innych (nie chcialo mi się poprawiać starych.)
Teraz jest OK, plik ppw.csv
ma nast strukturę:
rok;nrt;interwencje;in;zng;zngn;zp;zpn;znk;znkn;wypadki;wn;zabici;zn;ranni;rn;d1;d7
coś co się kończy na `n' to liczba tego co jest w kolumnie poprzedniej, np zn
to
liczba dni tygodnia dla kolumny zabici. Generalnie kolumny kończące się na `n'
zawierają 7 :-)
Kolumna d1
to pierwszy dzień tygodnia a kolumna d7
ostatni.
maxY <- max (d$zabici) pz <- ggplot(d, aes(x= as.factor(nrt), y=zabici )) + geom_bar(fill="steelblue", stat="identity") + xlab(label="") + ggtitle("Wypadki/zabici (Polska/2020)", subtitle="policja.pl/pol/form/1,Informacja-dzienna.html")
W sumie agregacja jest niepotrzebna, bo można ją zrobić na poziomie R
używając funkcji stat_summary
:
pw <- ggplot(d, aes(x= week, y=wypadki)) + stat_summary(fun.y = sum, geom="bar", fill="steelblue") + scale_x_date( labels = date_format("%y/%m"), breaks = "2 months") + xlab(label="") + #theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) + ggtitle("Wypadki (Polska/2018--2020)", subtitle="policja.pl/pol/form/1,Informacja-dzienna.html")
albo najpierw agregując dane a potem wykreślając wykres szeregu zagregowanego. Drugi sposób pozwala na przykład na dodanie linii oznaczających poziomy zagregowanego zjawiska/etykiety słupków w sposób `inteligentny'. Dodajemy etykiety (z numerem tygodnia) tylko dla słupków poniżej/powyżej Q1/Q3:
## agregowanie do danych tygodniowych kolumn ranni, zabici, wypadki dw <- d %>% group_by ( YrWeek) %>% summarise_at ( vars(ranni,zabici,wypadki), sum ) ## Obliczanie mediany i kwartyli median.zw <- median(dw$zabici) q1.zw <- quantile(dw$zabici, 1/4) q3.zw <- quantile(dw$zabici, 3/4) ## YrWeekZZ to numer tygodnia, dla tygodni w których liczba wypadków jest typowa ## numer jest pusty (żeby nie był drukowany); taki trick może jest lepszy dw$YrWeekZZ <- substr(dw$YrWeek,4,5) dw$YrWeekZZ[ (dw$zabici > q1.zw) & (dw$zabici < q3.zw) ] <- "" pz.2 <- ggplot(dw, aes(x= YrWeek, y=zabici)) + geom_bar(stat="identity", fill = "steelblue") + geom_text(data=dw, aes(label=sprintf("%s", YrWeekZZ), x=YrWeek, y= zabici), vjust=-0.9, size=3 ) + geom_hline(yintercept=median.zw, linetype="solid", color = "violet", size=.4) + geom_hline(yintercept=q3.zw, linetype="solid", color = "red", size=.4) + geom_hline(yintercept=q1.zw, linetype="solid", color = "red", size=.4) + xlab(label="rok/tydzień") + ylab(label="zabici") + scale_x_discrete(breaks=c("18/01", "18/10", "18/20", "18/30", "18/40", "19/01", "19/10", "19/20", "19/30", "19/40", "20/01", "20/10")) + # labels = c("/18/01", "18/10", "18/20", "")) ## tutaj niepotrzebne ggtitle("Wypadki/zabici (Polska/2018--2020)", subtitle="Linie poziomie: q1/me/q3 (źródło: policja.pl/pol/form/1,Informacja-dzienna.html)")
Stare Pi usiłuję zamienić na kamerę w ogródku:
## co to za wersja Pi? less /proc/cpuinfo ... Model : Raspberry Pi Model A Rev 2
Problemem jest zasięg WiFi (od routera do kamery jest jakieś 8 metrów przez szybę, żadnych murów). Do tego Pi ma tylko jedno złącze USB co utrudnia sprawę, bo ujawnia się złośliwość przedmiotów martwych (ZPM): z dwóch moich rezerwowych hubów żaden nie działa z tym konkretnym RPi. Podłączam hub z PCeta--ten działa...
Żeby było lepiej (z łącznością) zakupiłem TP-Link TL-WN722N (v2), ale okazało się że akurat wersja 2 nie jest rozpoznawana przez fabryczny Raspbian. Pech albo ZPM2.
Najpierw usiłowałem zainstalować stosowny sterownik
wg wskazówek ze strony https://github.com/lwfinger/rtl8188eu
:
git clone https://github.com/lwfinger/rtl8188eu.git sudo apt-get install raspberrypi-kernel-headers sudo ln -s /usr/src/linux-headers-$(uname -r) /lib/modules/$(uname -r)/build cd rtl8188eu/ make all sudo make install sudo reboot
Nie działa dalej, a moduł się skompilował, zainstalował i nawet jest ładowany...
Zadziałał ten przepis (https://www.raspberrypi.org/forums/viewtopic.php?t=250911#p1532103
):
You can download the driver from http://downloads.fars-robotics.net/wifi-drivers/8188eu-drivers/
.
Choose the driver that matches the output of command uname -a
for the correct kernel version:
uname -a # Linux aisara 4.19.97+ #1294 Thu Jan 30 13:10:54 GMT 2020 armv6l GNU/Linux mkdir Temp && cd Temp wget http://downloads.fars-robotics.net/wifi-drivers/8188eu-drivers/8188eu-4.19.97-1294.tar.gz tar -zxvf 8188eu-4.19.97-1294.tar.gz ./install.sh
BTW Installing my driver will disable the built in r8188eu driver. To re-enable the built in driver you will need to run the following commands:
sudo rm /etc/modprobe.d/8188eu.conf sudo rm /lib/modules/4.19.69-v7l+/kernel/drivers/net/wireless/8188eu.ko sudo depmod -a
Miałem już ze starych czasów puszkę z Raspberry w środku, z kamerą przyczepioną do ścianki, zrobioną dziurą na obiektyw i kablem doprowadzającym zasilanie przez piny P2+P6 (też tak można). Zasilacz w domu, poza oknem już 5V. Na wszelki wypadek -- żeby mi kogoś nie zabiło przez przypadek (a nawet czegoś, bo teraz czasy takie że jakby dzika poraziło, to też by była afera) -- wolałem nie ciągnąć 230V z mieszkania. Dlatego też zasilanie jest po zwykłym dwużyłowym kablu a nie przez USB (bo tak mi się wydawało prościej).
Znowu coś nie tak od pierwszego strzału (ZPM3):
mmal: mmal_component_create_core: could not create component 'vc.ril.camera' (1)
Ze starych czasów pamiętam, że problem może być z połączeniem/uszkodzeniem (taśmy na przykład). Zmieniłem kamerę i taśmę na inną. Działa. Musiała być uszkodzone widocznie.
Kamerą będę fotografował mój blok od tyłu, że tak powiem. Z drzewkiem brzoskiwni na pierwszym planie. Zdjęcia wysyłał na Twittera co 3 godziny na przykład:
0 6,9,12,15,18,21 * * * /home/pi/bin/mk1photo.sh
Na zdjęciu widać czujnik ruchu, bo faktycznie takowy też dokleiłem do obudowy. Że niby miał zdjęcia robić jak ruch wykryje, ale ponieważ to w zasadzie nie działało, więc teraz ten czujnik wprawdzie jest, ale nie podłączony. Atrapa...
Dane pierwotne: Center for Systems Science and
Engineering (CSSE/Johns Hopkins University)
https://github.com/CSSEGISandData/COVID-19
(także słynna
wizualizacja:
https://gisanddata.maps.arcgis.com/apps/opsdashboard/index.html#/bda7594740fd40299423467b48e9ecf6
.)
European Centre for Disease Prevention and Control
https://ecdc.europa.eu/en/geographical-distribution-2019-ncov-cases
.
Dane agregowane od innych: Our World in
Data/Coronavirus Source Data/WHO Situation Reports
https://ourworldindata.org/coronavirus-source-data
.
Są też tzw. dane w czasie rzeczywistym:
https://worldometers.info/coronavirus/
, ale ich
wiarygodność jest podejrzana, bo w przeciwieństwie do tych wyżej
opisanych nie wiadomo jak są zbierane i/lub
skąd agregowane (więc nie ma klikalnego linku).
Na stronie
https://ourworldindata.org/coronavirus-source-data
są dane nt liczby przypadków/zgonów z powodu zarażenia wirusem covid19, których źródłem
są 'Raporty Sytuacyjne WHO'
(https://www.who.int/emergencies/diseases/novel-coronavirus-2019/situation-reports/
).
Raporty są w formacie PDF więc bezpośrednio nie można korzystać
z publikowanych tam danych.
No ale uprzejmi ludzie z ourworldindata.org
już te raporty zamienili na csv i są one gotowe do pobrania:
wget -N https://covid.ourworldindata.org/data/full_data.csv
Za pomocą prostych skryptów Perla modyfikuję plik full_data.csv
tak, żeby
poszczególne kolumny zawierały:
date;id;country;newc;newd;totalc;totald
(data, ISO-kod
kraju, nazwa-kraju, nowe-przypadki, nowe-zgony, wszystkie-przypadki,
wszystkie-zgony)
Dla wybranych krajów rysują wykresy liniowe (wykorzystując R):
library("dplyr") library("ggplot2") library("ggpubr") ## today <- Sys.Date() tt<- format(today, "%d/%m/%Y") d <- read.csv("covid19.csv", sep = ';', header=T, na.string="NA"); d <- d %>% filter(as.Date(date, format="%Y-%m-%d") > "2020-02-15") %>% as.data.frame ## c1 <- c('ITw', 'DEw', 'ESw', 'UKw', 'FRw', 'DKw', 'SEw') # date;id;country;newc;newd;totalc;totald d1 <- d %>% filter (id %in% c1) %>% as.data.frame t1 <- d1 %>% group_by(id) %>% summarise(cc = sum(newc, na.rm=T), dd=sum(newd, na.rm=T)) lab1c <- toString(paste (sep=" = ", t1$id, t1$cc)) lab1d <- toString(paste (sep=" = ", t1$id, t1$dd)) str(d1) pc1 <- ggplot(d1, aes(x= as.Date(date, format="%Y-%m-%d"), y=newc)) + geom_line(aes(group = id, color = id), size=.8) + xlab(label="") + theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) + ggtitle(sprintf("COVID19: new confirmed cases (%s)", tt), subtitle=sprintf("Total: %s\n%s", lab1c, surl)) pd1 <- ggplot(d1, aes(x= as.Date(date, format="%Y-%m-%d"), y=newd)) + geom_line(aes(group = id, color = id), size=.8) + xlab(label="") + theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) + ggtitle(sprintf ("COVID19: deaths (%s)", tt), subtitle=sprintf("Total: %s\n%s", lab1d, surl)) c2 <- c('PLw', 'CZw', 'SKw', 'HUw', 'ROw', 'BGw', 'ELw') d2 <- d %>% filter (id %in% c2) %>% as.data.frame t2 <- d2 %>% group_by(id) %>% summarise(cc = sum(newc, na.rm=T), dd=sum(newd, na.rm=T)) str(d2) lab2c <- toString(paste (sep=" = ", t2$id, t2$cc)) lab2d <- toString(paste (sep=" = ", t2$id, t2$dd)) pc2 <- ggplot(d2, aes(x= as.Date(date, format="%Y-%m-%d"), y=newc)) + geom_line(aes(group = id, color = id), size=.8) + theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) + xlab(label="") + ggtitle(sprintf("COVID19: new confirmed cases (%s)", tt), subtitle=sprintf("Total: %s\n%s", lab2c, surl)) pd2 <- ggplot(d2, aes(x= as.Date(date, format="%Y-%m-%d"), y=newd)) + geom_line(aes(group = id, color = id), size=.8) + theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) + xlab(label="") + scale_y_continuous(breaks=c(1,2,3,4,5,6,7,8,9)) + ggtitle(sprintf ("COVID19: deaths (%s)", tt), subtitle=sprintf("Total: %s\n%s", lab2d, surl)) p1 <- ggarrange(pc1,pd1, ncol=2, nrow=1) p2 <- ggarrange(pc2,pd2, ncol=2, nrow=1) ggsave(plot=p1, "Covid19_1w.png", width=15) ggsave(plot=p2, "Covid19_2w.png", width=15)
Zatem: Liczba przypadków/zgonów z powodu zarażenia wirusem covid19 na podstawie danych ourworldindata.org/WHO:
UE, która jest powszechnie krytykowana, że nic nie robi w sprawie, okazuje się że coś tobi -- też udostępnia jakieś dane w temacie (https://www.ecdc.europa.eu/en/publications-data/download-todays-data-geographic-distribution-covid-19-cases-worldwide). Wprawdzie nie jest określone skąd te dane pochodzą, ale sądząc po ich zawartości źródło jest to samo (WHO).
wget -N https://www.ecdc.europa.eu/sites/default/files/documents/COVID-19-geographic-disbtribution-worldwide-2020-03-15.xls\ -O covid19.csv
Ponieważ dane, że tak powiem, unijne są w formacie xls
zamieniam
je na csv
, wykorzystując
do tego LibreOffice:
## zamień wszystkie pliki z bieżącego katalogu na csv ze ; (59) jako znakiem separacji: soffice --convert-to csv:"Text - txt - csv (StarCalc)":59,,0,1,1 --outdir . *.xls
Rysuję wykresy liniowe zmodyfikowanym z dokładnością do pliku z danymi R-skryptem. Wyniki są prawie takie same. Może bym nawet nie zwrócił uwagi, że się różnią gdyby nie podejrzane załamanie liczby przypadków dla Włoch dla 15.03.2020 (z 2,5 tys na 90).
Drążąc temat wyrysowałem wykresy dla wybranych czterech krajów w dwóch wariantach danych
(dane z ourworldindata.org oznaczone literką w
). W szczególności i niestety
Włochy 15/3/2020 odnotowały ponad 3497 nowych przypadków
a nie 90 jak podano w bazie Unijnej. Są też mniejsze różnice w innych miejscach:
Wszystko to robione jest automatem co pobiera/zamienia/rysuje/wstawia na githuba
(https://github.com/hrpunio/Nafisa/tree/master/Covid19
)
oraz wysyła na twittera (https://twitter.com/tprzechlewski
). Automat działa
na RaspberryPi zresztą...
Dla przypomnienia moja nowa stacja to konsola WiFi o symbolu
DP1500 kupiona u Niemca (czyli w sklepie froggit.de): ze standardowym
zestawem czujników + czujnik pyłu zawieszonego DP200 (też
Froggita). BTW DP1500 to to samo co GW1000 firmy Ecowitt.
Zaś DP 200 Froggita to WH41 Ecowitta
(http://www.ecowitt.com/wifi_weather/83.html
), albo nawet
PM25 firmy AcuWeather
(https://www.ambientweather.com/ampm25.html
)
Zamiast konsoli kupiłem na OLX 5 calowy ekran (110 PLN) i podłączyłem do Raspberry Pi (bezproblemowo).
# konfiguruję (logowanie bez hasła) Boot Options -> Desktop/CLI -> Desktop Autologin/Desktop GUI # # start chromium-browser załadowanie pliku vi ~/.config/lxsession/LXDE-pi/autostart ## wpisuję @/usr/bin/chromium-browser --disable-restore-session-state \ file:///var/www/html/DP1500_live.html
Ekran jest dotykowy. Jak chcę obejrzeć dane pogodowe to dotykam i się wyświetla.
Że jednak brakowało konsoli, to w końcu jednak kupiłem Froggit HP 1000SE PRO (która jest klonem Ecowitt HP 2551). Froggit sprzedaje konsolę jako ersatz, bez niezbędnego specjalnego czujnika temperatury/wilgotności/ciśnienia. Sama stacja jest bowiem dumb -- nic nie mierzy. Ja miałem czujnik temperatury/wilgotności DP50, ale on nie mierzy ciśnienia. Musiałem dokupić oddzielnie ten specjalny, bo się nie dogadałem a dokumentacja/informacja jest w tym aspekcie mało jasna.
Konfigurowanie HP 1000SE PRO to już nie był żaden problem, bo stacja ma aż 8 klawiszy. Połączyłem ją z routerem a potem z dwoma serwisami: ecowitt.net (www.ecowitt.net; oglądanie wymaga zarejestrowania się w sieci ecowitt.net) oraz z WOW ( wow.metoffice.gov.uk)
Kiedyś kupiłem też inny czujnik pyłu zawieszonego (Nova SDS011) i też go podłączyłem do raspberry (bezproblemowo). Teraz przykręciłem go na ścianie obok DP 200 z zamiarem porównania odczytów.
Reasumując posiadam: starą stację WH 2080 działającą od +10 lat.
Nową stację DP1500/WH1000SE PRO oraz czujnik SDS101.
WH 2080 podłączona jest przez kabel USB
do Sheevaplug (to też zaszłość, planuję docelowo przejść na Raspberry)
a obsługiwania jest przez pywws
.
DP1500/WH1000SE podłączona jest przez router, przy czym DP1500 skonfigurowana
jest na serwer lokalny/weewx, a WH1000SE wysyła dane od razu do WOW/ecowitt.net.
Całość kosztował ponad 2500 PLN (wliczając starą stację WH2080). W sumie mogłem kupić HP1000SE PRO + DP200 za 330 EUR, a kupiłem na raty powyższe plus dodatkowy czujnik DP50/DP200 za 450 EUR. Przepłaciłem...