Kaszebe Runda 2020 (pandemiczna)

Mijamy Betlejem

Zdecydowałem się na start na średnim dystansie. Trochę poruta, ale 200 km to dla mnie ciut za dużo zwłaszcza w tym roku. Na stronie jest 1164 zgłoszeń z czego 1010 opłaconych. Do tego: każdy wchodzący do biura zobowiązany jest do poddania się pomiarowi temperatury, zasłonięcia ust i nosa i utrzymania dystansu społecznego w wyznaczonych strefach bezpieczeństwa. Biorąc pod uwagę, że te pomiary i strefy mogą spowodować tłok i bałagan na starcie, zarządzam pobudkę na 4:00 rano.

No więc pojechaliśmy z Elką 4:30 do Kościerzyny żeby na spokojnie wystartować z pierwszą grupą o 7:00 (średni dystans czyli 125km oficjalnie ale chyba w rzeczywistości 5km mniej). Przyjechaliśmy na miejsce jakoś tak 5:30. Generalnie pustki. Z tą temperaturą to lipa była--nikt nie mierzy. Ale za to start kawałek od biura. Zanim się przebrałem, przyczepiłem numer, znalazłem start itd. to w sumie już było 20 minut przed siódmą.

Na starcie tłoku nie ma. Pojechałem jak planowałem z pierwszą grupą--nie wiem czy nas 15 było, no może. Po chwili już było tak na oko z 10. W tym składzie dojechaliśmy do pierwszego bufetu, bo ten wyścig ma taką oto specyfikę że są po drodze bufety gdzie można się pożywić regionalnymi specjałami. Tym razem z tych specjałów to był tylko smalec+ogórek konserwowy+jajecznica. Słabo ale, że nie jadłem śniadania ani nic nie wziąłem do jedzenia, to zjadłem trochę tej jajecznicy + kawałek chleba i ogórek (słony) Większość z naszej grupy na tym bufecie została (konsumując?) a dalej pojechaliśmy w trzech.

Następny przystanek też smalec+ogórek, tyle że kawy nie ma. Pojechaliśmy już dalej we dwóch. Trzeci bufet olaliśmy antycypując że tam też smalec. Zresztą tłum tam był, bo trasa już była wspólna od pewnego czasu z dystansem 55km. BTW dzięki temu zrobiło się trochę mniej nudniej, bo co chwila kogoś wyprzedzaliśmy. Ostatni bufet stajemy bo trzeba coś jednak zjeść. Znowu smalec, więc chwilę tylko tu stoimy i jedziemy dalej. Wreszcie za Betlejem (jest taka wieś na Kaszubach) dogoniliśmy gościa o porównywalnym z naszym tempie. Więc dalej we trzech ciagle mijając tych zmagających się z 55km.

Meta tak z zaskoczenia przyszła, bo patrząc na GPSa myślałem że jeszcze trochę mam do przejechania. A tu za zakrętem jakaś tablica świetlna się pojawiła. To chyba meta jest? W rzeczy samej...

Podium nazwijmy to (gdzie medale dawali) było w jeszcze jednym miejscu, a makaron jeszcze w innym dawali. Że widziałem kilkunastu ludzi, którzy wyglądali na takich co jadą 120km, a 3/4 trasy przejechaliśmy w dwóch/trzech to myślałem, że frekwencja katastrofa a ja tym samym osiągnąłem światowy wynik :-) ale nie -- okazało się, że w kategorii jechało prawie 200, a ja skończyłem na 21 miejscu. Całkiem nieźle anyway...

Po wyścigu i zjedzeniu makaronu wróciłem do domu DK20. Razem wyszło zatem 120+65 = 185km czyli prawie tyle co max-dystans w Kaszebe Runda, ale przyznam się że od Rębiechowa, to ledwo kręciłem pedałami i całe szczęście że od tego miejsca do domu to generalnie z góry jest. Kompletnie uszło ze mnie powietrze. Jednak dystans średni to była racjonalna decyzja...

Po długiej przerwie skorzystałem ze swojej kamery. Wyszło 90Gb filmu, który zamieniłem na 6Gb w niskiej rozdzielczości podzielone na trzy kawałki. Dodałem dane telemetryczne z pliku FIT za pomocą VirbEdit Garmina (próbowałem konkurencyjnego dashware, ale poległem na synchronizacji śladu z wideo; dla mnie to niewykonalne w dashware.)

No więc wyniki są tutaj (skopiowane ze strony organizatora i posortowane). Mój wynik konkretnie to: 27.67km/h i czas 4:20 Ślad GPS tutaj a film z trasy (w trzech częściach) tutaj tutaj oraz tutaj.

Zmarli/zakażeni i zmarli/1mln

W mediach podają albo to co było wczoraj, albo razem od początku świata. Zwłaszcza to łącznie od początku ma niewielką wartość, bo są kraje gdzie pandemia wygasa a są gdzie wręcz przeciwnie. Co z tego, że we Włoszech liczba ofiar to 35 tysięcy jak od miesięcy jest to 5--25 osób dziennie czyli tyle co PL, gdzie zmarło do tej pory kilkanaście razy mniej chorych. Podobnie w Chinach, gdzie cała afera się zaczęła, zmarło we wrześniu 15 osób. Media ponadto lubią duże liczby więc zwykle podają zmarłych, a nie zmarłych na 1mln ludności. Nawet w danych ECDC nie ma tego wskaźnika zresztą. Nie ma też wskaźnika zgony/przypadki, który mocno niefachowo i zgrubie pokazywałby skalę zagrożenia (im większa wartość tym więcej zarażonych umiera)

Taka filozofia stoi za skryptem załączonym niżej. Liczy on zmarłych/1mln ludności oraz współczynnik zgony/przypadki dla 14 ostatnich dni. Wartości przedstawia na wykresach: punktowym, histogramie oraz pudełkowym. W podziale na kontynenty i/lub kraje OECD/pozostałe. Ponieważ krajów na świecie jest ponad 200, to nie pokazuje danych dla krajów w których liczba zmarłych jest mniejsza od 100, co eliminuje małe kraje albo takie, w których liczba ofiar COVID19 też jest mała.

Wynik obok. Skrypt poniżej

#!/usr/bin/env Rscript
#
library("ggplot2")
library("dplyr")
#

## Space before/after \n otherwise \n is ignored
url <- "https://www.ecdc.europa.eu/en/publications-data/ \n download-todays-data-geographic-distribution-covid-19-cases-worldwide"
surl <- sprintf ("retrived %s from %s", tt, url)

## population/continent data:
c <- read.csv("ecdc_countries_names.csv", sep = ';',  header=T, na.string="NA" )

## COVID data
d <- read.csv("covid19_C.csv", sep = ';',  header=T, na.string="NA", 
   colClasses = c('factor', 'factor', 'factor', 'character', 'character', 'numeric', 'numeric'));

# today <- Sys.Date()
# rather newest date in data set:
today <- as.Date(last(d$date))
# 14 days ago
first.day <- today - 14
tt<- format(today, "%Y-%m-%d")
fd<- format(first.day, "%Y-%m-%d")

# select last row for every country
##t <- d %>% group_by(id) %>%  top_n(1, date) %>% as.data.frame
## top_n jest obsolete w nowym dplyr
t <- d %>% group_by(id) %>%  slice_tail(n=1) %>% as.data.frame

# cont. select id/totald columns, rename totald to alld
t <- t %>% select(id, totald) %>% rename(alld = totald) %>% as.data.frame
# join t to d:
d <- left_join(d, t, by = "id")
## extra column alld = number of deaths at `today'
##str(d)

# only countries with > 99 deaths:
d <- d %>% filter (alld > 99 ) %>% as.data.frame

## OECD membership data
o <- read.csv("OECD-list.csv", sep = ';',  header=T, na.string="NA" )

