-
-
Notifications
You must be signed in to change notification settings - Fork 1
/
_14-ex-ja.Rmd
173 lines (154 loc) · 7.11 KB
/
_14-ex-ja.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
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
```{asis 14-ex-asis1, message=FALSE}
回答するには、以下のパッケージをアタッチすることとする (他のパッケージも必要に応じてアタッチする)。
```
```{r 14-ex-e0, message=FALSE, warning=FALSE}
library(sf)
library(dplyr)
library(purrr)
library(terra)
library(osmdata)
library(spDataLarge)
```
E1. 100 m セル解像度の住民情報を含む csv ファイルをダウンロードしなさい (https://www.zensus2011.de/SharedDocs/Downloads/DE/Pressemitteilung/DemografischeGrunddaten/csv_Bevoelkerung_100m_Gitter.zip?__blob=publicationFile&v=3)。
解凍したファイルのサイズは 1.23 GB である。
このファイルを R に読み込むには、`readr::read_csv`を使うことができる。
16 GB の RAM を搭載したパソコンで 30 秒かかる。
`data.table::fread()` はさらに速く、`data.table()` クラスのオブジェクトを返す。
`dplyr::as_tibble()` を使用して、それを tibble に変換しなさい。
住民ラスタを作成し、セル解像度 1 km に集約し、クラスの平均値を用いて作成した住民ラスタ (`inh`) との差を比較しなさい。
```{r, 14-ex-e1, eval=FALSE}
# 粗い住民ラスタ (解像度 1 km)
#*******************************************
# 住民ラスタ (解像度は低い); これは
# 以前の演習の結果
data("census_de", package = "spDataLarge")
input = select(census_de, x = x_mp_1km, y = y_mp_1km, pop = Einwohner,
women = Frauen_A, mean_age = Alter_D, hh_size = HHGroesse_D)
input_tidy = dplyr::mutate(input, dplyr::across(.fns = ~ifelse(. %in% c(-1, -9), NA, .)))
input_ras = terra::rast(input_tidy, type = "xyz", crs = "EPSG:3035")
inh_coarse = input_ras$pop
# 再分類 クラスをクラス平均で住民に変換
rcl = matrix(c(1, 1, 125, 2, 2, 375, 3, 3, 1250, 4, 4, 3000, 5, 5, 6000,
6, 6, 8000), ncol = 3, byrow = TRUE)
inh_coarse = terra::classify(inh_coarse, rcl = rcl, right = NA)
# 詳細の住民ラスタ (解像度 100 m)
#******************************************
url =
paste0("https://www.zensus2011.de/SharedDocs/Downloads/DE/Pressemitteilung/",
"DemografischeGrunddaten/csv_Bevoelkerung_100m_Gitter.zip",
"?__blob=publicationFile&v=3")
# 詳細の住民ラスタをダウンロード
download.file(url = url, destfile = file.path(tempdir(), "census.zip"),
method = "auto", mode = "wb")
# ファイル名を表示
nms = unzip(file.path(tempdir(), "census.zip"), list = TRUE)
# csv ファイルのみ解凍
base_name = grep(".csv$", nms$Name, value = TRUE)
unzip(file.path(tempdir(), "census.zip"), files = base_name, exdir = tempdir())
# csv ファイルを読む
input = data.table::fread(file.path(tempdir(), base_name)) |>
dplyr::as_tibble()
input = select(input, x = starts_with("x_mp_1"),
y = starts_with("y_mp_1"), inh = Einwohner)
# -1 から -9 を NA に設定
input = dplyr::mutate(input,
dplyr::across(.fns = ~ifelse(. %in% c(-1, -9), NA, .)))
# テーブルをラスタに変換 (x と y はセルの中央値)
inh_fine = terra::rast(input, type = "xyz", crs = "EPSG:3035")
# 注: inh_fine は、ラスタセルごとの住民数
# なお、粗いラスタでは 1km の時は住民数の平均値であった
# 粗いラスタと詳細なラスタを比較
#******************************************
# 粗いラスタの解像度に集約
inh_fine = terra::aggregate(
inh_fine, fact = terra::res(inh_coarse)[1] / terra::res(inh_fine)[1],
fun = sum, na.rm = TRUE)
# オリジンは同じはず
terra::origin(inh_fine) = terra::origin(inh_coarse)
# 比較する
summary(inh_fine - inh_coarse)
plot(inh_fine - inh_coarse)
plot(abs(inh_fine - inh_coarse) > 1000)
# Berlin などの大都市では偏差が最大
terra::global((abs(inh_fine - inh_coarse) > 1000), fun = "sum", na.rm = TRUE)
# 18,121 のセルで偏差が > 1000 住民
terra::global((abs(inh_fine - inh_coarse) > 5000), fun = "sum", na.rm = TRUE)
# 338 のセルで偏差が > 5000
```
E2. 仮に、自転車店が主に高齢者に電動自転車を販売していたとしよう。
それに応じて年齢ラスタを変更し、残りの分析を繰り返し、その変化を元の結果と比較しなさい。
```{r, 14-ex-e2, eval=FALSE}
# 前の演習で `input_ras` をすでに作成済みと仮定
# 必要なデータをアタッチ
data("metro_names", "shops", package = "spDataLarge")
# 高齢者は電動自転車を選択すると仮定している
# よって、高齢者の多い地域のラスタセルに
# 重み付けをする
rcl_pop = matrix(c(1, 1, 127, 2, 2, 375, 3, 3, 1250,
4, 4, 3000, 5, 5, 6000, 6, 6, 8000),
ncol = 3, byrow = TRUE)
rcl_women = matrix(c(1, 1, 3, 2, 2, 2, 3, 3, 1, 4, 5, 0),
ncol = 3, byrow = TRUE)
# 最高齢の人がいるセルクラス (3 から 5)
# に、最大の重み
rcl_age = matrix(c(1, 1, 1, 2, 2, 1, 3, 5, 3),
ncol = 3, byrow = TRUE)
rcl_hh = rcl_women
rcl = list(rcl_pop, rcl_women, rcl_age, rcl_hh)
reclass = input_ras
for (i in 1:terra::nlyr(reclass)) {
reclass[[i]] = terra::classify(x = reclass[[i]], rcl = rcl[[i]], right = NA)
}
names(reclass) = names(input_ras)
# ここからの解析は本の通り
# sf オブジェクト metros に駅名を付与
#************************************
metro_names = dplyr::pull(metro_names, city) |>
as.character() |>
{\(x) ifelse(x == "Velbert", "Düsseldorf", x)}() |>
{\(x) gsub("ü", "ue", x)}()
pop_agg = terra::aggregate(reclass$pop, fact = 20, fun = sum, na.rm = TRUE)
pop_agg = pop_agg[pop_agg > 500000, drop = FALSE]
polys = pop_agg |>
terra::patches(directions = 8) |>
terra::as.polygons() |>
sf::st_as_sf()
metros = polys |>
dplyr::group_by(patches) |>
dplyr::summarize()
metros$metro_names = metro_names
# shop/poi 密度ラスタを作成
#*******************************
shops = sf::st_transform(shops, sf::st_crs(reclass))
# poi ラスタを作成
poi = terra::rasterize(x = shops, y = reclass, field = "osm_id", fun = "length")
# 再分類行列を作成
int = classInt::classIntervals(values(poi), n = 4, style = "fisher")
int = round(int$brks)
rcl_poi = matrix(c(int[1], rep(int[-c(1, length(int))], each = 2),
int[length(int)] + 1), ncol = 2, byrow = TRUE)
rcl_poi = cbind(rcl_poi, 0:3)
# 再分類
poi = terra::classify(poi, rcl = rcl_poi, right = NA)
names(poi) = "poi"
# 人口ラスタを削除、poi ラスタを追加
reclass = reclass[[names(reclass) != "pop"]] |>
c(poi)
# 適している位置を探索
#****************************
# 合計点を計算
result = sum(reclass)
# Berlin における自転車店出店適地を見る
berlin = metros[metro_names == "Berlin", ]
berlin_raster = terra::crop(result, berlin)
# summary(berlin_raster)
# berlin_raster
berlin_raster = berlin_raster > 9
berlin_raster[berlin_raster == 0] = NA
# プロットする
leaflet::leaflet() |>
leaflet::addTiles() |>
leaflet::addRasterImage(raster::raster(berlin_raster), colors = "darkgreen", opacity = 0.8) |>
leaflet::addLegend("bottomright", colors = c("darkgreen"),
labels = c("potential locations"), title = "Legend")
```