Weblog Tomasza Przechlewskiego [Zdjęcie T. Przechlewskiego]


scrum
random image [Photo gallery]
Zestawienie tagów
1-wire | 18b20 | 1wire | 2140 | 3rz | adamowicz | afera | alsamixer | amazon | amber | amman | anniversary | antypis | apache | api | applebaum | arm | armenia | astronomy | asus | atom.xml | awk | aws | bachotek | bakłażan | balcerowicz | balta | banan | bash | batumi | berlin | białowieża | białystok | bibtex | bieszczady | biznes | blogger | blogging | blosxom | bme280 | bono | borne-sulinowo | breugel | bt747 | budapeszt | budyń | bursztyn | campagnolo | canon | cedewu | chaos | chello | chiller | chillerpl | chown | christophe dominici | chujowetaśmy | ciasto | cmentarz | contour | coronavirus | covi19 | covid | covid19 | cron | css | csv | cukinia | curl | cycling | d54250wykh | dbi | debian | dejavu | dhcp | dht22 | dia | docbook | dom | dp1500 | ds18b20 | duda | dulkiewicz | dyndns | dynia | ebay | economy | ecowitt | ekonomia | elka | elm | emacs | emacs23 | english | ep | erasmus | erasmusplus | ess | eu | eurostat | excel | exif | exiftool | f11 | fc | fc11 | fc15 | fc29 | fc5 | fc8 | fedora | fedora21 | fenix | ffmpeg | finepix | firefox | flickr | folau | fontforge | fontspec | fonty | food | fop | forms | foto | france | francja | fripp | froggit | fuczki | fuji | fuse | gammu | garden | garmin | gas | gawk | gazwyb | gdańsk | gdynia | gender | geo | geocoding | georgia | gft | ggplot | git | github | gmail | gmaps | gnokii | gnus | google | google apps script | googlecl | googleearth | googlemaps | gotowanie | gphoto | gphoto2 | gps | gpsbabel | gpsphoto | gpx | gpx-viewer | greasemonkey | gruzja | grzyby | gus | gw1000 | haldaemon | handbrake | hhi | historia | history | hitler | holocaust | holokaust | hp1000se | hpmini | humour | iblue747 | ical | iiyama | ikea | imagemagick | imap | inkscape | inne | internet | j10i2 | javascript | jhead | jordania | k800i | kajak | kamera | karob | kibbeh | kleinertest | kml | kmobiletools | knuth | kociewie kołem | kod | kolibki | komorowski | konwersja | krutynia | krynki | kuchnia | kurski | kłamstwo | latex | latex2rtf | latex3 | lcd | legend | lenny | lesund | lewactwo | lgbt-folly | liban | liberation | linksys | linux | lisp | lisrel | litwa | lizbona | logika | ltr | lubowla | lwp | lwów | m2wś | malta | mapquest | mapsource | maradona | marchew | marimekko | marvell | math | mathjax | mazury | mbank | mediolan | mencoder | mevo | mh17 | michalak | michlmayr | microsoft | monitor | mp4box | mplayer | ms | msc | mssql | msw | mswindows | mtkbabel | museum | muzyka | mymaps | mysql | mz | nafisa | nanopi | natbib | navin | nekrolog | neo | neopi | netbook | niemcy | niemieckie zbrodnie | nikon | nmea | nowazelandia | nuc | nxml | oauth | oauth2 | obituary | ocr | odessa | okular | olympus | ooffice | ooxml | opera | osm | otf | otftotfm | other | overclocking | ozbekiston | panoramio | paryż | pdf | pdfpages | pdftex | pdftk | pedophilia | perl | photo | photography | picasa | picasaweb | pim | pine | pis | pit | plotly | pls | plugin | po | podcast | podlasie | podróże | pogoda | politics | polityka | polsat | portugalia | postęp | powerpoint | połtawa | prelink | problem | propaganda | pseudointeligencja | pstoedit | putin | python | pywws | r | r1984 | radio | random | raspberry | raspberry pi | raspberrypi | raspbian | refugees | relaxng | ridley | router | rower | rowery | rpi | rsync | rtf | ruby | rugby | rumunia | russia | rwc | rwc2007 | rwc2011 | rwc2019 | rzym | samba | sds011 | selenium | sem | senah | sernik | sheevaplug | sienkiewicz | signature | sks | skype | skytraq | smoleńsk | sqlite | srtm | sshfs | ssl | staszek wawrykiewicz | statistcs | statistics | stats | statystyka | stix | stretch | supraśl | suwałki | svg | svn | swanetia | swornegacie | szwajcaria | słowacja | tbilisi | terrorism | tesseract | tex | texgyre | texlive | thunderbird | tomato | totalnaopozycja | tourism | tramp | trang | transylwania | truetype | trzaskowski | ttf | turcja | turkey | turystyka | tusk | tv | tv5monde | tweepy | twitter | tykocin | typetools | ubuntu | uchodźcy | udev | ue | ukraina | umap | unix | upc | updmap | ups | utf8 | uzbekistan | varia | video | vienna | virb edit | virbedit | vostro | wammu | wdc | wdfs | weather | weathercloud | webcam | webdav | webscrapping | weewx | wh2080 | wiedeń | wikicommons | wilno | win10 | windows | windows8 | wine | wioślarstwo | word | wordpress | wrt54gl | ws1080 | wtyczka | wunderground | ww2 | www | wybory | wybory2015 | włochy | węgry | xemex | xetex | xft | xhtml | xine | xml | xmllint | xsd | xslt | xvidtune | youtube | yum | zaatar | zakopane | zakupy | zawodzie | zdf | zdrowie | zgarden | zgony | zprojekt | łeba | świdnica | żywność
Archiwum
O stronie
wykorzystywany jest blosxom plus następujące wtyczki: tagging, flatarchives, rss10, lastbuilddatexhtmlmime. Niektóre musiałem dopasować nieco do swoich potrzeb. Więcej o blosxom jest tutaj
Subskrypcja
RSS 1.0
Tygodniowe dane nt. zgonów z GUS






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)

url | Mon, 23/11/2020 17:26 | tagi: , , , ,
COVID19: współczynnik WZ dla powiatów województwa pomorskiego

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.

url | Sun, 08/11/2020 06:09 | tagi: , ,
Covid19 wg kontynentów


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

url | Mon, 02/11/2020 08:47 | tagi: , , ,