## Remove rows with NA in ID column (just to be sure)
## d <- d %>% filter(!is.na(id))

d$newc <- as.numeric(d$newc)
d$newd <- as.numeric(d$newd)

cat ("### Cleaning newc/newd: assign NA to negatives...\n")
d$newc[ (d$newc < 0) ] <- NA
d$newd[ (d$newd < 0) ] <- NA

## Filter last 14 days only 
d <- d %>% filter(as.Date(date, format="%Y-%m-%d") > fd) %>% as.data.frame
## Last day
last.obs <- last(d$date)
## Format lable firstDay--lastDay (for title):
last14 <- sprintf ("%s--%s", fd, last.obs)

## Sum-up cases and deaths
t <- d %>% group_by(id) %>% summarise( tc = sum(newc, na.rm=TRUE), td = sum(newd, na.rm=TRUE)) %>% as.data.frame 
## Add country populations
t <- left_join(t, c, by = "id")
## Add OECD membership
## if membership column is not NA = OECD member
t <- left_join(t, o, by = "id") 

t$access[ (t$access > 0) ] <- 'oecd'
t$access[ (is.na(t$access)) ] <- 'non-oecd'
str(t)

## Deaths/Cases ratio %%
t$totr <- t$td/t$tc * 100
## Deaths/1mln
t$tot1m <- t$td/t$popData2019 * 1000000

## PL row
dPL <- t[t$id=="PL",]
##str(dPL)
## Extract ratios for PL
totr.PL <- dPL$totr
totr.PL
tot1m.PL <- dPL$tot1m

## Set scales
pScale <- seq(0,max(t$totr, na.rm=T), by=1)
mScale <- seq(0,max(t$tot1m, na.rm=T), by=10)

## Graphs
## 1. deaths/cases ratio (dot-plot)
## color=as.factor(access)): draw OECD/nonOECD with separate colors
p1 <- ggplot(t, aes(x = reorder(country, totr) )) +
  ### One group/one color:
  ##geom_point(aes(y = totr, color='navyblue'), size=1) +
  ### Groups with different colors:
  geom_point(aes(y = totr, color=as.factor(access)), size=1) +
  xlab("country") +
  ylab("%") +
  ggtitle(sprintf("Covid19 death cases ratio (%s)", last14),
    subtitle="Only countries with 100 deaths and more | pale red line = PL level") +
  theme(axis.text = element_text(size = 4)) +
  theme(plot.title = element_text(hjust = 0.5)) +
  ##theme(legend.position="none") +
  scale_color_discrete(name = "Membership") +
  ## Indicate PL-level of  d/c ratio:
  geom_hline(yintercept = totr.PL, color="red", alpha=.25, size=1) +
  labs(caption=surl) +
  scale_y_continuous(breaks=pScale) +
  ##coord_flip(ylim = c(0, 15))
  coord_flip()

ggsave(plot=p1, file="totr_p14.png", width=10)

## 2. deaths/cases ratio (histogram)
p2 <- ggplot(t, aes(x = totr) ) + 
 geom_histogram(binwidth = 0.5, fill='navyblue', alpha=.5) +
 ###
 ylab("N") +
 xlab("%") +
 ggtitle(sprintf("Covid19 deaths/cases ratio (%s)", last14),
   subtitle="Only countries with 100 deaths and more | pale red line = PL level") +
 scale_x_continuous(breaks=pScale) +
 scale_y_continuous(breaks=seq(0, 40, by=2)) +
 geom_vline(xintercept = totr.PL, color="red", alpha=.25, size=1) +
 ##coord_cartesian(ylim = c(0, 30), xlim=c(0, 8))
 labs(caption=surl)

ggsave(plot=p2, file="totr_h14.png", width=10)

## 3. deaths/1m (dot-plot)
## color=as.factor(continent)): draw continents with separate colors
p3 <- ggplot(t, aes(x = reorder(country, tot1m) )) +
  geom_point(aes(y = tot1m, color=as.factor(continent)), size =1 ) +
  xlab("") +
  ylab("deaths") +
  ggtitle(sprintf("Covid19 deaths per 1 million (%s)", last14),
    subtitle="Only countries with 100 deaths and more | red pale line = PL level") +
  theme(axis.text = element_text(size = 6)) +
  geom_hline(yintercept = tot1m.PL, color="red", alpha=.25, size=1) +
  labs(caption=surl) +
  scale_color_discrete(name = "Continent") +
  scale_y_continuous(breaks=mScale) +
  coord_flip()

ggsave(plot=p3, file="totr_m14.png", height=10)

## 3. deaths/1m (box-plots for continents)
p4 <- ggplot(t, aes(x=as.factor(continent), y=tot1m, fill=as.factor(continent))) + 
 geom_boxplot() +
 ylab("deaths") +
 xlab("continent") +
 ggtitle(sprintf("Covid19 deaths per 1 million (%s)", last14),
  subtitle="Only countries with 100 deaths and more") +
 labs(caption=surl) +
 scale_y_continuous(breaks=mScale) +
 theme(axis.text = element_text(size = 6)) +
 theme(legend.position="none")

ggsave(plot=p4, file="totr_c14.png", width=10)

## Mean/Median/IQR per continent
ts <-  t %>% group_by(as.factor(continent)) %>% 
    summarise(Mean=mean(tot1m), Median=median(tot1m), Iqr=IQR(tot1m))
ts

Wartość zmarli/przypadki zależy ewidentnie od słynnej liczby testów. Taki Jemen wygląda, że wcale nie testuje więc ma mało przypadków. Z drugiej strony ci co dużo testują mają dużo przypadków, w tym słynne bezobjawowe więc też nie bardzo wiadomo co to oznacza w praktyce. Bo chyba się testuje po coś, a nie żeby testować dla samego testowania. Więc: dużo testów -- na mój niefachowy rozum -- to wczesne wykrycie, a mało to (w teorii przynajmniej) dużo nieujawnionych/nieizolowanych zarażonych, co powinno dawać w konsekwencji coraz większą/lawinowo większą liczbę przypadków, a co za tym idzie zgonów a tak nie jest.

Konkretyzując rysunek: Ameryka na czele liczby zgonów z medianą prawie 30/1mln, Europa poniżej 3,5/1mln a Azja coś koło 1,4/1mln (średnie odpowiednio 26,4/6,00/5,1.)

Przy okazji się okazało, że mam dplyra w wersji 0.8 co jest najnowszą wersją w moim Debianie. No i w tej wersji 0.8 nie ma funkcji slice_head() oraz slice_tail(), która jest mi teraz potrzebna. Z instalowaniem pakietów R w Linuksie to różnie bywa, ale znalazłem na stronie Faster package installs on Linux with Package Manager beta release coś takiego:

RStudio Package Manager 1.0.12 introduces beta support for the binary format of R packages. In this binary format, the repository's packages are already compiled and can be `installed' much faster -- `installing' the package amounts to just downloading the binary!

In a configured environment, users can access these packages using `install.packages' without any changes to their workflow. For example, if the user is using R 3.5 on Ubuntu 16 (Xenial):

> install.packages("dplyr", repos = "https://demo.rstudiopm.com/cran/__linux__/xenial/latest")

Zadziałało, ale (bo jest zawsze jakieś ale) powtórzenie powyższego na raspberry pi nie zadziała. Wszystko się instaluje tyle że architektury się nie zgadzają. Więc i tak musiałem wykonać

> install.packages("dplyr")
## sprawdzenie wersji pakietów
> sessionInfo()

Instalacja trwała 4 minuty i zakończyła się sukcesem.

Jeszcze jeden skrypt do wizualizacji danych nt COVID19

Jeszcze jeden skrypt do rysowania danych nt COVID19:

