library(composits)
library(ggplot2)
library(forecast)
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
library(tidyr)
library(stringr)
library(broom)
library(maptools)
#> Loading required package: sp
#> Checking rgeos availability: FALSE
#> Please note that 'maptools' will be retired by the end of 2023,
#> plan transition at your earliest convenience;
#> some functionality will be moved to 'sp'.
#> Note: when rgeos is not available, polygon geometry computations in maptools depend on gpclib,
#> which has a restricted licence. It is disabled by default;
#> to enable gpclib, type gpclibPermit()
As describe in the paper four dimension reduction methods are used in
the multivariate ensemble PCA, DOBIN, ICS and ICA. It is recomended that
users verified the scaling and centering options of these methods for
each particular example.
Univariate time series outliers
This example uses a univariate time series containing the daily gold
prices from the R package forecast.
gold2 <- forecast::na.interp(gold)
out <- uv_tsout_ens(gold2)
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
inds <- names(which(table(out$outliers) > 2))
ts_gold <- dplyr::as_tibble(gold2) %>% mutate(t = 1:length(gold2)) %>% rename(value = x)
ggplot(ts_gold, aes(x=t, y=value)) +
geom_line() +
geom_vline(xintercept =as.numeric(inds), color="red", alpha=0.8, size=0.5, linetype ='dashed') + ylab("Gold prices") +
theme_bw()
#> Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.
The red dashed vertical lines show the time points that have been
picked by 3 or more outlier detection methods.
Next we look at the quarterly production of woollen yarn in
Australia. This time series is also taken from the R package
forecast.
out <- uv_tsout_ens(woolyrnq)
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 59.5 days
inds <- names(which(table(out$outliers) > 2))
ts_wool <- dplyr::as_tibble(woolyrnq) %>% mutate(t = 1:length(woolyrnq)) %>% rename(value = x)
ggplot(ts_wool, aes(x=t, y=value)) +
geom_line() +
geom_vline(xintercept =as.numeric(inds), color="red", alpha=0.8, size=0.5, linetype ='dashed') + ylab("Woollen Yarn Production") +
theme_bw()
#> Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.
Multivariate time series outliers
This example includes EU stock market data from the package
datasets. It contains the daily closing prices of 4 European
stock indices: Germany DAX, Switzerland SMI, France CAC, and UK FTSE.
First we plot the data.
stpart <- EuStockMarkets[1:600, ]
stpart <- EuStockMarkets[1:600, ]
as_tibble(stpart) %>% mutate(t = 1:n()) %>%
pivot_longer(cols=1:4) %>%
ggplot2::ggplot( ggplot2::aes(x = t, y = value, color = name)) + ggplot2::geom_line() + ggplot2::theme_bw()
Then we find multivariate outliers. To find the multivariate
outliers, first we decompose the time series to univariate outliers by
using Principle Component Analysis (PCA), Independent Component Analysis
(ICA), DOBIN (Distance based Outlier BasIs using Neighbours) and ICS
(Invariant Coordinate Selection) decomposition methods. The
fast=TRUE option leaves out ICS decomposition. Then for a
selected number of components (default being 2) we find outliers using
the univariate time series ensemble. The outliers are given in the table
below.
out <- mv_tsout_ens(stpart, fast=TRUE)
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
out$outliers
#> Indices Total_Score Num_Coords Num_Methods DOBIN PCA ICA
#> 11 36 1.6928 3 3 0.5460645 0.6552774 0.4914581
#> 18 127 0.3120 1 1 0.3120000 0.0000000 0.0000000
#> 20 205 0.8160 3 2 0.2733508 0.4059739 0.1366754
#> 32 316 1.4400 3 2 0.4710592 0.6156465 0.3532944
#> 34 319 0.6880 2 2 0.1720000 0.0000000 0.5160000
#> 41 331 1.5008 3 3 0.2728727 0.4093091 0.8186182
#> 42 366 0.3120 1 1 0.0000000 0.3120000 0.0000000
#> 48 546 0.3120 1 1 0.3120000 0.0000000 0.0000000
#> forecast tsoutliers otsad anomalize Gap_Score_2
#> 11 0 1.248 0.320 0.1248 72
#> 18 0 0.312 0.000 0.0000 5
#> 20 0 0.624 0.192 0.0000 30
#> 32 0 1.248 0.192 0.0000 60
#> 34 0 0.624 0.064 0.0000 23
#> 41 0 1.248 0.128 0.1248 63
#> 42 0 0.312 0.000 0.0000 5
#> 48 0 0.312 0.000 0.0000 5
draw_table_html(out)
|
|
Indices
|
DOBIN
|
PCA
|
ICA
|
Num_Coords
|
forecast
|
tsoutliers
|
otsad
|
anomalize
|
Num_Methods
|
Gap_Score_2
|
Total_Score
|
|
11
|
36
|
0.55
|
0.66
|
0.49
|
3
|
0
|
1.25
|
0.32
|
0.12
|
3
|
72
|
1.69
|
|
41
|
331
|
0.27
|
0.41
|
0.82
|
3
|
0
|
1.25
|
0.13
|
0.12
|
3
|
63
|
1.5
|
|
32
|
316
|
0.47
|
0.62
|
0.35
|
3
|
0
|
1.25
|
0.19
|
0
|
2
|
60
|
1.44
|
|
20
|
205
|
0.27
|
0.41
|
0.14
|
3
|
0
|
0.62
|
0.19
|
0
|
2
|
30
|
0.82
|
|
34
|
319
|
0.17
|
0
|
0.52
|
2
|
0
|
0.62
|
0.06
|
0
|
2
|
23
|
0.69
|
|
18
|
127
|
0.31
|
0
|
0
|
1
|
0
|
0.31
|
0
|
0
|
1
|
5
|
0.31
|
|
42
|
366
|
0
|
0.31
|
0
|
1
|
0
|
0.31
|
0
|
0
|
1
|
5
|
0.31
|
|
48
|
546
|
0.31
|
0
|
0
|
1
|
0
|
0.31
|
0
|
0
|
1
|
5
|
0.31
|
The decomposed time series is shown in the figure below.
plot_decomposed_all(obj=out, X = stpart)
animate_ts_ensemble(out, X= stpart, max_frames = 100)
The decomposition above shows three decomposition methods. If you
want to plot only one method, you can use the
plot_decomposed function. The outlying time-points are
shown in the figure above and below.
plot_decomposed(obj=out, X = stpart, method = "pca")
Composite time series outliers
For this example we use daily mortality counts in Spain organised by
autonomous communities provided by the Spanish government. The dataset
records mortality counts from the 18th of April 2018 until the 31th of
July 2020. To make the data compositional we divide the autonomous
communities counts by the total count for each day.
Compositional data live on a simplex, because the sums add up to a
constant. Before analysing such data, it is important to make the data
unconstrained. We make the data unconstrained by using a coordinate
transformation called the null space coordinate transformation.
We use the multivariate outlier detection ensemble discussed above on
these unconstrained data. More details are available in our paper.
First let us load this data and plot it.
data('spanish_morte')
df <- spanish_morte[[1]]
uniq_dates <- spanish_morte[[2]]
df2 <- cbind.data.frame(uniq_dates, df)
as_tibble(df2) %>%
pivot_longer(cols=2:20) %>%
ggplot2::ggplot( ggplot2::aes(x = uniq_dates, y = value, color = name)) + ggplot2::geom_line() + ggplot2::theme_bw()
Next we find the outliers using compositional time series outlier
ensemble.
out <- comp_tsout_ens(df, fast=FALSE)
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
#> Warning in locate.outliers.oloop(y = y, fit = fit, types = types, cval = cval, :
#> stopped when 'maxit.oloop = 4' was reached
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
#> Converting from tbl_df to tbl_time.
#> Auto-index message: index = date
#> frequency = 7 days
#> trend = 91 days
out$outliers
#> Indices Total_Score Num_Coords Num_Methods DOBIN PCA ICS
#> 9 608 1.333810 3 3 0 0.2659128 0.6785748
#> 14 614 1.245822 3 2 0 0.4823768 0.4103209
#> 17 617 2.527133 3 4 0 0.5900087 1.0411267
#> 18 618 2.058266 3 4 0 0.7032293 0.7282780
#> 19 619 2.011199 3 3 0 0.6871481 0.7116240
#> 20 620 2.320759 3 3 0 0.9789607 0.7211629
#> 21 621 1.718309 3 3 0 0.5712493 0.7288769
#> 22 622 1.905486 3 3 0 0.7245766 0.7503857
#> 23 623 1.490212 3 2 0 0.5770036 0.4908127
#> 24 624 1.490212 3 2 0 0.5770036 0.4908127
#> 25 625 1.490212 3 2 0 0.5770036 0.4908127
#> 26 626 1.490212 3 2 0 0.5770036 0.4908127
#> 27 627 1.180651 3 2 0 0.5472661 0.2327586
#> 28 628 1.180651 3 2 0 0.5472661 0.2327586
#> 30 630 1.490212 3 2 0 0.5770036 0.4908127
#> 31 631 1.490212 3 2 0 0.5770036 0.4908127
#> 38 638 1.105713 3 3 0 0.4284070 0.2488990
#> ICA forecast tsoutliers otsad anomalize Gap_Score_2
#> 9 0.3893229 0.0000000 0.2751652 0.6842925 0.3743527 1
#> 14 0.3531239 0.0000000 0.0000000 0.6842925 0.5615291 0
#> 17 0.8959980 0.3095609 0.8254957 0.4561950 0.9358818 4
#> 18 0.6267591 0.6191218 0.2751652 0.2280975 0.9358818 3
#> 19 0.6124265 0.6191218 0.0000000 0.4561950 0.9358818 2
#> 20 0.6206358 0.9286827 0.0000000 0.4561950 0.9358818 3
#> 21 0.4181830 0.9286827 0.0000000 0.2280975 0.5615291 2
#> 22 0.4305234 0.9286827 0.0000000 0.2280975 0.7487054 2
#> 23 0.4223955 0.9286827 0.0000000 0.0000000 0.5615291 1
#> 24 0.4223955 0.9286827 0.0000000 0.0000000 0.5615291 1
#> 25 0.4223955 0.9286827 0.0000000 0.0000000 0.5615291 1
#> 26 0.4223955 0.9286827 0.0000000 0.0000000 0.5615291 1
#> 27 0.4006261 0.6191218 0.0000000 0.0000000 0.5615291 0
#> 28 0.4006261 0.6191218 0.0000000 0.0000000 0.5615291 0
#> 30 0.4223955 0.9286827 0.0000000 0.0000000 0.5615291 1
#> 31 0.4223955 0.9286827 0.0000000 0.0000000 0.5615291 1
#> 38 0.4284070 0.0000000 0.2751652 0.4561950 0.3743527 0
draw_table_html(out, uniq_dates)
|
|
Indices
|
DOBIN
|
PCA
|
ICS
|
ICA
|
Num_Coords
|
forecast
|
tsoutliers
|
otsad
|
anomalize
|
Num_Methods
|
Gap_Score_2
|
Total_Score
|
|
17
|
2020-03-19
|
0
|
0.59
|
1.04
|
0.9
|
3
|
0.31
|
0.83
|
0.46
|
0.94
|
4
|
4
|
2.53
|
|
20
|
2020-03-22
|
0
|
0.98
|
0.72
|
0.62
|
3
|
0.93
|
0
|
0.46
|
0.94
|
3
|
3
|
2.32
|
|
18
|
2020-03-20
|
0
|
0.7
|
0.73
|
0.63
|
3
|
0.62
|
0.28
|
0.23
|
0.94
|
4
|
3
|
2.06
|
|
19
|
2020-03-21
|
0
|
0.69
|
0.71
|
0.61
|
3
|
0.62
|
0
|
0.46
|
0.94
|
3
|
2
|
2.01
|
|
22
|
2020-03-24
|
0
|
0.72
|
0.75
|
0.43
|
3
|
0.93
|
0
|
0.23
|
0.75
|
3
|
2
|
1.91
|
|
21
|
2020-03-23
|
0
|
0.57
|
0.73
|
0.42
|
3
|
0.93
|
0
|
0.23
|
0.56
|
3
|
2
|
1.72
|
|
23
|
2020-03-25
|
0
|
0.58
|
0.49
|
0.42
|
3
|
0.93
|
0
|
0
|
0.56
|
2
|
1
|
1.49
|
|
24
|
2020-03-26
|
0
|
0.58
|
0.49
|
0.42
|
3
|
0.93
|
0
|
0
|
0.56
|
2
|
1
|
1.49
|
|
25
|
2020-03-27
|
0
|
0.58
|
0.49
|
0.42
|
3
|
0.93
|
0
|
0
|
0.56
|
2
|
1
|
1.49
|
|
26
|
2020-03-28
|
0
|
0.58
|
0.49
|
0.42
|
3
|
0.93
|
0
|
0
|
0.56
|
2
|
1
|
1.49
|
|
30
|
2020-04-01
|
0
|
0.58
|
0.49
|
0.42
|
3
|
0.93
|
0
|
0
|
0.56
|
2
|
1
|
1.49
|
|
31
|
2020-04-02
|
0
|
0.58
|
0.49
|
0.42
|
3
|
0.93
|
0
|
0
|
0.56
|
2
|
1
|
1.49
|
|
9
|
2020-03-10
|
0
|
0.27
|
0.68
|
0.39
|
3
|
0
|
0.28
|
0.68
|
0.37
|
3
|
1
|
1.33
|
|
14
|
2020-03-16
|
0
|
0.48
|
0.41
|
0.35
|
3
|
0
|
0
|
0.68
|
0.56
|
2
|
0
|
1.25
|
|
27
|
2020-03-29
|
0
|
0.55
|
0.23
|
0.4
|
3
|
0.62
|
0
|
0
|
0.56
|
2
|
0
|
1.18
|
|
28
|
2020-03-30
|
0
|
0.55
|
0.23
|
0.4
|
3
|
0.62
|
0
|
0
|
0.56
|
2
|
0
|
1.18
|
|
38
|
2020-04-09
|
0
|
0.43
|
0.25
|
0.43
|
3
|
0
|
0.28
|
0.46
|
0.37
|
3
|
0
|
1.11
|
The table above shows the outlying dates and the associated scores.
The outlying dates correspond to the COVID-19 outbreak. We can see the
decomposed univariate time series in the following plot and the spike
around mid March in 2020.
plot_decomposed_all(obj=out, X = df)
# animate_ts_ensemble(out, X= df, max_frames = 1)
animate_ts_ensemble(out, X= df, max_frames = 100)
Next we apportion the outlying scores back to the autonomous
communities. This is another coordinate transformation from the outlier
score space to the original coordinate space.
apportioned <- apportion_scores_comp(out)
apportioned1 <- cbind.data.frame(colnames(df), apportioned$scores_out)
colnames(apportioned1)[1] <- 'reg'
colnames(apportioned1)[2:dim(apportioned1)[2]] <- paste(uniq_dates[out$outliers$Indices])
apportioned1
#> reg 2020-03-10 2020-03-16 2020-03-19 2020-03-20 2020-03-21 2020-03-22
#> 1 AN 0.32273901 0.36379383 0.51008233 0.49953318 0.47993907 0.60880806
#> 2 AR 0.05842446 0.03628315 0.08867115 0.06241029 0.06258093 0.06317523
#> 3 AS 0.04788701 0.04917141 0.06750365 0.06540414 0.06301702 0.07871691
#> 4 IB 0.03272778 0.03216793 0.04512506 0.04295105 0.04148237 0.05114173
#> 5 CN 0.14526591 0.12197341 0.21225648 0.18119642 0.17716649 0.20367105
#> 6 CB 0.26052347 0.15993285 0.37453350 0.27308470 0.27268160 0.27533277
#> 7 CL 0.19281874 0.13759115 0.27541814 0.21735882 0.21481594 0.23154029
#> 8 CM 0.22278322 0.18650323 0.29992750 0.26216600 0.25592254 0.29698531
#> 9 CT 0.44674165 0.46908523 0.80444355 0.72269075 0.70008671 0.84875176
#> 10 VC 0.16692874 0.15239173 0.21248857 0.19921218 0.19299631 0.23387763
#> 11 EX 0.09255991 0.06235281 0.13107098 0.10071354 0.09991861 0.10514685
#> 12 GA 0.20325082 0.17883181 0.28854412 0.25500124 0.24836658 0.29200226
#> 13 MD 0.30787102 0.30000515 0.34201508 0.32519694 0.31104019 0.40414800
#> 14 MC 0.14413758 0.10689759 0.20590698 0.16572030 0.16334743 0.17895361
#> 15 NC 0.03592833 0.02410593 0.05383122 0.04077918 0.04049675 0.04235427
#> 16 PV 0.11483314 0.08101374 0.16073198 0.12674209 0.12531959 0.13467530
#> 17 RI 0.07900978 0.04783895 0.11512343 0.08361032 0.08357671 0.08379774
#> 18 CE 0.09986108 0.06115920 0.14469427 0.10574578 0.10561809 0.10645791
#> 19 ML 0.05886064 0.03691069 0.08569440 0.06325983 0.06308599 0.06422932
#> 2020-03-23 2020-03-24 2020-03-25 2020-03-26 2020-03-27 2020-03-28 2020-03-29
#> 1 0.46571844 0.48281480 0.43515850 0.43515850 0.43515850 0.43515850 0.39368428
#> 2 0.06231721 0.06054775 0.04340074 0.04340074 0.04340074 0.04340074 0.01772839
#> 3 0.06562172 0.06597069 0.05881726 0.05881726 0.05881726 0.05881726 0.04958263
#> 4 0.04378383 0.04371913 0.03847824 0.03847824 0.03847824 0.03847824 0.03127780
#> 5 0.17849365 0.18031010 0.14590067 0.14590067 0.14590067 0.14590067 0.10166734
#> 6 0.27601969 0.27492000 0.19130654 0.19130654 0.19130654 0.19130654 0.07633869
#> 7 0.21874553 0.21865557 0.16458211 0.16458211 0.16458211 0.16458211 0.09107462
#> 8 0.27374503 0.27025637 0.22308917 0.22308917 0.22308917 0.22308917 0.15434195
#> 9 0.61709752 0.66343808 0.56110469 0.56110469 0.56110469 0.56110469 0.48511442
#> 10 0.21482720 0.20972006 0.18228610 0.18228610 0.18228610 0.18228610 0.13812317
#> 11 0.10224505 0.10177798 0.07458442 0.07458442 0.07458442 0.07458442 0.03704501
#> 12 0.25602619 0.25682812 0.21391287 0.21391287 0.21391287 0.21391287 0.15675908
#> 13 0.41204478 0.37991149 0.35885653 0.35885653 0.35885653 0.35885653 0.28667247
#> 14 0.16655965 0.16665746 0.12786746 0.12786746 0.12786746 0.12786746 0.07535663
#> 15 0.03956592 0.04010252 0.02883474 0.02883474 0.02883474 0.02883474 0.01427100
#> 16 0.12962985 0.12873521 0.09690604 0.09690604 0.09690604 0.09690604 0.05249515
#> 17 0.08317318 0.08335407 0.05722341 0.05722341 0.05722341 0.05722341 0.02200588
#> 18 0.10565943 0.10571727 0.07315667 0.07315667 0.07315667 0.07315667 0.02904578
#> 19 0.06291979 0.06308748 0.04415138 0.04415138 0.04415138 0.04415138 0.01867997
#> 2020-03-30 2020-04-01 2020-04-02 2020-04-09
#> 1 0.39368428 0.43515850 0.43515850 0.23441794
#> 2 0.01772839 0.04340074 0.04340074 0.02732577
#> 3 0.04958263 0.05881726 0.05881726 0.02664851
#> 4 0.03127780 0.03847824 0.03847824 0.01683874
#> 5 0.10166734 0.14590067 0.14590067 0.07443229
#> 6 0.07633869 0.19130654 0.19130654 0.10336609
#> 7 0.09107462 0.16458211 0.16458211 0.08348470
#> 8 0.15434195 0.22308917 0.22308917 0.09505268
#> 9 0.48511442 0.56110469 0.56110469 0.38183808
#> 10 0.13812317 0.18228610 0.18228610 0.06774689
#> 11 0.03704501 0.07458442 0.07458442 0.03761599
#> 12 0.15675908 0.21391287 0.21391287 0.10150948
#> 13 0.28667247 0.35885653 0.35885653 0.10009500
#> 14 0.07535663 0.12786746 0.12786746 0.06419777
#> 15 0.01427100 0.02883474 0.02883474 0.01673485
#> 16 0.05249515 0.09690604 0.09690604 0.04688627
#> 17 0.02200588 0.05722341 0.05722341 0.03213524
#> 18 0.02904578 0.07315667 0.07315667 0.04032321
#> 19 0.01867997 0.04415138 0.04415138 0.02444588
Next we plot this using the map of Spain.
shapefile_ccaa <- readShapePoly("ComunidadesAutonomas_ETRS89_30N/Comunidades_Autonomas_ETRS89_30N.shp", proj4string=CRS("+proj=longlat"))
#> Warning: readShapePoly is deprecated; use rgdal::readOGR or sf::st_read
data_ccaa <- tidy(shapefile_ccaa)
#> Regions defined for each Polygons
nombres_ccaa <- tibble(shapefile_ccaa$Texto) %>%
mutate(id = as.character(seq(0, nrow(.)-1)))
reg <- c("AN", "AR", "AS", "PM", "CN", "CB", "CL", "CM",
"CT", "VC", "EX", "GA", "MD", "MU", "NA", "PV",
"LO", "CE", "CE")
nombres_ccaa <- tibble(reg) %>%
mutate(id = as.character(seq(0, nrow(.)-1)))
data_ccaa_mapa <- data_ccaa %>%
left_join(nombres_ccaa, by = "id")
df3 <- apportioned1 %>%
mutate(reg = str_replace(reg, "RI", "LO")) %>%
mutate(reg = str_replace(reg, "NC", "NA")) %>%
mutate(reg = str_replace(reg, "MC", "MU")) %>%
pivot_longer(-reg, names_to = "date", values_to = "score")
sp_data <- inner_join(data_ccaa_mapa, df3, by="reg") %>%
filter(reg != "CN")
sp_data %>%
ggplot() +
geom_polygon(aes(x= long, y = lat, group=group, fill=score),
color="black", size=0.2) +
theme_bw() +
scale_fill_distiller(palette = "YlOrBr", direction = 1) +
xlab("") +
ylab("") +
coord_fixed() +
facet_wrap(~date, ncol=5) +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank())
See our website or
our paper (Kandanaarachchi et al. 2020)
for more examples.