Jednego dnia zmarli Maradona, Christophe Dominici oraz Zenon Plech. A przynajmniej jednego dnia media o tym poinformowały. Plech to lokalna polska gwiazda żużla. Dominici gwiazda francuskiego rugby, uczestnik słynnego meczu z Nową Zelandią w półfinale pucharu świata w 1999r. Maradona to wiadomo i (prawie) każdy czuł się w obowiązku poinformować o tym wydarzeniu swoją publiczność, nawet tacy co mają po 10 folowersów na Twitterze.

Dla mnie Maradona to oszust przede wszystkim i prymityw (w przeciwieństwie do Pelego na przykład). Każdy piłkarz pewnie ma na sumieniu faul w ważnym meczu, typu bramka ze spalonego, ale nie każdy popełnił faul aż tak oczywisty i bezczelny. Że to było dawno to większość nie kojarzy, ale wszyscy na stadionie będący dostatecznie blisko widzieli rękę, bo jak nie zobaczyć czegoś tak ordynarnego za wyjątkiem sędziego głównego + liniowego. Tak się dziwnie złożyło...

W piłce kopanej na turniejach o Puchar Świata jest (była?) w ogóle taka IMO patologia, że sędziują ludzie często słabi, bo z całego świata. Na co dzień sędziuje taki as ligę w Gwatemali dajmy na to, a raz w życiu będzie sędziował Niemcy vs Meksyk. Ten konkretny mecz sędziowali Tunezyjczyk + Bułgar (liniowy). Inne ćwierćfinały sędziowali wtedy: DDR-owiec, Rumun i Kolumbijczyk. No same mocne ligi na co dzień...

El mano to po hiszpańsku ręka jest właśnie. Na zdjęciu rzekomo kontrowersyjna sytuacja. Ona by była kontrowersyjna gdyby wepchnął piłkę w tłoku, albo byłaby wątpliwość czy ręką czy barkiem. Tak bywa... Ale tutaj? Dwóch ludzi skacze do piłki, i wszyscy widzę że ten co ma 160 w kapeluszu przeskakuje bramkarza. No ludzie...

Szybkie bo za pomocą google apps script (GAS).

Nie miałem do wczoraj świadomości że można automatyzować różne rzeczy korzystając z narzędzia pn Google Apps Script. Konkretnie zaczynamy tu. A nawet konkretniej https://script.google.com/home/my. Oczywiście trzeba być zalogowanym żeby link zadziałał.

No więc kliknąłem nowy projekt a potem zmodyfikowałem przykład z samouczka bo zawierał wszystko co potrzeba:

// Create a new form, then add a checkbox question, a multiple choice question,
// a page break, then a date question and a grid of questions.
function convertToForm(){
var form = FormApp.create('New Form')
form
.setTitle('Badanie wizerunku i marki)
.setDescription('Badamy wizerunek i markę')

 form.addPageBreakItem()
    .setTitle('Część 1/3');               
                
form.addCheckboxItem()
   .setTitle('What condiments would you like on your hot dog?')
   .setChoiceValues(['Ketchup', 'Mustard', 'Relish'])

form.addMultipleChoiceItem()
    .setTitle('Do you prefer cats or dogs?')
    .setChoiceValues(['Cats','Dogs'])
    .showOtherOption(true);
  
form.addPageBreakItem()
    .setTitle('Część 2/3');

form.addDateItem()
    .setTitle('When were you born?');

form.addGridItem()
    .setTitle('Rate your interests')
    .setRows(['Cars', 'Computers', 'Celebrities'])
    .setColumns(['Boring', 'So-so', 'Interesting']);
  
form.addScaleItem()
     .setTitle('Scale item')
     .setLabels('zły', 'dobry')

form.addTextItem()
     .setTitle('Uwagi')
     .setRequired(true)

form.addMultipleChoiceItem()
    .setTitle('Płeć')
    .setChoiceValues(['K','M'])
    .setRequired(true);

form.addTextItem()
     .setTitle('Wiek')
     .setRequired(true)

Logger.log('Published URL: ' + form.getPublishedUrl());
Logger.log('Editor URL: ' + form.getEditUrl());
}

Następnie Save. Uruchamiamy przez Run/Run function. Jeżeli nie będzie błędów pod adresem https://docs.google.com/forms/u/0/ pojawi się formularz. Jak coś nie pasuje, to można resztę doklikać. Długie formularze w ułamku tego czasu który jest potrzebny do ręcznego wklikania da się zrobić...

Jak wiadomo w wyniku działań pewnego upierdliwego 19 latka okazało się, że sprawozdawczość MinZdrowia/GiS to farsa. Zadeklarowano w urzędzie wzw/ z tym nowe otwarcie: stary system zamknięto (he, he) nowy wystartował 23.11.2020. Przy okazji dodano do puli zakażonych 22594 przypadków wcześniej nie uwzględnionych (około 2,5%). Zatem oficjalnie 23-11-2020 było 10139 zakażeń na stronie MZ; ale wszyscy inni (ECDC na przykład) wpisują 32733 (czyli 10139 + 22594) no bo jak inaczej. W sumie to nawet ma sens, bo przecież baza odpowiada zakażeniem zarejestrowanym a nie, że tak powiem nabytym. Skoro je zarejestrowano 23-11-2020 to tyle ma być...

Pozostaje problem jak to rozliczyć na województwa. No więc ja chciałem rozliczyć tak, że do zakażeń z 23-11-2020 dodam te 22594 w proporcji do wielkości zakażeń z 23-11-2020. Ściągnąłem nowe-dobre dane ze strony wykaz-zarazen-koronawirusem i zacząłem od zsumowania zakażeń w województwach. Się okazało, że jest to inna liczba niż 10139. Kuźwa. W pliku z 24-11-2020 to samo...

No więc miało być lepiej a wyszło tak samo. Uprawnione jest w tym momencie pytanie czy ludzie nie potrafiący liczyć muszą akurat walczyć z pandemią i pracować w GiS/MZ? Załączam rysunki, żeby nie było że zmyślam. Liczenie zakażeń/zgonów wg województw na razie zawieszam skoro nie bardzo wiadomo co liczyć...

GUS się wychylił niespodziewanie z dużą paczką danych nt. zgonów. Dane są tygodniowe, w podziale na płcie, regiony w klasyfikacji NUTS oraz 5 letnie grupy wiekowe.

Dane są udostępnione w formacie XSLX w dość niepraktycznej z punktu widzenia przetwarzania strukturze (kolumny to tygodnie/wiersze to różne kategorie agregacji: płeć, wiek, region), który zamieniłem na CSV o następującej prostej 7 kolumnowej strukturze:

year;sex;week;date;age;geo;value

W miarę oczywiste jest, że year to rok, sex to płeć, week to numer tygodnia, date to data pierwszego dnia tygodnia (poniedziałek), geo to identyfikator obszaru a value liczba zgonów odpowiadająca kolumn 1--6. Ten plik jest podzielony na lata bo w całości zajmuje circa 200Mb. Umieściłem go tutaj.

Skrypt też w R wymodziłem co wizualizuje zgony wg grup wieku oraz województw. Ponieważ kombinacji płeć/wiek/region są setki, moje wykresy dotyczą zgonów ogółem/kobiet/mężczyzn w podziale na grupy wiekowe oraz ogółem w podziale na województwa. Każdy wykres zawiera dwa szeregi: liczbę zgonów w 2020 roku oraz średnią liczbę zgonów z lat 2015--2019. Ponadto jest wykres z jedną krzywą: procent liczony dla stosownych tygodni jako liczba zgonów w 2020 przez średnią 5 letnią z lat 2015--2019. Ten wykres występuje też w wariancie skróconym: tylko 6 ostatnich tygodni, co pozwala dodać do punktów wartości liczbowe (które nie zachodzą na siebie).

library("ggplot2")
library("dplyr")
library("scales")
library("ggthemes")
library("ggpubr")
library("tidyr")

picWd <- 12
spanV <- 0.5
GUS.url <- "https://stat.gov.pl/obszary-tematyczne/ludnosc/ludnosc/zgony-wedlug-tygodni,39,2.html"
NIKW.url <- "(c) NI-KW @ github.com/knsm-psw/GUS_mortality"
NIKW <- sprintf ("%s | %s", GUS, NIKW.url)

z <- read.csv("PL-mortality-2015.csv", sep = ';',  header=T, na.string="NA" )
lastO <- last(z$date)
lastT <- last(z$week)

nuts <- c('PL21', 'PL22', 'PL41', 'PL42', 'PL43', 'PL51', 'PL52', 'PL61', 'PL62',
'PL63', 'PL71', 'PL72', 'PL81', 'PL82', 'PL84', 'PL91', 'PL92')

### Ogółem
z00 <- z %>% filter ( sex == 'O'  & geo == 'PL' ) %>% as.data.frame

z0 <- z00 %>% filter ( year >= 2015  & year < 2020 ) %>% as.data.frame
z1 <- z00 %>% filter ( year == 2020 ) %>% as.data.frame

## średnie w okresie 1 -- (n-1)
zz0 <- z0 %>% group_by(age,week) %>% summarise( year = 't19',
  vv = mean(value, na.rm=TRUE)) %>% as.data.frame
zz1 <- z1 %>% group_by(age,week) %>% summarise( year = 't20',
  vv = mean(value, na.rm=TRUE)) %>% as.data.frame
### Połącz
zz1 <- bind_rows(zz0, zz1)

farbe19 <- '#F8766D'
farbe20 <- '#00BFC4'

p1 <- ggplot(zz1, aes(x=week, y=vv, color=year)) +
 geom_smooth(method="loess", se=F, span=spanV, size=.4) +
 geom_point(size=.4, alpha=.5) +
 facet_wrap( ~age, scales = "free_y") +
 xlab(label="") +
 ylab(label="") +
 ##theme_nikw()+
 theme(plot.subtitle=element_text(size=9), legend.position="top")+
 scale_color_manual(name="Rok: ", labels = c("średnia 2015--2019", "2020"), values = c("t19"=farbe19, "t20"=farbe20 )) +
 ggtitle("Zgony wg grup wiekowych (PL/Ogółem)", subtitle=sprintf("%s | ostatni tydzień: %s", NIKW, lastO) )

ggsave(plot=p1, "zgony_PL_by_age_O.png", width=picWd)

### M ###
z00 <- z %>% filter ( sex == 'M'  & geo == 'PL' ) %>% as.data.frame

z0 <- z00 %>% filter ( year >= 2015  & year < 2020 ) %>% as.data.frame
z1 <- z00 %>% filter ( year == 2020 ) %>% as.data.frame

## średnie w okresie 1 -- (n-1)
zz0 <- z0 %>% group_by(age,week) %>% summarise( year = 't19',
  vv = mean(value, na.rm=TRUE)) %>% as.data.frame
zz1 <- z1 %>% group_by(age,week) %>% summarise( year = 't20',
  vv = mean(value, na.rm=TRUE)) %>% as.data.frame
### Połącz
zz1 <- bind_rows(zz0, zz1)

p2 <- ggplot(zz1, aes(x=week, y=vv, group=year, color=year)) +
 geom_smooth(method="loess", se=F, span=spanV, size=.4) +
 geom_point(size=.4, alpha=.5) +
 facet_wrap( ~age, scales = "free_y") +
 xlab(label="") +
 ylab(label="") +
 ##theme_nikw()+
 ##labs(caption=source) +
 theme(plot.subtitle=element_text(size=9), legend.position="top")+
 scale_color_manual(name="Rok: ", labels = c("średnia 2015--2019", "2020"),
   values = c("t19"=farbe19, "t20"=farbe20 )) +
 ggtitle("Zgony wg grup wiekowych (PL/Mężczyźni)", subtitle=sprintf("%s | ostatni tydzień: %s", NIKW, lastO) )

ggsave(plot=p2, "zgony_PL_by_age_M.png", width=picWd)


### K #########################################
z00 <- z %>% filter ( sex == 'K'  & geo == 'PL' ) %>% as.data.frame

z0 <- z00 %>% filter ( year >= 2015  & year < 2020 ) %>% as.data.frame
z1 <- z00 %>% filter ( year == 2020 ) %>% as.data.frame

## średnie w okresie 1 -- (n-1)
zz0 <- z0 %>% group_by(age,week) %>% summarise( year = 't19',
  vv = mean(value, na.rm=TRUE)) %>% as.data.frame
zz1 <- z1 %>% group_by(age,week) %>% summarise( year = 't20',
  vv = mean(value, na.rm=TRUE)) %>% as.data.frame
### Połącz
zz1 <- bind_rows(zz0, zz1)

p3 <- ggplot(zz1, aes(x=week, y=vv, group=year, color=year)) +
 geom_smooth(method="loess", se=F, span=spanV, size=.4) +
 geom_point(size=.4, alpha=.5) +
 facet_wrap( ~age, scales = "free_y") +
 xlab(label="") +
 ylab(label="") +
 ##theme_nikw()+
 ##labs(caption=source) +
 theme(plot.subtitle=element_text(size=9), legend.position="top")+
 scale_color_manual(name="Rok: ", labels = c("średnia 2015--2019", "2020"),
   values = c("t19"=farbe19, "t20"=farbe20 )) +
 ggtitle("Zgony wg grup wiekowych (PL/Kobiety)", subtitle=sprintf("%s | ostatni tydzień: %s", NIKW, lastO) )

ggsave(plot=p3, "zgony_PL_by_age_K.png", width= picWd)

### ogółem wg województw #####################################
n <- read.csv("nuts.csv", sep = ';',  header=T, na.string="NA" )
## dodaj nazwy
z <- left_join(z, n, by='geo')

## wiek razem
z00 <- z %>% filter ( sex == 'O' & geo %in% nuts & age == 'OGÓŁEM') %>% as.data.frame

z0 <- z00 %>% filter ( year >= 2015  & year < 2020 ) %>% as.data.frame
z1 <- z00 %>% filter ( year == 2020 ) %>% as.data.frame

## średnie w okresie 1 -- (n-1)
zz0 <- z0 %>% group_by(name,week) %>%
 summarise( year = 't19', vv = mean(value, na.rm=TRUE)) %>% as.data.frame
 zz1 <- z1 %>% group_by(name,week) %>%
  summarise( year = 't20', vv = mean(value, na.rm=TRUE)) %>% as.data.frame
### Połącz
zz1 <- bind_rows(zz0, zz1)

lastWeek <- last(zz1$week)
firstWeek <- lastWeek - 6

zz1 <- zz1 %>% filter ( week >= firstWeek  ) %>% as.data.frame
print(zz1)

p4 <- ggplot(zz1, aes(x=week, y=vv, group=year, color=year)) +
 geom_smooth(method="loess", se=F, span=spanV, size=.4) +
 geom_point(size=.4, alpha=.5) +
 facet_wrap( ~name, scales = "free_y") +
 xlab(label="") +
 ylab(label="") +
 ##theme_nikw()+
 ##labs(caption=source) +
 theme(plot.subtitle=element_text(size=9), legend.position="top")+
 scale_color_manual(name="Rok: ", labels = c("średnia 2015--2019", "2020"),
   values = c("t19"=farbe19, "t20"=farbe20 )) +
 ggtitle("Zgony wg województw* (PL/Ogółem)", 
   subtitle=sprintf("*wg klasyfikacji NUTS stąd mazowieckie/stołeczne | %s | ostatni tydzień: %s", NIKW, lastO))

ggsave(plot=p4, "zgony_PL_by_woj_O.png", width=picWd)

## jako %% w średniej w poprzednich 5 lat

zz1 <- zz1 %>% spread(year, vv)

zz1$yy <- zz1$t20 / zz1$t19 * 100

p5 <- ggplot(zz1, aes(x=week, y=yy), color=farbe20) +
 geom_smooth(method="loess", se=F, span=spanV, size=.4, color=farbe20) +
 geom_point(size=.4, alpha=.5) +
 facet_wrap( ~name, scales = "fixed") +
 xlab(label="nr tygodnia") +
 ylab(label="%") +
 theme(plot.subtitle=element_text(size=9), legend.position="top")+
 scale_color_manual(name="Rok 2020: ", labels = c("% 2020/(średnia 2015--2015)"),
   values = c("yy"=farbe20 )  ) +
 ggtitle("Zgony wg województw* (PL/Ogółem)", 
   subtitle=sprintf("*wg klasyfikacji NUTS stąd mazowieckie/stołeczne | %s | ostatni tydzień: %s", NIKW, lastO))

ggsave(plot=p5, "zgony_PL_by_woj_P.png", width=picWd)

zz1 <- zz1 %>% filter ( week >= firstWeek  ) %>% as.data.frame

p6 <- ggplot(zz1, aes(x=week, y=yy), color=farbe20) +
 geom_smooth(method="loess", se=F, span=spanV, size=.4, color=farbe20) +
 geom_point(size=.4, alpha=.5) +
 geom_text(aes(label=sprintf("%.1f", yy)), vjust=-1.25, size=1.5) +
 facet_wrap( ~name, scales = "fixed") +
 xlab(label="nr tygodnia") +
 ylab(label="%") +
 theme(plot.subtitle=element_text(size=9), legend.position="top")+
 scale_color_manual(name="Rok 2020: ", labels = c("% 2020/(średnia 2015--2015)"),
   values = c("yy"=farbe20 )  ) +
   ggtitle(sprintf("Zgony wg województw* (PL/Ogółem) tygodnie: %i--%i (%i tydzień zaczyna się %s)",
     firstWeek, lastWeek, lastWeek, lastO), 
   subtitle=sprintf("*wg klasyfikacji NUTS stąd mazowieckie/stołeczne | %s", NIKW))

ggsave(plot=p6, "zgony_PL_by_woj_P6.png", width=picWd)

Pojechaliśmy w tą sobotę z MK zbierać bursztyn do Mikoszewa. Ponieważ miałem zajęcia, to dotarliśmy tam dopiero około 17:00 czyli było już ciemno. Okazało się że jesteśmy trochę późno, bo wysyp był 24h wcześniej (ale chyba niezbyt obfity)...

My zastaliśmy masę patyków na brzegu, a w nich miliony kawałków bursztynu, tyle że takie nalewkowe i mniejsze. Przez dwie godziny nazbieraliśmy tego prawie 600g (ja tym razem więcej bo 320g a MK 260g.)

Wskaźnik zapadalności (dalej WuZet) dla ostatnich 14 dni po naszemu w oryginale zaś 14-day notification rate of newly reported COVID-19 cases per 100 000 population (Data on 14-day notification rate of new COVID-19 cases and deaths). Nasz wspaniały rząd walczący z COVIDem przyjął w swoim czasie `nową strategię' (od tego czasu już kilka razy odnawianą), a w niej był podział na strefy zielona/żółta/czerwona, definiowane odpowiednio jako wartości WuZet poniżej 6 (zielona); 6--12 (żółta) oraz powyżej 12 (czerwona). Dla Sopotu na przykład, który oficjalnie ma około 35 tys mieszkańców, do wejście do czerwonej strefy wystarczały zatem około 4 zakażenia (w ciągu dwóch tygodni). To wszysto wydaje się dziś śmieszne jak ostatnio zakażeń dziennie potrafi być na poziomie trzy-tygodniowej dawki...

Parę dni temu mój bank danych nt. COVID19 uzupełniłem o dane dotyczące liczby zakażeń w powiatach województwa pomorskiego. Akurat WSSE w Gdańsku takie dane w sposób w miarę porządny publikuje i da się je względnie łatwo odzyskać z raportów publikowanych i archiwizowanych (brawo--na to nie wpadł nawet Minister w MZ) pod adresem http://www.wsse.gda.pl/.

library("dplyr")
library("tidyr")
library("ggplot2")
m1unit <- 100000
# pomorskie_by_powiat.csv pobrane z www.wsse.gda.pl
# format: data;powiat;nc 
d <- read.csv("pomorskie_by_powiat.csv", sep = ';',  header=T, na.string="NA" )

# wartości skumulowane (po powiatach dlatego tak dziwnie)
# replace_na jest potrzebne bo cumsum nie obsługuje NA
e <- d %>% group_by(powiat) %>% dplyr::mutate(tc = cumsum(replace_na(nc, 0)),
tc1m = cumsum(replace_na(nc, 0)) / pop * m1unit
)
## dzień ostatni
day00 <- as.Date(last.obs)
## dwa tygodnie przed ostatnim
day14 <- as.Date(last.obs) - 14
## 4 tygodnie przed ostatnim
day28 <- as.Date(last.obs) - 28

## Stan na dzień ostatni
e00 <- e %>% filter (as.Date(dat) == day00 ) %>%
 group_by(powiat) %>% as.data.frame

## BTW Żeby było dziwniej zapis
## e0 <- e %>% group_by(powiat) %>%  slice_tail(n=1) %>% as.data.frame
## Daje inny wynik, liczby te same, ale porządek wierszy jest inny

## Stan na dzień dwa tygodnie przed
e14 <- e %>% filter (as.Date(dat) == day14 ) %>%
 group_by(powiat) %>% as.data.frame

e14 <- e %>% filter (as.Date(dat) == day14 ) %>%
 group_by(powiat) %>% as.data.frame

## Stan na dzień 4 tygodnie przed
e28 <- e %>% filter (as.Date(dat) == day28 ) %>%
 group_by(powiat) %>% as.data.frame

## nowe zakażenia w tygodniach 3/4
c28 <- e14$tc - e28$tc

## nowe zakażenie w tygodniach 2/1
c14 <- e00$tc - e14$tc
## To samo co c14/c28 ale w przeliczeniu
## na 100 tys:
c28m1 <- e14$tc1m - e28$tc1m
c14m1 <- e00$tc1m - e28$tc1m
## Dynamika zmiana WuZet w dwóch ostatnich okresach
## tj WuZet12 względem WuZet34
#d14v28 <- (c14m1 - c28m1) / c28m1 * 100
d14v28 <- (c14m1/c28m1) * 100
##
## Można sobie teraz c14m1/d14v28 na wykresach przestawić

Na dzień 7 listopada wyniki były takie (jeżeli się nie rąbnąłem w powyższym kodzie):

Na dzień 7 listopada wyniki były takie (jeżeli się nie rąbnąłem w powyższym kodzie):

sprintf("%10.10s = %6.2f | %6.2f | %6.2f - %6.2f | %4i | %4i | %4i (%4i)",
   e00$powiat, d14v28, c14m1, e00$tc1m, e28$tc1m, e00$tc, e14$tc, e28$tc, e00$pop )
 [1] "    Gdynia = 229.15 | 577.74 | 1129.08 - 299.22 | 2781 | 1358 |  737 (246306)"
 [2] "    Gdańsk = 149.50 | 416.37 | 1045.98 - 351.10 | 4856 | 2923 | 1630 (464254)"
 [3] "    Słupsk = 228.26 | 803.59 | 1387.42 - 231.78 | 1269 |  534 |  212 (91465)"
 [4] "     Sopot = 144.90 | 583.03 | 1404.21 - 418.80 |  513 |  300 |  153 (36533)"
 [5] "  tczewski = 437.50 | 905.98 | 1323.60 - 210.53 | 1534 |  484 |  244 (115896)"
 [6] "   gdański = 399.21 | 889.61 | 1353.71 - 241.26 | 1543 |  529 |  275 (113983)"
 [7] "  kartuski = 197.07 | 855.49 | 1846.22 - 556.63 | 2471 | 1326 |  745 (133841)"
 [8] "  bytowski = 268.71 | 998.31 | 1693.33 - 323.50 | 1340 |  550 |  256 (79134)"
 [9] " malborski = 329.61 | 923.16 | 1408.21 - 204.97 |  900 |  310 |  131 (63911)"
[10] "     pucki = 225.35 | 766.35 | 1545.69 - 439.26 | 1309 |  660 |  372 (84687)"
[11] "wejherowsk = 150.71 | 396.16 |  971.92 - 312.90 | 2078 | 1231 |  669 (213803)"
[12] "starogardz = 216.36 | 744.52 | 1388.93 - 300.31 | 1776 |  824 |  384 (127868)"
[13] " chojnicki = 266.33 | 813.36 | 1311.04 - 192.29 | 1275 |  484 |  187 (97251)"
[14] "  sztumski = 309.52 | 619.84 |  979.83 - 159.73 |  411 |  151 |   67 (41946)"
[15] " kwidzyńsk = 251.34 | 563.39 |  973.35 - 185.80 |  812 |  342 |  155 (83423)"
[16] " kościersk = 293.30 | 786.88 | 1392.60 - 337.43 | 1007 |  438 |  244 (72311)"
[17] "nowodworsk = 263.21 | 777.35 | 1273.30 - 200.61 |  457 |  178 |   72 (35891)"
[18] "   słupski = 244.23 | 514.50 |  969.24 - 244.08 |  957 |  449 |  241 (98737)"
[19] "  lęborski = 172.27 | 618.09 | 1135.18 - 158.29 |  753 |  343 |  105 (66333)"
[20] " człuchows = 268.29 | 388.16 |  571.65 -  38.82 |  324 |  104 |   22 (56678)"

Dynamika WuZeta: Gdynia/Gdańsk/Sopot = 129.1%, 149.5% i 144.9%; wartość maksymalna 437% (tczewski); wartość minimalna 150% (wejherowski). Najnowsze wartości współczynnika WuZet: Gdynia/Gdańsk/Sopot= 577.7, 416.4 oraz 583.0; wartość maksymalna 998,3 (bytowski); wartość minimalna 388.1 (człuchowski). Dane i skrypty są tutaj.

Oryginalny plik prognozy

Wycięty fragment

`Podglądówka' całej prognozy

Są na stronie http://www.meteo.pl/ udostępnione 120h prognozy falowania Bałtyku w postaci rysunków. Na tych rysunkach wysokość fal jest odwzorowana kolorem od czerwonego do jasno niebieskiego (najmniejsze). Kierunek fal jest odwzorowany strzałkami. Są też dostępne dane historyczne.

Wszystko jest łatwe do pobrania, bo nazwa pliku historycznego to na przykład wavehgt_0W_2015122200_015-00.png czyli wavehgt_0W_YYYYMMDD00_0HH-00.png (gdzie HH to 03/09/15/21). Pliki prognoz mają też nazwy wg schematu: wavehgt_0W_YYYYMMDD00_HHH-00.png, tyle że HHH się zmienia w zakresie 12--120 co trzy (godziny).

Pobieram wgetem 8 plików 12--120. Teraz trzeba ustalić jakie są kolory na obrazku i wysłać komunikat (analiza strzałek to beznadzieja sprawa, nawet się nie zabieram.)

Zaczynam od wycięcia interesującego mnie fragmentu Bałtyku:

## 8 plików 12--120
convert wavehgt_0W_YYYYMMDD00_HHH-00.png-00.png -crop 65x86+280+575 \
 wavehgt_0W_YYYYMMDD00_HHH-00.png-00D.png
## Wypisz kolory z rysunku:
convert wavehgt_0W_YYYYMMDD00_HHH-00.png-00D.png \
  -define histogram:unique-colors=true -format %c histogram:info:- | \
  getMaxFala.pl 
## Łączę w jeden:
montage PLIKI... -tile 2x4 -border 4 -geometry 480  S_YYYMMDD.png

Plik getMaxFala.pl zwraca kolor odpowiadający za najwyższą falę. Szczęśliwie tych kolorów za dużo nie jest, więc sprawdzanie jest enumeratywne

Mając najwyższe fale w horyzoncie 12--120 ustalamy max z tych 8 liczb i wypisujemy komunikat (Perl):

print "Max: $maxFF expected (in 120h time window)\n";
print "Details: ";
for $f (sort keys %faleMax) { print "$f = $faleMax{$f} /" }
print "\nAmber likelihood: ";
if ($maxFF < 3) { print "NONE\n" } 
elsif ($maxFF < 4) { print "VERY SMALL\n" }
elsif ($maxFF < 5) { print "SMALL\n" }
elsif ($maxFF < 5) { print "MEDIUM\n" }
elsif ($maxFF < 6) { print "LARGE\n" }
else  { print "Amber: HUGE"}
print "http://pinkaccordions.homelinux.org/fale/forecast/S_${yyyymmdd}.png\n";
print "Bye...\n";

Prawdę powiedziawszy polecenia convert/mogrify/wget też są `wbudowane' w skrypt Perla. Całe zadanie realizuje jeden skrypt, który w efekcie wypisuje na ekran powyższy komunikat:

  ## Fragment skryptu
  $GETMAXFALA='/home/pi/bin/getMaxFala.pl';
  ## ... ##
for $h (@Hours) {
    $hrNo = sprintf "%03.3i", $h;
    $url= "$URL/${yyyymmdd}00/wavehgt_0W_${yyyymmdd}00_${hrNo}-00.png";
    ###print STDERR "$url\n";
    system ("wget $url -O W_${yyyymmdd}00_${hrNo}-00.png");
    system ("convert W_${yyyymmdd}00_${hrNo}-00.png -crop 280x330+125+350 -fill red" .
        " -annotate +140+360 '$monthNo/${dayNo}+${hrNo}' -pointsize 24 -fill blue" .
        " W_${yyyymmdd}00_${hrNo}-00_C.png");
    #### Miniatura zatoki
    system ("convert W_${yyyymmdd}00_${hrNo}-00.png -crop 65x86+280+575 W_${yyyymmdd}00_${hrNo}-00_D.png");
    $maxFala = `convert "W_${yyyymmdd}00_${hrNo}-00_D.png" -define histogram:unique-colors=true
       -format %c histogram:info:- | $GETMAXFALA`;
    chomp($maxFala); 
    $files .= "W_${yyyymmdd}00_${hrNo}-00_C.png ";
    $files_S .= "W_${yyyymmdd}00_${hrNo}-00_D.png ";
    $faleMax{"+${hrNo}"}=$maxFala;
}

system ("montage $files -tile 2x4 -border 4 -geometry 480  W_${yyyymmdd}.png" );
system ("montage $files_S -tile 2x4 -border 4 -geometry 480  S_${yyyymmdd}.png" );
$maxFF = max(values(%faleMax)); chomp($maxFF);
  ## .. ##

Ten komunikat (w potoku) czyta inny skrypt i wysyła listy do zainteresowanych.

Skrypt rysujący wykres łącznej liczby zakażeń i zgonów wg kontynentów. Przy okazji podjąłem nieśmiajłą próbę zapanowania nad (jednolitym) wyglądem. Udając, że istnieje coś takiego jak Narodowy Instytut Korona Wirusa stworzyłem dla niego wizualizację, którą dodaje się do wykresu w postaci jednego polecenia. To jedno polecenie to funkcja theme_nikw:

#!/usr/bin/env Rscript
#
library("ggplot2")
library("dplyr")
library("scales")
library("ggthemes")
library("ggpubr")
#
theme_nikw <- function(){
 theme_wsj() %+replace%
 theme(
  panel.grid.minor = element_line(size = 0.25, linetype = 'solid', colour = "cornsilk3"),
  panel.grid.major = element_line(size = 0.25, linetype = 'solid', colour = "cornsilk3"),
  ##
  axis.text.x  = element_text(size = 6 ),
  axis.text.y  = element_text(size = 6 ),
  ## https://stackoverflow.com/questions/14379737/how-can-i-make-xlab-and-ylab-visible-when-using-theme-wsj-ggthemes
  axis.title  = element_text(family="sans", size=6)
  ## Poniższe natomiast nie działa:
  #axis.title.x  = element_text(size = 6 ),
  #axis.title.y  = element_text(size = 6 ),
  ## margin=margin(r=, t=, b = 3, l=, unit="pt")
  plot.title=element_text(family="sans", size=14, hjust=0, margin=margin(b = 3, unit="pt")),
  plot.subtitle=element_text(family="sans", size=8, hjust=0),
  legend.title=element_text(family="sans", size=8),
  plot.caption = element_text(family="sans", size = 6)
  )
}

Uprzedzając wydarzenia, kolejna funkcja notin służy do tego co funkcja in tyle, że zamiast zwracać wartości pasujące, zwraca te które nie pasują:

## https://stackoverflow.com/questions/34444295/how-to-specify-does-not-contain-in-dplyr-filter
`%notin%` = function(x,y) !(x %in% y)

Teraz ustawiam różne wartości globalne (options służy do wyłączenia zapisu liczb w trybie scientic notation).

options(scipen = 999) ## albo options(scipen = 999, digits=3)
note <- "(c) NI-KW @ github.com/knsm-psw/NI-KW"
today <- Sys.Date()
tt <- format(today, "%Y-%m-%d")
spanV <- 0.25
mainBreaks <- "4 weeks"
mainColor <- "deeppink"
loessColor <- "deeppink"

## Przed/po \n musi być odstęp inaczej nie działa
## url <- "https://www.ecdc.europa.eu/en/publications-data/ \n download-todays-data-geographic-distribution-covid-19-cases-worldwide"
url <- "www.ecdc.europa.eu/en/publications-data/download-todays-data-geographic-distribution-covid-19-cases-worldwide"
surl <- sprintf ("Data retrived from %s", url)

Czytanie zbiorów danych. Dane dotyczące zakażeń pochodzą ze strony ECDC (European Center for Disease Prevention and Control). Dane dotyczące ludności też z tej strony tyle, że wydłubałem je raz i trzymam w oddzielnym pliku zamiast powielać niepotrzebnie w pliku z danymi o COVID19 (jak to robi ECDC)

## info o wielkości populacji i kontynencie
c <- read.csv("ecdc_countries.csv", sep = ';',  header=T, na.string="NA" )
## dane dzienne/cały świat (ECDC)
d <- read.csv("covid19_C.csv", sep = ';',  header=T, na.string="NA", 
   colClasses = c('factor', 'factor', 'factor', 'character', 'character', 'numeric', 'numeric'));

Złączenie zbiorów w oparciu o kolumnę id:

d <- left_join(d, c, by = "id")

Parę rzeczy trzeba uporządkować. Po pierwsze usunąć wpisy zawierające kraj bez ISO-kodu (coś jest nie tak z transformacją XSLX na CSV, muszę to sprawdzić). Po drugie czasami liczby zakażeń/śmierci są ujemne bo dane zostają zrewidowane. Nie jest to najlepsze rozwiązanie, ale żeby przynajmniej na pierwszy rzut oka wykres nie wyglądał absurdalnie to zamieniam liczby ujemne na NA:

# czyszczenie
# tylko kraje z ID
d <- d %>% filter(!is.na(id))
# zamień na liczby (był problem z czytaniem CSV bezpośrednio)
d$newc <- as.numeric(d$newc)
d$newd <- as.numeric(d$newd)

# Dane bezsensowne zamień na NA
# Zdarzają się ujemne bo kraje zmieniają klasyfikację/przeliczają
d$newc[ (d$newc < 0) ] <- NA
d$newd[ (d$newd < 0) ] <- NA

Liczę wartości na 1/mln. Obliczam dzienne wartości zakażeń/śmierci na 1mln wg kontynentów

## zgony na 1mln
d$tot1m <- d$totald/d$popData2019 * 1000000
## Oblicz sumy dla kontynentów
## Kontynentów jest więcej niż 5:-) w zbiorze (jakieś śmieci)
continents <- c('Africa', 'America', 'Asia', 'Europe', 'Oceania')
## ogółem wg 
cont <- d %>% filter (continent %in% continents) %>% group_by(date,continent) %>% summarise(
  tc = sum(newc, na.rm=TRUE),
  td = sum(newd, na.rm=TRUE),
  tc1m = sum(newc, na.rm=TRUE)/sum(popData2019, na.rm=TRUE) * 1000000,
  td1m = sum(newd, na.rm=TRUE)/sum(popData2019, na.rm=TRUE) * 1000000
  ) %>% as.data.frame
## str(tc)

## Z ciekawości
noncont <- d %>% filter (continent %notin% continents) %>% as.data.frame
print(noncont)

Wykreślam wykresy. Punkty przestawiają dane dzienne. Krzywa loess trend. Źródło danych jest w podtytule. W napisie umieszczanym pod wykresem jest informacja o autorze i link do repozytorium githuba (polecenie labs(caption=note))

pd1 <- ggplot(cont, aes(x= as.Date(date, format="%Y-%m-%d"), color = continent, y=tc)) +
 geom_point(aes(y=tc, color = continent), size=.4, alpha=.5) +
 geom_smooth(method="loess", se=F, span=spanV, aes(fill=continent, y=tc, group=continent)) +
 xlab(label="") + ylab(label="cases") + 
 scale_x_date( labels = date_format("%m/%d"), breaks = mainBreaks) +
 theme_nikw() +
 labs(caption=note) +
 ggtitle(sprintf ("COVID19: cases (%s)", last.obs), subtitle=sprintf("%s", surl))

pd2 <- ggplot(cont, aes(x= as.Date(date, format="%Y-%m-%d"), , color = continent, y=td)) +
 geom_point(aes(y=td, color = continent), size=.4, alpha=.5) +
 geom_smooth(method="loess", se=F, span=spanV, aes(fill=continent, y=td, group=continent)) +
 xlab(label="") + ylab(label="deaths") + 
 scale_x_date( labels = date_format("%m/%d"), breaks = mainBreaks) +
 theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
 theme_nikw() +
 labs(caption=note) +
 ggtitle(sprintf ("COVID19: deaths (%s)", last.obs), subtitle=sprintf("%s", surl))

ggsave(plot=pd1, file="tc_by_c.png", width=10)
ggsave(plot=pd2, file="td_by_c.png", width=10)

pd1m <- ggplot(cont, aes(x= as.Date(date, format="%Y-%m-%d"), color = continent, y=tc1m)) +
 geom_point(aes(y=tc1m, color = continent), size=.4, alpha=.5) +
 geom_smooth(method="loess", se=F, span=spanV, aes(y=tc1m )) +
 xlab(label="") + ylab(label="cases") + 
 scale_x_date( labels = date_format("%m/%d"), breaks = mainBreaks) +
 theme_nikw() +
 labs(caption=note) +
 ggtitle(sprintf ("COVID19: cases per 1mln (%s)", last.obs), subtitle=sprintf("%s", surl))

pd2m <- ggplot(cont, aes(x= as.Date(date, format="%Y-%m-%d"), color = continent, y=td1m)) +
 geom_point(aes(y=td1m, color = continent), size=.4, alpha=.5) +
 geom_smooth(method="loess", se=F, span=spanV, aes(y=td1m )) +
 xlab(label="") + ylab(label="deaths") + 
 scale_x_date( labels = date_format("%m/%d"), breaks = mainBreaks) +
 theme_nikw() +
 labs(caption=note) +
 ggtitle(sprintf ("COVID19: deaths per 1mln(%s)", last.obs), subtitle=sprintf("%s", surl))

ggsave(plot=pd1m, file="tc1m_by_c.png", width=10)
ggsave(plot=pd2m, file="td1m_by_c.png", width=10)

Wszystko. Skrypt jest wykonywany automatycznie o ustalonej godzinie, po czym wykresy wysyłane są na twitterowe konto "Instytutu".

Taka książka jest pn 'Koniec PiSu' (Mrozowski rozmawia z Kamińskim). Całkiem nadająca się do czytania BTW, bo tak jak nie cenię autorów, to w tym przypadku nie poszli w tanie hejterstwo. Książka z 2012 roku, okazawszy się raczej marnym proroctwem, stała się przedmiotem śmichów-chichów, ale wygląda że -- być może -- co się odwlecze to nie uciecze.

Od razu zaznaczę, że od zawsze popieram PiS a młodszym przypomnę, że do roku 2015 była to beznadziejna sprawa--zawsze wpiernicz i nigdy nie miałem wątpliwości na kogo głosować. Do teraz...

Teraz bym nie zagłosował, co więcej nawet gdyby wybory się odbyły za trzy lata to na dziś wygląda że też bym nie zagłosował. Nie to że mi się spodobał p. Budka czy p. Hołownia. Nie, nie głosowałbym wcale. Moi idole przegięli pałę...

Po pierwsze i najważniejsze: nie potrafię wytłumaczyć co konkretnie chciał osiągnąć Jarosław Kaczyński proponując ustawę w/s ochrony praw zwierząt. Nie można z piątku na sobotę powiedzieć ludziom i to ludziom ze środowiska, bez którego nie byłoby zwycięstwa partii PiS w ostatnich wyborach: Mamy w dupie czy wzieliście kredyty na 5 lat czy na 15 lat, macie rok na to żeby się zwijać. Nie mogę tego nazwać inaczej niż bezczelnością/butą/oderwaniem się od rzeczywistości. Nie znajduję też innego wytłumaczenia, że taki pomysł w ogóle się pojawił, jak uwiąd intelektualny JK...

Ustawa oprócz tego, że szkodzi interesom ekonomicznym wsi (i nie tylko) jest głupia w każdym praktycznie aspekcie, oprócz tego, pod którym jest szerzej znana. Bo znana jest jako futerkowa czyli zakazująca hodowli zwierząt futerkowych, tyle że akurat zakaz hodowli zwierząt futerkowych byłby OK (pod warunkiem oczywiście wyrównania strat ekonomicznych hodowcom). To co uważam za największy idiotyzm `ustawy futerkowej', to zakaz uboju rytualnego. Lewacka logika 101: multikulti ale bez halal. Bardzo was lubimy, ale tego/tamtego wam nie wolno. Wygląda że na starość p. Bóg odebrał rozum prezesowi. Nie jest naszą sprawą ocena/zakazy tego co 1/3 ludzkości uważa za normalne.

Po drugie: tzw. rekonstrukcja rządu. Zapowiadana, odwlekana, wreszcie zrealizowana w jakimś dziwnym trybie. Świetnych ministrów usunięto, zastępując ich jeszcze lepszymi (ironia); podpisano uroczyste i kuriozalne porozumienie. Kuriozalne bo przecież nie programowe--program był ogłoszony w 2019 a my zawsze dotrzymujemy słowa. Więc w TV pokazano uroczyste porozumienie w/s podziału stanowisk. Tzn. pokazano jak podpisują, co podpisują już nie pokazano. No wspaniale... Bezczelność i buta #2.

A że p. Bóg karze za głupotę i bezczelność to ledwo poprzydzielano fotele nowym/starym (złotousty-Gowin) zaczęła się pandemia. Gołym okiem widać panikę i brak planów tych nowych wspaniałych jeszcze lepszych od starych wspaniałych. W lecie kiedy był czas na przygotowanie się (choćby w trybie symulacji -- co zrobić jak będzie druga fala epidemii) to rządzący zajęci byli `rekonstrukcją'/ustawą futerkową. Teraz improwizują.

Jeden tylko przykład: do 29 października zapowiadają w telewizjach, że nie będzie zakazu wstępu na cmentarze; 30 października premier oznajmia że będzie zakaz. Kłamie przy tym bezczelnie, że mieli nadzieję, że pandemia wygaśnie (na jakiej podstawie? Że się cud Boski zdarzy?) Potem widocznie komuś w KPRM się przypomina, że na lodzie została kolejna grupa przedsiębiorców a premier nic im nie obiecał, więc wieczorem naprędce wydają uzupełniający komunikat, że im się też zrekompensuje, ale w jaki sposób konkretnie, to będzie później (bo teraz to sami nie wiemy?). Następnego dnia jest komunikat w/s tego później: Agencja Restrukturyzacji i Modernizacji Rolnictwa będzie skupować kwiaty (są już żarty, że w grudniu będą odkupować choinki i karpie). Ponieważ większość tych przedsiębiorców to drobni handlarze prowadzący działalność nierejestrowaną/ryczałt/handel obwoźny, to ciekawe jak Agencja ustali komu/ile się należy. Kto bezczelny dostanie, kto uczciwy obejdzie się smakiem?

Jak to nie jebnie to będzie cud-boski. Ale jak ten cud się zrealizuje to będzie tylko trwanie, ten rząd już niczego nie zrealizuje, więc nie wiadomo co lepsze: trwanie czy upadek.

Nawiasem i na koniec: akurat protesty po wyroku TK ws aborcji eugenicznej (wbrew wielkim nadziejom antyPiSu) mają -- moim zdaniem -- niewielkie szanse na sukces. Gołym okiem widać głupotę protestujących (wulgarność, atakowanie kościołów i podobnych obiektów, które należy zostawić w spokoju a demonstrować gdzie indziej, bezsensowne blokowanie ruchu na ulicach itd/itp). Sama szefowa, bliżej mi nieznana 41-letnia Marta Wypierdalać-Lampart wygląda na osobę mocno intelektualnie ograniczoną za to o dużych ambicjach (nieciekawy mix jak wiadomo). Sam fakt że w wieku 41-lat niewiele osiągnęła wskazuje, że coś z nią jest nie tak. A co konkretnie to niedługo się okaże, jak to wszystko pizdnie. Ot KOD/Kijowski season 2020, skończą podobnie...

Dopisek z ostatniej chwili: ustawa futerkowa z powodu której: 1) grożono rozpadem koalicji; 2) usunięto ze stanowiska ministra rolnictwa; 3) zawieszono iluś tam posłów, idzie do kibla. Czyli była do dupy. Czyli `nieomylny/genialny' prezes...?