#!/usr/bin/env Rscript
# Przekazywanie argumentów z wiersza poleceń
# np.: Rscript --vanilla c19.R --iso PL -clean
library("optparse")
#
library("ggplot2")
library("dplyr")
library("scales")
library("ggpubr")
#
# parametr wygładzania (loess)
spanV <- 0.25
# UWAGA: przed/po \n musi być odstęp inaczej nie działa
surl <- "https://www.ecdc.europa.eu/en/publications-data/ \n download-todays-data-geographic-distribution-covid-19-cases-worldwide"
c0 <- 'PL'

option_list <- list(
  make_option(c("-i", "--iso"), action="store", type="character", default=c0, help="country ISO2 code"),
  make_option(c("-c", "--clean"), action="store_true", default=T, help="extra clean data?")
); 
 
opt_parser <- OptionParser(option_list=option_list);
opt <- parse_args(opt_parser);

c0 <- opt$iso
dataClean <- opt$clean

# wczytanie danych
# date;id;country;newc;newd;totalc;totald
d <- read.csv("covid19_C.csv", sep = ';',  header=T, na.string="NA", 
   colClasses = c('factor', 'factor', 'factor', 'character', 'character', 'numeric', 'numeric'));

str(d)

d$newc <- as.numeric(d$newc)
d$newd <- as.numeric(d$newd)

## Dane zawierają wartości ujemne co jest bez sensu
## z opcją --clean te ujemne wartości są zamieniane na NA
if ( dataClean ) {
  cat ("### Cleaning newc/newd: assign NA to negatives...\n")
  d$newc[ (d$newc < 0) ] <- NA
  d$newd[ (d$newd < 0) ] <- NA
}

## Współczynniki zmarli/zakażeni
d$newr <- d$newd/d$newc * 100
d$totr <- d$totald/d$totalc * 100

## Wartości > 50% są zamieniane na NA (zwykle >50% wynika z błędnych danych)
if ( dataClean ) {
  cat ("### Cleaning newc/newd: assign NA to newr/totr higher than 50...\n")
  d$newr[ (d$newr > 50) ] <- NA
  d$totr[ (d$totr > 50) ] <- NA
}

## Pomiń obserwacje wcześniejsze niż 15/02
d <- d %>% filter(as.Date(date, format="%Y-%m-%d") > "2020-02-15") %>% as.data.frame

d0 <- d %>% filter (id == c0) %>% as.data.frame
t0 <- d0 %>% group_by(id) %>%  summarise(cc = sum(newc, na.rm=T), dd=sum(newd, na.rm=T))

lab0c <- toString(paste (sep=" = ", t0$id, t0$cc))
lab0d <- toString(paste (sep=" = ", t0$id, t0$dd))

## koniecznie dodać na.rm=T bo inaczej zwraca NA (jeżeli znajdzie NA)
maxCC <- max (d0$newc, na.rm=T)
maxDD <- max (d0$newd, na.rm=T)
maxRR <- max (d0$totr, na.rm=T)

last.obs <- last(d0$date)
first.date <- first(d0$date)
fstDay <- as.Date(first.date)
last.totr <- last(d0$totr)
max.newr <- max(d0$newr, na.rm=T)

## Przykład dodania 14 dni do daty
## srcDay <- as.Date(first.date) +14
## https://stackoverflow.com/questions/10322035/r-adding-days-to-a-date
srcDay <- as.Date(last.obs)

## Nazwa pliku wynikowego
## c19_ISO2_DATA.png, gdzie DATA jest datą ostatniej obserwacji
## np.: c19_SE_2020-09-16.png
c0o <- sprintf ("c19_%s_%s.png", c0, last.obs)

## Rysunek1: nowe przypadki
pc0 <- ggplot(d0, aes(x= as.Date(date, format="%Y-%m-%d"), y=newc)) + 
 geom_point(aes(group = id, color = id), size=.8) +
 geom_smooth(method="loess", se=F, span=spanV) +
 theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
 scale_x_date( labels = date_format("%m/%d"), breaks = "2 weeks") +
 annotate("text", x = fstDay, y = 0.95 * maxCC, 
  label = sprintf("Total: %i cases", t0$cc), hjust = 0, vjust=0,
  alpha=0.3, color='steelblue', size=6) +
  xlab(label="") +
 ## Nie drukuj legendy
 theme(legend.position="none") +
 ggtitle(sprintf("%s: new confirmed cases (%s)", t0$id, last.obs))

## Rysunek2: nowe zgony
pd0 <- ggplot(d0, aes(x= as.Date(date, format="%Y-%m-%d"), y=newd)) + 
 geom_point(aes(group = id, color = id), size=.8) +
 geom_smooth(method="loess", se=F, span=spanV) +
 theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
 annotate("text", x = fstDay, y = 0.95 * maxDD, label = sprintf("Total: %i deaths", t0$dd), hjust = 0, vjust=0,
  alpha=0.3, color='steelblue', size=6) +
 scale_x_date( labels = date_format("%m/%d"), breaks = "2 weeks") +
 xlab(label="") +
 theme(legend.position="none") +
 ggtitle(sprintf ("%s: deaths (%s)", t0$id, last.obs))

## Rysunek3: nowe zgony/przypadki *100%
pr0 <- ggplot(d0, aes(x= as.Date(date, format="%Y-%m-%d"), y=newr)) + 
 geom_point(aes(group = id, color = id), size=.8) +
 geom_smooth(method="loess", se=F, span=spanV) +
 theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
 annotate("text", x = fstDay, y = 0.95 * max.newr, label = sprintf("Maximum: %.2f %%", max.newr), hjust = 0, vjust=0,
  alpha=0.3, color='steelblue', size=6) +
 scale_x_date( labels = date_format("%m/%d"), breaks = "2 weeks") +
 xlab(label="") +
 ylab(label="%") +
 theme(legend.position="none") +
 ggtitle(sprintf ("%s: deaths/cases %% (%s)", t0$id, last.obs) )

## Rysunek4: łączne zgony/przypadki *100%
prt0 <- ggplot(d0, aes(x= as.Date(date, format="%Y-%m-%d"), y=totr)) + 
 geom_point(aes(group = id, color = id), size=.8) +
 geom_smooth(method="loess", se=F, span=spanV) +
 theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
 annotate("text", x = fstDay, y = 0.95 * maxRR, label = sprintf("Average: %.2f %%", last.totr), hjust = 0, vjust=0,
  alpha=0.3, color='steelblue', size=6) +
 scale_x_date( labels = date_format("%m/%d"), breaks = "2 weeks") +
 xlab(label="") +
 ylab(label="%") +
 theme(legend.position="none") +
 annotate("text", x = srcDay, y = 0, label = surl, hjust = 1, alpha=.3, size=3) +
 ggtitle(sprintf ("%s: total deaths/cases %% (%s)", t0$id, last.obs) )

p00 <- ggarrange(pc0,pd0, pr0, prt0, ncol=2, nrow=2)

ggsave(plot=p00, c0o, width=15)

Użycie:

Rscript --vanilla c19.R --iso PL

albo:

for i in 'AU' 'BR' 'IN' 'US' 'ES' 'SE' 'PL' 'DE' 'GB'; 
 Rscript --vanilla c19.R --iso $i
done

Przykładowy wynik dla Hiszpanii

Jak widać jakość danych jest katastrofalna: pojawiające się liczne zera to w rzeczywistości brak danych. Zwykle sobota/niedziela zero a potem sruuuu 30 tysięcy bo za trzy dni razem. Wszyscy są zmęczeni w tej Hiszpanii pandemią i nawet nie chce im się danych podsyłać do ECDC?

Rosyjskie dane dotyczące urodzin i zgonów

Reuters donosi, że While Russia has confirmed the world's fourth largest tally of coronavirus cases, it has a relatively low death toll from the associated disease, COVID-19[...] But data released by the Rosstat State Statistics Service on Sept. 4 show there were 57,800 excess deaths between May and July, the peak of the outbreak. The figure was calculated by comparing fatalities over those three months in 2020 with the average number of May-July deaths between 2015 and 2019. The excess total is more than three times greater than the official May-July COVID-19 death toll of 15,955. Co mnie zmotywowało do poszukania danych źródłowych.

