-
Notifications
You must be signed in to change notification settings - Fork 22
/
Copy pathV23_map_grid.Rmd
92 lines (73 loc) · 2.23 KB
/
V23_map_grid.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
```{r include=FALSE}
knitr::opts_chunk$set(cache = T, warning = F, message = F,
class.output = "output", out.width='100%',
fig.asp = 0.5, fig.align = 'center')
library(tidyverse)
```
## Mapping data with grid
```{r}
library(sf)
```
### Loading Taiwan map
```{r}
TW.island <- st_read("data/shapefiles/COUNTY_MOI_1090820.shp") %>%
st_transform(3826) %>%
mutate(id = row_number())
```
### Building grid
```{r}
# Defining grid size
grid.extent <-
matrix(c(-50000, 2920000, # (Xmin, Ymax)
610000, 2920000, # (Xmax, Ymax)
610000, 2420000, # (Xmax, Ymin)
-50000, 2420000, # (Xmin, Ymin)
-50000, 2920000), # (Xmin, Ymax)
byrow = TRUE, ncol = 2) %>%
list() %>% # convert to list for st_polygon()
st_polygon() %>% # generate polygon
st_sfc(crs = 3826) # convert format and crs
# plot(grid.extent)
# Generating grid
Grid.sys <-
st_make_grid(grid.extent,
n = c(132, 100), # Resolution of grids
crs = 3826, # crs: TWD97 121
what = 'polygons') %>% # output format: polygon
st_sf('geometry' = ., data.frame('ID' = 1:length(.))) # convert to sf with id
# st_transform(3826) # assigning crs again ?
plot(Grid.sys)
```
```{r}
Grid.TW <-
Grid.sys[subset(TW.island),]
plot(Grid.TW)
```
### loading data
```{r}
president_vote <- readxl::read_xlsx('data/president.xlsx') %>%
mutate(total = chu + tsai + song) %>%
mutate(chu_ratio = chu / total,
tsai_ratio = tsai / total,
song_ratio = song / total,
tsai_chu_ratio = tsai / chu)
```
### Merging data
```{r}
tw_info <- TW.island %>%
st_set_geometry(NULL) %>%
left_join(president_vote, by=c("COUNTYNAME"="county"))
```
```{r}
# TW_info <- sf::st_intersects(Grid.TW, TW.island) # creat a data.frame of IDs in IBA for 1km grid
grid_id <- sapply(st_intersects(Grid.TW, TW.island), function(z) if (length(z)==0) NA_integer_ else z[1])
Grid.TW <- Grid.TW %>%
mutate(grid_id = grid_id) %>%
left_join(tw_info, by=c("grid_id"="id"))
```
```{r}
Grid.TW %>%
ggplot(aes(fill = tsai_ratio)) + geom_sf(lwd = 0.1, color="black") +
scale_fill_continuous(high="#2EFF71", low="blue") +
theme_void()
```