Weblog Tomasza Przechlewskiego

Wględne tempo wzrostu (koronowirusa)

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

url | Tue, 31/03/2020 05:16 | tagi: coronavirus, covid19, r

Dane Eurostatu nt zgonów/urodzeń

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)

url | Wed, 25/03/2020 08:51 | tagi: eurostat, inne, r, statystyka

Dzienne dane dot. wypadków drogowych w Polsce

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)")

url | Wed, 25/03/2020 07:29 | tagi: inne, r, statystyka

TP-Link TL-WN722N v2

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

Kamera

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...

url | Sun, 22/03/2020 04:16 | tagi: raspberry, raspbian, rpi

Dane nt #Covid19 podsumowanie

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).

url | Tue, 17/03/2020 05:12 | tagi: coronavirus, covid19, inne

Dane nt rozwoju epidemii covid19

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ą...

url | Mon, 16/03/2020 07:47 | tagi: coronavirus, covid19

Moja nowa stacja pogody: podsumowanie

HP1000SE_PRO.jpg

Ersatz konsola

SDS011

DP200

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...

url | Sun, 08/03/2020 06:56 | tagi: dp1500, froggit, hp1000se, pogoda, python, raspberrypi, rpi, sds011, weather, weewx