Okazało się że dość prosto jest je pobrać ze strony Urzędu Statystycznego Federacji Rosyjskiej (https://eng.gks.ru/ a konkretnie https://showdata.gks.ru/finder/) Ponieważ okazało się to aż tak proste, to pobrałem nie tylko zgony ale także urodzenia. W okresie 2015--2020 (miesięcznie).

Dane są w formacie XSLX. Kolumny to miesiące wiersze poszczególne regiony i ogółem Federacja. Eksportuję do CSV używając LibreOffice i wybieram wiersz z danymi dla całej federacji. W rezultacie mam trzy wiersze: rok-miesiąc, urodzenia i zgony. Obracam moim super skryptem do transpozycji plików CSV:

month;urodzenia;zgony
2015-01-01;149269.99;174722.99
2015-02-01;144591.99;156737
2015-03-01;160974.99;175809.99
...

Trzy wykresy liniowe przedstawiające dynamikę zgonów, urodzin i przyrostu naturalnego (czyli różnicy między urodzinami a zgonami)

library("dplyr")
library("ggplot2")
library("scales")
spanV <- 0.5

## przyrost naturalny
d <- read.csv("UZRT.csv", sep = ';', header=T, na.string="NA");
d$diff <- d$urodzenia - d$zgony

pz <- ggplot(d, aes(x= as.Date(month), y=zgony)) + 
 geom_line(color="steelblue", size=.6, alpha=.5) +
 geom_point(color="steelblue", size=.8) +
 ## Trend 
 geom_smooth(method="loess", se=F, span=spanV, color="red3") +
 ## Skala na osi X-ów
 scale_x_date( labels = date_format("%y/%m"), breaks = "4 months") +
 xlab(label="") +
 ylab(label="deaths") +
  geom_vline(xintercept = as.Date("2015-07-01"), alpha=.25, size=1) +
  geom_vline(xintercept = as.Date("2016-07-01"), alpha=.25, size=1) +
  geom_vline(xintercept = as.Date("2017-07-01"), alpha=.25, size=1)+
  geom_vline(xintercept = as.Date("2018-07-01"), alpha=.25, size=1) +
  geom_vline(xintercept = as.Date("2019-07-01"), alpha=.25, size=1) +
  geom_vline(xintercept = as.Date("2020-07-01"), alpha=.25, size=1) +
 labs(caption="https://showdata.gks.ru/finder/") +
 ggtitle("Russian Federation: Deaths (ths)",  
  subtitle= sprintf("red line: loess with span=%.2f ; gray vertical: july", spanV)) 

pu <- ggplot(d, aes(x= as.Date(month), y=urodzenia)) + 
 geom_line(color="steelblue", size=.6, alpha=.5) +
 geom_point(color="steelblue", size=.8) +
 geom_smooth(method="loess", se=F, span=spanV, color="red3") +
 scale_x_date( labels = date_format("%y/%m"), breaks = "4 months") +
 xlab(label="") +
 ylab(label="births") +
 labs(caption="https://showdata.gks.ru/finder/") +
  geom_vline(xintercept = as.Date("2015-08-01"), alpha=.25, size=1) +
  geom_vline(xintercept = as.Date("2016-08-01"), alpha=.25, size=1) +
  geom_vline(xintercept = as.Date("2017-08-01"), alpha=.25, size=1)+
  geom_vline(xintercept = as.Date("2018-08-01"), alpha=.25, size=1) +
  geom_vline(xintercept = as.Date("2019-08-01"), alpha=.25, size=1) +
  geom_vline(xintercept = as.Date("2020-08-01"), alpha=.25, size=1) +
 ggtitle("Russian Federation: Births (ths)",
  subtitle= sprintf("red line: loess with span=%.2f ; gray vertical: august", spanV)) 

pp <- ggplot(d, aes(x= as.Date(month), y=diff)) + 
 geom_line(color="steelblue", size=.6, alpha=.5) +
 geom_point(color="steelblue", size=.8) +
 geom_smooth(method="loess", se=F, span=spanV, color="red3") +
 scale_x_date( labels = date_format("%y/%m"), breaks = "4 months") +
  geom_vline(xintercept = as.Date("2015-05-01"), alpha=.25, size=1) +
  geom_vline(xintercept = as.Date("2016-05-01"), alpha=.25, size=1) +
  geom_vline(xintercept = as.Date("2017-05-01"), alpha=.25, size=1)+
  geom_vline(xintercept = as.Date("2018-05-01"), alpha=.25, size=1) +
  geom_vline(xintercept = as.Date("2019-05-01"), alpha=.25, size=1) +
  geom_vline(xintercept = as.Date("2020-05-01"), alpha=.25, size=1) +
 xlab(label="") +
 ylab(label="balance") +
 labs(caption="https://showdata.gks.ru/finder/") +
 ggtitle("Russian Federation: natural balance (ths)", 
  subtitle= sprintf("red line: loess with span=%.2f ; gray vertical: may", spanV))

Wykres liniowy dla miesięcy jednoimiennych (MJ):

library(ggplot2)
library(dplyr)
library(tidyr)
library(scales)
options(dplyr.print_max = 1e9)
size0 <- 1.2
size1 <- 0.6
size2 <- 0.6

tt <- read.csv("UZRT.csv", sep = ';',  header=T, na.string="NA");
tt$rok <- substr(tt$month,  1, 4)

tt$yyyymmdd <- as.Date(sprintf ("%s-%s", substr(tt$month,  6, 7),
substr(tt$month,  9, 10) ), format="%m-%d")

## Obliczenie przyrostu naturalnego
tt$diff <- tt$urodzenia - tt$zgony

## Obliczenie 
tt.yy2020 <- tt %>% filter(rok==2020) %>% as.data.frame
tt.yy2019 <- tt %>% filter(rok==2019) %>% as.data.frame
tt.yy2018 <- tt %>% filter(rok==2018) %>% as.data.frame
tt.yy2017 <- tt %>% filter(rok==2017) %>% as.data.frame
tt.yy2016 <- tt %>% filter(rok==2016) %>% as.data.frame
tt.yy2015 <- tt %>% filter(rok==2015) %>% as.data.frame

pf <- ggplot() +
    ggtitle("Russian Federation: births 2015--2020") +
    geom_line( data = tt.yy2020, mapping = aes(x=yyyymmdd, y = urodzenia,  colour = '2020'), alpha=.5, size=size0) +
    geom_point( data = tt.yy2020, mapping = aes(x=yyyymmdd, y = urodzenia,  colour = '2020'), alpha=.5, size=size0) +
    geom_line( data = tt.yy2019, mapping = aes(x=yyyymmdd, y = urodzenia,  colour = '2019'), alpha=.5, size=size1) +
    geom_point( data = tt.yy2019, mapping = aes(x=yyyymmdd, y = urodzenia,  colour = '2019'), alpha=.5, size=size0) +
    geom_line( data = tt.yy2018, mapping = aes(x=yyyymmdd, y = urodzenia,  colour = '2018'), alpha=.5, size=size1) +
    geom_point( data = tt.yy2018, mapping = aes(x=yyyymmdd, y = urodzenia,  colour = '2018'), alpha=.5, size=size1) +
    geom_line( data = tt.yy2017, mapping = aes(x=yyyymmdd, y = urodzenia,  colour = '2017'), alpha=.5, size=size1) +
    geom_point( data = tt.yy2017, mapping = aes(x=yyyymmdd, y = urodzenia,  colour = '2017'), alpha=.5, size=size1) +
    ###
    geom_line( data = tt.yy2016, mapping = aes(x=yyyymmdd, y = urodzenia,  colour = '2016'), alpha=.5, size=size1) +
    geom_point( data = tt.yy2016, mapping = aes(x=yyyymmdd, y = urodzenia,  colour = '2016'), alpha=.5, size=size1) +
    geom_line( data = tt.yy2015, mapping = aes(x=yyyymmdd, y = urodzenia,  colour = '2015'), alpha=.5, size=size2) +
    geom_point( data = tt.yy2015, mapping = aes(x=yyyymmdd, y = urodzenia,  colour = '2015'), alpha=.5, size=size1) +
    scale_x_date( labels = date_format("%y/%m"), breaks = "2 months") +
    labs(caption="https://showdata.gks.ru/finder/") +
    ylab(label="ths") +
    xlab(label="") +
    labs(colour = "") +
    theme(legend.position="top") +
    theme(legend.text=element_text(size=10));

qf <- ggplot() +
    ggtitle("Russian Federation: deaths 2015-2020") +
    geom_line( data = tt.yy2020, mapping = aes(x=yyyymmdd, y = zgony,  colour = '2020'), alpha=.5, size=size0) +
    geom_point( data = tt.yy2020, mapping = aes(x=yyyymmdd, y = zgony,  colour = '2020'), alpha=.5, size=size0) +
    geom_line( data = tt.yy2019, mapping = aes(x=yyyymmdd, y = zgony,  colour = '2019'), alpha=.5, size=size1) +
    geom_point( data = tt.yy2019, mapping = aes(x=yyyymmdd, y = zgony,  colour = '2019'), alpha=.5, size=size0) +
    geom_line( data = tt.yy2018, mapping = aes(x=yyyymmdd, y = zgony,  colour = '2018'), alpha=.5, size=size1) +
    geom_point( data = tt.yy2018, mapping = aes(x=yyyymmdd, y = zgony,  colour = '2018'), alpha=.5, size=size1) +
    geom_line( data = tt.yy2017, mapping = aes(x=yyyymmdd, y = zgony,  colour = '2017'), alpha=.5, size=size1) +
    geom_point( data = tt.yy2017, mapping = aes(x=yyyymmdd, y = zgony,  colour = '2017'), alpha=.5, size=size1) +
    ###
    geom_line( data = tt.yy2016, mapping = aes(x=yyyymmdd, y = zgony,  colour = '2016'), alpha=.5, size=size1) +
    geom_point( data = tt.yy2016, mapping = aes(x=yyyymmdd, y = zgony,  colour = '2016'), alpha=.5, size=size1) +
    geom_line( data = tt.yy2015, mapping = aes(x=yyyymmdd, y = zgony,  colour = '2015'), alpha=.5, size=size2) +
    geom_point( data = tt.yy2015, mapping = aes(x=yyyymmdd, y = zgony,  colour = '2015'), alpha=.5, size=size2) +
    scale_x_date( labels = date_format("%y/%m"), breaks = "2 months") +
    labs(caption="https://showdata.gks.ru/finder/") +
    ylab(label="ths") +
    xlab(label="") +
    labs(colour = "") +
    theme(legend.position="top") +
    theme(legend.text=element_text(size=10));

qf <- ggplot() +
    ggtitle("Russian Federation: natural balance 2015-2020") +
    geom_line( data = tt.yy2020, mapping = aes(x=yyyymmdd, y = diff,  colour = '2020'), alpha=.5, size=size0) +
    geom_point( data = tt.yy2020, mapping = aes(x=yyyymmdd, y = diff,  colour = '2020'), alpha=.5, size=size0) +
    geom_line( data = tt.yy2019, mapping = aes(x=yyyymmdd, y = diff,  colour = '2019'), alpha=.5, size=size1) +
    geom_point( data = tt.yy2019, mapping = aes(x=yyyymmdd, y = diff,  colour = '2019'), alpha=.5, size=size0) +
    geom_line( data = tt.yy2018, mapping = aes(x=yyyymmdd, y = diff,  colour = '2018'), alpha=.5, size=size1) +
    geom_point( data = tt.yy2018, mapping = aes(x=yyyymmdd, y = diff,  colour = '2018'), alpha=.5, size=size1) +
    geom_line( data = tt.yy2017, mapping = aes(x=yyyymmdd, y = diff,  colour = '2017'), alpha=.5, size=size1) +
    geom_point( data = tt.yy2017, mapping = aes(x=yyyymmdd, y = diff,  colour = '2017'), alpha=.5, size=size1) +
    ###
    geom_line( data = tt.yy2016, mapping = aes(x=yyyymmdd, y = diff,  colour = '2016'), alpha=.5, size=size1) +
    geom_point( data = tt.yy2016, mapping = aes(x=yyyymmdd, y = diff,  colour = '2016'), alpha=.5, size=size1) +
    geom_line( data = tt.yy2015, mapping = aes(x=yyyymmdd, y = diff,  colour = '2015'), alpha=.5, size=size2) +
    geom_point( data = tt.yy2015, mapping = aes(x=yyyymmdd, y = diff,  colour = '2015'), alpha=.5, size=size2) +
    scale_x_date( labels = date_format("%y/%m"), breaks = "2 months") +
    labs(caption="https://showdata.gks.ru/finder/") +
    ylab(label="ths") +
    xlab(label="") +
    labs(colour = "") +
    theme(legend.position="top") +
    theme(legend.text=element_text(size=10));

Wreszcie wykres kombinowany liniowo-słupkowy dla miesięcy jednoimiennych (MJ). Linie to 1) średnia liczby zgonów w latach 2015--2020, 2) średnia liczba zgonów w latach 2015--2020 plus/minus odchylenie standardowe liczby zgonów (dla MJ); 3) liczba zgonów w roku 2020. Wykres słupkowy to różnica między liczbą zgonów w roku 2020 a maksymalną liczbą zgonów w latach 2015-2019 ORAZ liczba zgonów z powodu COVID19 (z bazy ECDC). Ten skrypt korzysta z danych z pliku covid_ru.csv, który powstał przez zagregowanie danych dziennych dostępnych ze strony ECDC (dla Federacji Rosyjskiej) oraz policzenie średnich/odchyleń dla jednoimiennych miesięcy z danych w pliku UZRT.csv

