Monkeypox Map
Dhihram Tenrisau
2022-08-13
Packages
library(rgdal)
library(RCurl)
library(tidyverse)
library(dplyr)
library(broom)
library(sp)
library(ggmap)
library(viridis)
library(scales)
library(mapsf)
library(sf)
library(RColorBrewer)
library(plotly)Data Preparation
Import the Map
In this project, I will use the shp file from here
setwd("C:/Users/dhihr/Downloads/template paskhas")
dunia <- readOGR(dsn = 'world-administrative-boundaries' , layer = 'world-administrative-boundaries')## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\dhihr\Downloads\template paskhas\world-administrative-boundaries", layer: "world-administrative-boundaries"
## with 256 features
## It has 8 fieldssummary(dunia@data)## status color_code region iso3
## Length:256 Length:256 Length:256 Length:256
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
## continent name iso_3166_1_ french_shor
## Length:256 Length:256 Length:256 Length:256
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :characterUsing Our World in Data Dataset
This data will use the dataset from Our World in Data
x <- getURL("https://raw.githubusercontent.com/owid/monkeypox/main/owid-monkeypox-data.csv")
MPX <- read.csv(text = x)
MPX <- rename(MPX, name = location)
MPX[is.na(MPX)] = 0
summary(MPX)## name date new_cases new_cases_smoothed
## Length:4004 Length:4004 Min. : 0.00 Min. : 0.000
## Class :character Class :character 1st Qu.: 0.00 1st Qu.: 0.000
## Mode :character Mode :character Median : 0.00 Median : 0.570
## Mean : 19.34 Mean : 17.977
## 3rd Qu.: 1.00 3rd Qu.: 4.035
## Max. :1998.00 Max. :1014.140
## total_cases new_cases_per_million total_cases_per_million
## Min. : 1.0 Min. : 0.0000 Min. : 0.000
## 1st Qu.: 3.0 1st Qu.: 0.0000 1st Qu.: 0.270
## Median : 13.0 Median : 0.0000 Median : 1.116
## Mean : 507.9 Mean : 0.2948 Mean : 7.000
## 3rd Qu.: 101.2 3rd Qu.: 0.0300 3rd Qu.: 6.949
## Max. :38709.0 Max. :54.5170 Max. :121.970
## new_cases_smoothed_per_million new_deaths new_deaths_smoothed
## Min. :0.0000 Min. :0.000000 Min. :0.000000
## 1st Qu.:0.0000 1st Qu.:0.000000 1st Qu.:0.000000
## Median :0.0410 Median :0.000000 Median :0.000000
## Mean :0.2196 Mean :0.005994 Mean :0.003147
## 3rd Qu.:0.2590 3rd Qu.:0.000000 3rd Qu.:0.000000
## Max. :7.9050 Max. :2.000000 Max. :1.000000
## total_deaths new_deaths_per_million total_deaths_per_million
## Min. : 0.0000 Min. :0.00e+00 Min. :0.0000000
## 1st Qu.: 0.0000 1st Qu.:0.00e+00 1st Qu.:0.0000000
## Median : 0.0000 Median :0.00e+00 Median :0.0000000
## Mean : 0.1006 Mean :4.77e-05 Mean :0.0001843
## 3rd Qu.: 0.0000 3rd Qu.:0.00e+00 3rd Qu.:0.0000000
## Max. :12.0000 Max. :5.60e-02 Max. :0.0560000
## new_deaths_smoothed_per_million
## Min. :0.000e+00
## 1st Qu.:0.000e+00
## Median :0.000e+00
## Mean :2.248e-06
## 3rd Qu.:0.000e+00
## Max. :1.000e-03head(MPX)## name date new_cases new_cases_smoothed total_cases
## 1 Andorra 2022-07-25 1 0 1
## 2 Andorra 2022-07-26 2 0 3
## 3 Andorra 2022-07-27 0 0 3
## 4 Andorra 2022-07-28 0 0 3
## 5 Andorra 2022-07-29 0 0 3
## 6 Andorra 2022-07-30 0 0 3
## new_cases_per_million total_cases_per_million new_cases_smoothed_per_million
## 1 12.653 12.653 0
## 2 25.306 37.958 0
## 3 0.000 37.958 0
## 4 0.000 37.958 0
## 5 0.000 37.958 0
## 6 0.000 37.958 0
## new_deaths new_deaths_smoothed total_deaths new_deaths_per_million
## 1 0 0 0 0
## 2 0 0 0 0
## 3 0 0 0 0
## 4 0 0 0 0
## 5 0 0 0 0
## 6 0 0 0 0
## total_deaths_per_million new_deaths_smoothed_per_million
## 1 0 0
## 2 0 0
## 3 0 0
## 4 0 0
## 5 0 0
## 6 0 0Cleaning the Data
In this data 0 cases in a country will be changed to 0, in order to make the country uncolored
agregat_MPX <- MPX %>% group_by(name) %>% summarise(total_kasus = sum(new_cases))
agregat_MPX$name <- recode(agregat_MPX$name, "United States" = "United States of America",
"United Kingdom" = "U.K. of Great Britain and Northern Ireland")
agregat_MPX[is.na(agregat_MPX)] = 0
agregat_MPX[agregat_MPX==0] <- NA
agregat_MPX <- agregat_MPX %>% arrange(desc(total_kasus))
peta_mpx <- merge(dunia, agregat_MPX)
names(peta_mpx)## [1] "name" "status" "color_code" "region" "iso3"
## [6] "continent" "iso_3166_1_" "french_shor" "total_kasus"Data Check
Checking the agregate files from cleaning
head(agregat_MPX)## # A tibble: 6 × 2
## name total_kasus
## <chr> <dbl>
## 1 World 38709
## 2 United States of America 12636
## 3 Spain 5792
## 4 Germany 3213
## 5 U.K. of Great Britain and Northern Ireland 3201
## 6 Brazil 3184Combine shp file in two methods
So do you can measure, which country was the biggest number in monkeypox infection?
mycolours <- brewer.pal(5, "OrRd")
spplot(peta_mpx,"total_kasus", par.settings = list(axis.line = list(col ="transparent")),
main = "Case MPX", cuts = 4, col ="gray", col.regions = mycolours)
Epicurve of World Cases
In this section, I made the Epidemiology Curve of monkeypox cases, based on daily cases and 7DMA cases.
world_time <- filter(MPX, name == "World")
world_time$date <- as.Date(world_time$date)
x2 <- ggplot(world_time) +
geom_bar(aes(x=date, y=new_cases),stat="identity", fill="mediumpurple", alpha = 0.3)+
geom_line(aes(x=date, y=new_cases_smoothed),stat="identity",color="red", size = 1.5)+
labs(title= "Total dan Tren Global Monkeypox",
x="Tanggal",y="Total Kasus") + scale_x_date(date_breaks = "14 day", date_labels = "%d %b") + theme_minimal() + theme(plot.title = element_text(size=15, hjust=0.5))
x2 + geom_text(x=as.Date("2022-08-11"), y=1100, label="7DMA Cases", color = "red", size = 4) +
geom_text(x=as.Date("2022-08-11"), y=2000, label="Daily Cases", color = "mediumpurple", size = 4) 
Method 2, using mapsf
setwd("C:/Users/dhihr/Downloads/template paskhas/world-administrative-boundaries")
mtq <- st_read("world-administrative-boundaries.shp")## Reading layer `world-administrative-boundaries' from data source
## `C:\Users\dhihr\Downloads\template paskhas\world-administrative-boundaries\world-administrative-boundaries.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 256 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -58.49861 xmax: 180 ymax: 83.6236
## Geodetic CRS: WGS 84mtq <- left_join(mtq, agregat_MPX)## Joining, by = "name"mtq$total_kasus[is.na(mtq$total_kasus)] = 0
mf_map(x = mtq)
head(mtq)## Simple feature collection with 6 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -25.36056 ymin: -18.07492 xmax: 115.3603 ymax: 58.08326
## Geodetic CRS: WGS 84
## status color_code region iso3 continent
## 1 Member State SMR Southern Europe SMR Europe
## 2 Member State SYR Western Asia SYR Asia
## 3 Member State LVA Northern Europe LVA Europe
## 4 Member State CPV Western Africa CPV Africa
## 5 Member State ZMB Eastern Africa ZMB Africa
## 6 Member State BRN South-Eastern Asia BRN Asia
## name iso_3166_1_ french_shor total_kasus
## 1 San Marino SM Saint-Marin 0
## 2 Syrian Arab Republic SY République arabe syrienne 0
## 3 Latvia LV Lettonie 3
## 4 Cape Verde CV Cabo Verde 0
## 5 Zambia ZM Zambie 0
## 6 Brunei Darussalam BN Brunéi Darussalam 0
## geometry
## 1 MULTIPOLYGON (((12.40913 43...
## 2 MULTIPOLYGON (((42.35562 37...
## 3 MULTIPOLYGON (((27.37206 57...
## 4 MULTIPOLYGON (((-24.36556 1...
## 5 MULTIPOLYGON (((32.9404 -9....
## 6 MULTIPOLYGON (((115.0291 4....Method 3, using ggplot
ggplot(data = mtq) + geom_sf()
ggplot(data = mtq) + geom_sf(color = "black", fill = "lightgreen")
x <- ggplot(data = mtq) + geom_sf(aes(fill = total_kasus), color = "snow") + scale_fill_gradient(name = "Total Kasus", low = "deepskyblue", high = "royalblue", na.value="gray") +
theme_minimal() + theme(axis.text.x=element_blank())
ggplotly(x)