library(ggplot2)
library(dplyr)
library(tidyr)
library(scales)
options(dplyr.print_max = 1e9)

#
tt <- read.csv("covid_ru.csv", sep = ';',  header=T, na.string="NA");
tt$diff <- tt$urodzenia - tt$zgony

cols <- c("mean19"="navyblue","mean+sd19"="#3591d1", "mean-sd19"="#3591d1",
  "deaths20"='brown4', "diff19"="red", "c19d"='cyan')

pf <- ggplot(data=tt) +
    ggtitle("Russian Federation: deaths 2015--2020", 
             subtitle='mean19/sd19: mean/sd 2015--19; diff19: 2020 - max(2015-019); deaths20: deaths 2020; c19d: C19 deaths') +
    geom_line(mapping = aes(x=yyyymmdd, y = mean,  colour = 'mean19'), alpha=.5, size=.8) +
    geom_line(mapping = aes(x=yyyymmdd, y = mean + sd,  colour = 'mean+sd19'), alpha=.5, size=.8) +
    geom_line(mapping = aes(x=yyyymmdd, y = mean - sd,  colour = 'mean-sd19'), alpha=.5, size=.8) +
    geom_text(mapping = aes(x=yyyymmdd, label=difflabel, y=diff), vjust=-1.0, size=3 ) +
    ##geom_line(mapping = aes(x=yyyymmdd, y = c19d,  colour = 'c19d'), alpha=.5, size=1.2) +
    labs(caption = "Source: https://showdata.gks.ru/finder/; https://www.ecdc.europa.eu/en/covid-19-pandemic") +
    ##
    geom_bar(mapping = aes(x=yyyymmdd, y = diff, fill = 'diff19' ), color='brown', stat="identity", alpha=.25) +
    geom_bar(mapping = aes(x=yyyymmdd, y = c19d, fill = 'c19d' ), color ='navyblue', stat="identity", alpha=.25) +
    ###
    geom_line(mapping = aes(x=yyyymmdd, y = deaths,  colour = 'deaths20'), alpha=.5, size=1.2) +
    scale_x_date( labels = date_format("%y/%m"), breaks = "2 months") +
    ylab(label="ths") +
    xlab(label="") +
    ##labs(colour = "Legend") +
    scale_colour_manual(name="Lines: ", values=cols) +
    scale_fill_manual(name="Bars: ", values=cols) +
    theme(legend.position="right") +
    theme(legend.text=element_text(size=12));

Skrypt pomocniczy liczący średnie 2015--2019 dla miesięcy jednoimiennych itp rzeczy

library(dplyr)
library(tidyr)
## Agreguje zgodny wg miesięcy jednoimiennych (MJ) dla lat 2015-19
## drukuje średnią/sd/max dla MJ
## Drukuje z2020 - max(z2015--2019)

tt <- read.csv("UZRT.csv", sep = ';',  header=T, na.string="NA");

tt$rok <- substr(tt$month,  1, 4)
tt$yyyymmdd <- as.Date(sprintf ("%s-%s", substr(tt$month,  6, 7), substr(tt$month,  9, 10) ), format="%m-%d")

tt$diff <- tt$urodzenia - tt$zgony

tt.yy2020 <- tt %>% filter(rok==2020) %>% as.data.frame
tt.yy2019 <- tt %>% filter(rok==2019) %>% as.data.frame
tt.yy2018 <- tt %>% filter(rok==2018) %>% as.data.frame
tt.yy2017 <- tt %>% filter(rok==2017) %>% as.data.frame
tt.yy2016 <- tt %>% filter(rok==2016) %>% as.data.frame
tt.yy2015 <- tt %>% filter(rok==2015) %>% as.data.frame

## max zgony
tt.yy2020$max.so.far <-  pmax (tt.yy2019$zgony, tt.yy2018$zgony,
  tt.yy2017$zgony, tt.yy2016$zgony, tt.yy2015$zgony)

## średnia
tt.yy2020$mean <-  (tt.yy2019$zgony + tt.yy2018$zgony
  + tt.yy2017$zgony + tt.yy2016$zgony + tt.yy2015$zgony)/5

## odchylenie std:
tt.yy2020$sqmean <-  (tt.yy2019$zgony^2 + tt.yy2018$zgony^2 +
   tt.yy2017$zgony^2 + tt.yy2016$zgony^2 + tt.yy2015$zgony^2)/5
tt.yy2020$sq <- sqrt(tt.yy2020$sqmean - tt.yy2020$mean^2)

tt.yy2020$diff20 <- tt.yy2020$zgony - tt.yy2020$max.so.far

sprintf ("%s;%.2f;%.2f;%.2f;%.2f", 
    tt.yy2020$yyyymmdd, tt.yy2020$diff20, tt.yy2020$sqmean, tt.yy2020$sq, tt.yy2020$max.so.far )

Ciasto z cukinii

Wariant ciasta z dynią. Wczoraj zrobiłem i wszyscy zachwyceni.

Generalnie to wykonałem przepis ze strony https://beszamel.se.pl/ciasta-z-orzechami/niedzielne-ciasto-z-cukinii,7928/, ale ze zmianami. Oryginał to 3 szklanki startej cukinii + 3 szklanki mąki + 3 jajka + 1 szklanka oleju + 1,5 szklanki cukru + 2 łyżeczki cynamonu + 2 łyżeczki cukru waniliowego + łyżeczka proszku do pieczenia + garść orzeczów włoskich. Wszystko się miesza i piecze w 220C przez 60 min.

Mój zmodyfikowany przepis zachowuje proporcje głównych składników: 2 szklanki cukinii (około 400g) + 2 szklanki mąki i karobu (2 łyżki karobu; reszta mąka). Cukinię ścieram a wodę odlewam więc trzeba znacząco więcej zetrzeć żeby wyszło 2 szklanki, bo cukinia wodnista jest. Dwa jajka skoro po dwie szklanki mąki/cukinii...

Cukru moim zdaniem za dużo jest, podobnie jak oleju. Pół szklanki oleju + trochę więcej niż pół szklanki cukru (ja sypię ksylitol dla zdrowotności). Czubata łyżeczka proszku do pieczenie i taka sama sody. Łyżeczka cynamonu. Orzechy ziemne/włoskie według uznania, czyli garść minimum, ale może być więcej. Dodaję garść rodzynek (ulubionych uzbeckich, takich małych czarno--granatowych.)

Cukinię mieszam z mąką i karobem. Jajka, olej i cukier mieszam intensywnie. Wlewam do cukinii z mąką dodaję proszek, sodę i bakalie. Konsystencja od razu mi wyszła prawidłowa i nie musiałem nic dolewać...

Piekę jak w oryginalnym przepisie 220C przez 60minut.

Wycieczka na Podlasie

Knyszyn

Pojechaliśmy pod koniec sierpnia na krótką wycieczkę na Podlasie. Wszystkiego 5 dni, a konkretnie 24--28 sierpnia czyli poniedziałek--piątek.

Start w poniedziałek około 9:00 i już około 16:00 byliśmy w Białowieży. Tam zakwaterowanie w pokojach pn [Apartamenty] Pod gruszą (ul. Pałacowa 24)... Nie chciało mi się straszenie, ale się zmusiłem i pojechałem rowerem do Hajnówki i z powrotem. Elka z Jankiem poszli zaś na miasto. Wracając zagiąłem jeszcze za miasto i obejrzałem restaurację pn. Carska, która mieści się w dawnym budynku stacyjnym Białowieża Towarowa, jakieś dwa kilometry za centrum.

Następnego dnia rano przed śniadaniem rowerem do Narewki i z powrotem. Po śniadaniu poszliśmy ścieżką pn. Żebra Żubra. Że ścieżka się zaczyna (lub kończy zależy od miejsca gdzie się zaczęło) w rezerwacie pokazowym żubrów obejrzeliśmy i ten rezerwat (średnio warto, a wstęp płatny). Powrót na parking. Potem zwiedzanie skansenu--ciekawy, tyle że względnie niewielki. Przed obiadem obejrzeliśmy jeszcze cerkiew pw. św Mikołaja. Wstęp płatny, a zwiedzanie z przewodnikiem w grupach co 2 godziny. My weszliśmy o 14:00. Interesujące i warto...

Oglądamy jeszcze restaurację Carska (za miasteczkiem, ja już tam byłem wczoraj więc robię za przewodnika), ale bez zamiaru zjedzenia tam czegokolwiek. Miejsce dla snobów. Rekord to bliny z kawiorem za 270 PLN/porcja. Co ciekawe nie ma w Internetach zdjęć jak te super-bliny wyglądają ani opisu, że ktoś zamówił/zjadł i mu smakowało. Nie ma bo ten kto zamówił, żałuje, że wydał 3 stówy bez sensu? (że lokal też nie zamieszcza w tej sytuacji jakby zrozumiałe) Literalnie znalazłem jedną rekomendację i to dotyczącą tańszej wersji tych blinów czyli z wędzonym łososiem. Z kawiorem to nikt się nie chwali, że kupił i mu smakowało...

Obiad w barze pn. Leśna Dziupla. Kartacze + soljanka. Potem do domu a po południu oglądamy park. Plan na jutro: jedziemy do Białegostoku zwiedzając po drodze Kruszyniany, Krynki i Supraśl. Kupujemy sera carskiego w liczbie siedem sztuk (w tym trzy dla znajomej Elki w Gdańsku).

Rano wstaję i jadę rowerem do Kruszynian (przez Narewkę). Elka z Jankiem też tam jadą, tyle że samochodem. Spotykamy się po 11:00. Idziemy zwiedzać meczet i cmentarz. Kupa luda jest tutaj, parking pełen samochodów. Po zwiedzaniu idziemy na kawę do Tatarskiej Jurty, ale tam generalnie obiadowe menu, a że na obiad za wcześnie nic nie zamawiamy i jedziemy do Krynek.

Krynki dziś to dziura, ale kiedyś to było duże miasteczko w którym mieszkali w większości Żydzi. Podobno 13 tysięcy przed pierwszą wojną, potem już tylko pięć, teraz zaledwie dwa no i nie Żydów. Centrum przemysłu garbarskiego (kiedyś)... Z tych czasów świetności pozostały dwa budynki po synagogach i cmentarz żydowski. Oglądamy synagogę przy Kaukaskiej -- teraz gminny ośrodek kultury. Jest tam nawet coś w rodzaju izby pamięci, ale mocno skromne: kilkanaście zdjęć na ścianie + jakieś świeczniki itp... Za to cmentarz duży -- kilkaset nagrobków.

Następny przystanek Supraśl, do którego docieramy około 14:00. Idziemy do klasztoru a tu przerwa obiadowa. No to my też na obiad. Potem wracamy i zwiedzamy kompleks klasztorny (ciekawe). Główna cerkiew została odbudowana całkiem niedawno bo w czasie wojny została zburzona, i jeszcze nie jest poświęcona. Zresztą w czasie zwiedzania widzimy gościa, który maluje freski na ścianie... Wychodzimy i Elka spotyka koleżankę z pracy, która czeka na wejście z następną grupą. W efekcie i w euforii zapominamy, że mieliśmy jeszcze w planach zwiedzanie muzeum ikon...

Około 17:00 jesteśmy w Białymstoku u koleżeństwa KBS. Idziemy wieczorem na miasto. Jutro będziemy zwiedzać okolicę czyli głównie Tykocin...

W Tykocinie zaczynamy od synagogi potem idziemy na zamek. Zamek wygląda całkiem, całkiem, ale się okazuje, że to nie jest oryginał--odbudowany. W środku nie ma rewelacji ale p. przewodniczka opowiada historię miejsca naprawdę ciekawie i zajmująco, więc warto kupić bilet. Potem obiad w restauracji przy moście nad Narwią. Po obiedzie, że jeszcze wcześnie i głupio wracać do Białegostoku, jedziemy do Knyszyna.

Dojeżdżamy, ale co my teraz mamy zwiedzać? Miał być zamek tyle że go nie ma... Był ale już nie ma. Koleżeństwo z Elką poszło studiować wielki plan miasta w rynku w poszukiwaniu atrakcji, a ja wbijam do smartfona: cmentarz+żydowski. Jest... To jedziemy...

No i warto było. Podobno największy na Podlasiu. Zwiedzamy zatem dokładnie...

Piątek. Wracamy do domu, ale przed powrotem jedziemy do Supraśla raz jeszcze zwiedzić muzeum ikon. Bardzo ciekawe miejsce...

W domu jesteśmy około 19:00

Do pobrania/oglądania ślady kml ze zdjęciami; album ze zdjęciami.

Wydłubywanie danych nt/ COVID19 z tweetów MZ

MZ to Ministerstwo Zdrowia. Specyfiką PL jest brak danych publicznych nt. Pandemii.. Na stronie GIS nie ma nic, a stacje wojewódzkie publikują jak chcą i co chcą. Ministerstwo zdrowia z kolei na swojej stronie podaje tylko dane na bieżący dzień, zaś na amerykańskim Twitterze publikuje komunikaty na temat. To co jest na stronie znika oczywiście następnego dnia zastępowane przez inne cyferki. Do tego są tam tylko dzienne dane zbiorcze o liczbie zarażonych i zmarłych, (w podziale na województwa). Nie ma na przykład wieku, płci itp... To co jest na Twitterze z kolei ma formę tekstowego komunikatu postaci: Mamy 502 nowe i potwierdzone przypadki zakażenia #koronawirus z województw: małopolskiego (79), śląskiego (77), mazowieckiego (75)... [...] Z przykrością informujemy o śmierci 6 osób zakażonych koronawirusem (wiek-płeć, miejsce zgonu): 73-K Kędzierzyn-Koźle, 75-M Łańcut, 92-K Lipie, 72-M, 87-M, 85-K Gdańsk.

Czyli podają zarażonych ogółem i w podziale na województwa oraz dane indywidualne zmarłych w postaci płci i wieku oraz miejsca zgonu (miasta żeby było inaczej niż w przypadku podawania zakażeń.)

No więc chcę wydłubać dane postaci 92-K z tweetów publikowanych przez Ministerstwo Zdrowia. W tym celu za pomocą tweepy pobieram cały streamline (+3200 tweetów zaczynających się jeszcze w 2019 roku), ale dalej to już działam w Perlu, bo w Pythonie jakby mniej komfortowo się czuję. Po pierwsze zamieniam format json na csv:

use JSON;
use Data::Dumper;
use Time::Piece;
use open ":encoding(utf8)";
use open IN => ":encoding(utf8)", OUT => ":utf8";
binmode(STDOUT, ":utf8");

##  ID = tweeta ID; date = data; 
##  repid -- odpowiedź na tweeta o numerze ID
##  text -- tekst tweeta
print "id;date;repid;text\n";

while (<>) {
  chomp();
  $tweet =  $_;

  my $json = decode_json( $tweet );
  #print Dumper($json);
  $tid = $json->{"id"};
  $dat = $json->{"created_at"};
  ## Data jest w formacie rozwlekłym zamieniamy na YYYY-MM-DDTHH:MM:SS
  ## Fri Oct 04 14:48:25 +0000 2019
  $dat = Time::Piece->strptime($dat,
     "%a %b %d %H:%M:%S %z %Y")->strftime("%Y-%m-%dT%H:%M:%S");
  $rep = $json->{"in_reply_to_status_id"};
  $ttx = $json->{"full_text"}; $ttx =~ s/\n/ /g;
  ## Zamieniamy ; na , w tekście żeby użyć ; jako separatora
  $ttx =~ s/;/,/g; ####

  print "$tid;$dat;$rep;$ttx\n";

Komunikaty dłuższe niż limit Twittera są dzielone na kawałki, z których każdy jest odpowiedzią na poprzedni, np:

1298900644267544576;2020-08-27T08:30:20;1298900642522685440;53-M, 78-K i 84-K Kraków. Większość osób ...
1298900642522685440;2020-08-27T08:30:20;1298900640714948608;67-K Lublin (mieszkanka woj. podkarpackiego), 85-K Łańcut,...
1298900640714948608;2020-08-27T08:30:19;1298900639586680833;kujawsko-pomorskiego (24), świętokrzyskiego (18), opolskiego...
1298900639586680833;2020-08-27T08:30:19;;Mamy 887 nowych i potwierdzonych przypadków zakażenia #koronawirus z województw: ...

Czyli tweet 1298900639586680833 zaczyna, 1298900640714948608 jest odpowiedzią na 1298900639586680833, a 1298900642522685440 odpowiedzią na 1298900640714948608 itd. Tak to sobie wymyślili... W sumie chyba niepotrzebnie, ale w pierwszym kroku agreguję podzielone komunikaty w ten sposób, że wszystkie odpowiedzi są dołączane do pierwszego tweeta (tego z pustym polem in_reply_to_status_id):

## nextRef jest rekurencyjna zwraca numer-tweeta,
## który jest początkiem wątku
sub nextRef {
  my $i = shift;

  if ( $RR{"$i"} > 0 ) {  
    return ( nextRef( "$RR{$i}" ) );
  } else { return "$i" }
}

### ### ###
while (<>) { chomp();
   ($id, $d, $r, $t) = split /;/, $_;

   $TT{$id} = $t;
   $RR{$id} = $r;
   $DD{$id} = $d;
}

### ### ###
for $id ( sort keys %TT ) {  
   $lastId = nextRef("$id");
   $LL{"$id"} = $lastId;
   $LLIds{"$lastId"} = "$lastId"; 
}

### ### ###
for $id (sort keys %TT) {  
    ## print "### $DD{$id};$id;$LL{$id};$TT{$id}\n";  }
    $TTX{$LL{"$id"}} .= " " . $TT{"$id"};
    $DDX{$LL{"$id"}} .= ";" . $DD{"$id"};
}
### ### ###

for $i (sort keys %TTX) { 

    $dates = $DDX{"$i"};
    $dates =~ s/^;//; ## pierwszy ; jest nadmiarowy
    @tmpDat = split /;/, $dates;

    $dat_time_ = $tmpDat[0]; 
    ($dat_, $time_) = split /T/, $dat_time_;
    $ffN = $#tmpDat + 1;
    $collapsedTweet = $TTX{$i};
    print "$i;$dat_;$time_;$ffN;$collapsedTweet\n";
}

Zapuszczenie powyższego powoduje konsolidację wątków, tj. np. powyższe 4 tweety z 2020-08-27 połączyły się w jeden:

1298900639586680833;2020-08-27;08:30:19;4; Mamy 887 nowych i potwierdzonych przypadków
zakażenia #koronawirus z województw: małopolskiego (233), śląskiego (118), mazowieckiego (107),
[...] Liczba zakażonych koronawirusem: 64 689 /2 010 (wszystkie pozytywne przypadki/w tym osoby zmarłe).

Teraz wydłubujemy tylko tweety z frazą nowych i potwierdzonych albo nowe i potwierdzone:

## nowych i potwierdzonych albo nowe i potwierdzone
## (MZ_09.csv to CSV ze `skonsolidowanymi' wątkami)
cat MZ_09.csv | grep 'nowych i pot\|nowe i pot' > MZ_09_C19.csv
wc -l MZ_09_C19.csv
189 MZ_09_C19.csv

Wydłubanie fraz wiek-płeć ze skonsolidowanych wątków jest teraz proste:

  perl -e '
while (<>) {
 ($i, $d, $t, $c, $t) =  split /;/, $_;

  while ($t =~ m/([0-9]+-[MK])/g ) {
   ($w, $p) = split /\-/, $1;
   print "$d;$w;$p\n";
  }
}' MZ_09_C19.csv > C19D.csv
wc -l C19PL_down.csv
1738

Plik wykazał 1738 osób. Pierwsza komunikat jest z 16 kwietnia. Ostatni z 31. sierpnia. Pierwsze zarejestrowane zgony w PL odnotowano 12 marca (albo 13 nieważne). W okresie 12 marca --15 kwietnia zmarło 286 osób. Dodając 286 do 1738 wychodzi 2024. Wg MZ w okresie 12.03--31.08 zmarło 2039. Czyli manko wielkości 15 zgonów (około 0,5%). No trudno, nie chce mi się dociekać kto i kiedy pogubił tych 15...

Równie prostym skryptem zamieniam dane indywidualne na tygodniowe

#!/usr/bin/perl -w
use Date::Calc qw(Week_Number);

while (<>) {
  chomp();
  ##if ($_ =~ /age/) { next }  ##

  my ($d, $w, $p ) = split /;/, $_;
  my ($y_, $m_, $d_) = split /\-/, $d;
  my $week = Week_Number($y_, $m_, $d_);

  $DW{$week} += $w;
  $DN{$week}++;

  $DD{$d} = $week;
  $YY{$week} = 0;
  $YY_last_day{$week} = $d;

  ## wiek wg płci
  $PW{$p} += $w;
  $PN{$p}++;
}

for $d (keys %DD) {
    $YY{"$DD{$d}"}++; ## ile dni w tygodniu   
 }

print STDERR "Wg płci/wieku (ogółem)\n";

for $p (keys %PW) {
  $s = $PW{$p}/$PN{$p};
  printf STDERR "%s %.2f %i\n", $p, $s, $PN{$p};
  $total += $PN{$p};
}

print STDERR "Razem: $total\n";
print "week;deaths;age;days;date\n";

for $d (sort keys %DW) {
  if ($YY{$d} > 2 ) {## co najmniej 3 dni 
    $s = $DW{$d}/$DN{$d};
    printf "%s;%i;%.2f;%i;%s\n", $d, $DN{$d}, $s, $YY{$d}, $YY_last_day{$d};
  }
}

Co daje ostatecznie (week -- numer tygodnia w roku; meanage -- średni wiek zmarłych; deaths -- liczba zmarłych w tygodniu; days -- dni w tygodniu; date -- ostatni dzień tygodnia):

week;deaths;meanage;days;date
16;55;77.07;4;2020-04-19
17;172;75.09;7;2020-04-26
18;144;77.29;7;2020-05-03
19;123;76.46;7;2020-05-10
20;126;76.40;7;2020-05-17
21;71;76.37;7;2020-05-24
22;68;78.12;7;2020-05-31
23;93;75.73;7;2020-06-07
24;91;75.93;7;2020-06-14
25;109;77.24;7;2020-06-21
26;83;75.06;7;2020-06-28
27;77;74.09;7;2020-07-05
28;55;76.91;7;2020-07-12
29;54;77.33;7;2020-07-19
30;48;76.52;7;2020-07-26
31;60;74.88;7;2020-08-02
32;76;77.17;7;2020-08-09
33;71;73.11;7;2020-08-16
34;77;75.61;7;2020-08-23
35;79;74.33;7;2020-08-30

To już można na wykresie przedstawić:-)

library("dplyr")
library("ggplot2")
library("scales")
##
spanV <- 0.25
d <- read.csv("C19D_weekly.csv", sep = ';',  header=T, na.string="NA")

first <- first(d$date)
last <- last(d$date)
period <- sprintf ("%s--%s", first, last)
d$deaths.dailymean <- d$deaths/d$days

cases <- sum(d$deaths);
max.cases <- max(d$deaths)

note <- sprintf ("N: %i (source: twitter.com/MZ_GOV_PL)", cases)

pf <- ggplot(d, aes(x= as.Date(date), y=meanage)) +
 geom_bar(position="dodge", stat="identity", fill="steelblue") +
 scale_x_date( labels = date_format("%m/%d"), breaks = "2 weeks") +
 scale_y_continuous(breaks=c(0,5,10,15,20,25,30,35,40,45,50,55,60,65,70,75,80)) +
 xlab(label="week") +
 ylab(label="mean age") +
 ggtitle(sprintf ("Mean age of COVID19 fatalities/Poland/%s", period), subtitle=note ) 

note <- sprintf ("N: %i (source: twitter.com/MZ_GOV_PL)", cases)
pg <- ggplot(d, aes(x= as.Date(date), y=deaths.dailymean)) +
 geom_bar(position="dodge", stat="identity", fill="steelblue") +
 scale_x_date( labels = date_format("%m/%d"), breaks = "2 weeks") +
 scale_y_continuous(breaks=c(0,5,10,15,20,25,30,35,40,45,50)) +
 xlab(label="week") +
 ylab(label="daily mean") +
 ggtitle(sprintf("Daily mean number of COVID19 fatalities/Poland/%s", period), subtitle=note )

ggsave(plot=pf, file="c19dMA_w.png")
ggsave(plot=pg, file="c19dN_w.png")

Wynik tutaj:

Średnia Wg płci/wieku/ogółem. Ogółem -- 75,9 lat (1738 zgonów) w tym: mężczyźni (914 zgonów) -- 73.9 lat oraz kobiety (824) -- 78.4 lat. Stan pandemii na 31.08 przypominam. Apogeum 2 fali... Średnia wieku w PL (2019 rok) 74,1 lat (M) oraz 81,8 lat (K).