_09-ex-ja.Rmd

```{r 09-ex-e0, message=FALSE}
library(sf)
library(terra)
library(dplyr)
library(spData)
```

ここでの演習は、新しくオブジェクト `africa` を使用する。
これは、**spData** のデータセット `world` と `worldbank_df` から、以下のように作成する。

```{r 08-mapping-41, warning=FALSE, include=TRUE}
library(spData)
africa = world |> 
  filter(continent == "Africa", !is.na(iso_a2)) |> 
  left_join(worldbank_df, by = "iso_a2") |> 
  select(name, subregion, gdpPercap, HDI, pop_growth) |> 
  st_transform("ESRI:102022") |> 
  st_make_valid() |> 
  st_collection_extract("POLYGON")
```

**spDataLarge** のデータセット `zion` と `nlcd` も使用する。

```{r 08-mapping-42, results='hide', include=TRUE}
zion = read_sf((system.file("vector/zion.gpkg", package = "spDataLarge")))
nlcd = rast(system.file("raster/nlcd.tif", package = "spDataLarge"))
```

E1. **graphics** (ヒント: `plot()`) と **tmap** パッケージ (ヒント: `tm_shape(africa) + ...`) を使って、Africa 全土の人間開発指数 (`HDI`) の地理的分布を示す地図を作成しなさい。

- それぞれの長所を経験に基づいて 2 つ挙げなさい。
- 他の地図作成パッケージを 3 つ挙げ、それぞれの利点を挙げなさい。
- ボーナス: これら 3 つの他のパッケージを使って、さらに 3 つのアフリカの地図を作りなさい。

```{r}
# graphics
plot(africa["HDI"])
# tmap
remotes::install_github("r-tmap/tmap")
library(tmap)
tm_shape(africa) + 
  tm_polygons("HDI")
# ggplot
library(ggplot2)
ggplot() +
  geom_sf(data = africa, aes(fill = HDI))
# ggplotly
library(plotly)
g = ggplot() +
  geom_sf(data = africa, aes(fill = HDI))
ggplotly(g)
# mapsf
library(mapsf)
mf_map(x = africa, var = "HDI", type = "choro")
```

E2. 前の演習で作成した **tmap** を拡張して、凡例に 3 つのビンを設定しなさい: "High" (0.7 を超える `HDI`)、"Medium" (0.55 と 0.7 の間の `HDI`)、"Low" (0.55 を下回る `HDI`)。

- ボーナス: 例えば、凡例のタイトル、クラスラベル、色パレットを変更することで、マップの美観を改善しなさい。

```{r}
library(tmap)
tm_shape(africa) + 
  tm_polygons("HDI",
              fill.scale = tm_scale_intervals(breaks = c(0, 0.55, 0.7, 1),
                                              labels = c("Low", "Medium", "High"),
                                              values = "-viridis"),
              fill.legend = tm_legend(title = "Human Development Index")) 
```

E3. `africa` の小地域を地図上に表示しなさい。
デフォルトの色パレットと凡例のタイトルを変更しなさい。
次に、この地図と前の練習で作成した地図を組み合わせて、一つのプロットし統合しなさい。

```{r}
asubregions = tm_shape(africa) +
  tm_polygons("subregion",
              fill.scale = tm_scale_categorical(values = "Set3"),
              fill.legend = tm_legend(title = "Subregion:"))
ahdi = tm_shape(africa) + 
  tm_polygons("HDI",
              fill.scale = tm_scale_intervals(breaks = c(0, 0.55, 0.7, 1),
                                              labels = c("Low", "Medium", "High"),
                                              values = "-viridis"),
              fill.legend = tm_legend(title = "Human Development Index:")) 
tmap_arrange(ahdi, asubregions)
```

E4. Zion 国立公園の土地被覆マップを作成しなさい。

- 土地被覆カテゴリの認識に合わせてデフォルトの色を変更
- 縮尺バーと北矢印を追加し、両方の位置を変更して地図の美観を向上
- ボーナス: ザイオン国立公園のユタ州との位置関係を示す挿入地図を追加 (ヒント: ユタを表すオブジェクトは `us_states` データセットから抽出できる)。

```{r}
tm_shape(nlcd) +
  tm_raster(col.scale = tm_scale_categorical(values = c("#495EA1", "#AF5F63", "#EDE9E4",
                                                        "#487F3F", "#EECFA8", "#A4D378",
                                                        "#FFDB5C", "#72D593"), levels.drop = TRUE)) +
  tm_scalebar(bg.color = "white", position = c("left", "bottom")) +
  tm_compass(bg.color = "white", position = c("right", "top")) +
  tm_layout(legend.position = c("left", "top"), legend.bg.color = "white")
```

```{r}
# Bonus
utah = subset(us_states, NAME == "Utah")
utah = st_transform(utah, st_crs(zion))

zion_region = st_bbox(zion) |> 
  st_as_sfc()

main = tm_shape(nlcd) +
  tm_raster(col.scale = tm_scale_categorical(values = c("#495EA1", "#AF5F63", "#EDE9E4",
                                                        "#487F3F", "#EECFA8", "#A4D378",
                                                        "#FFDB5C", "#72D593"), levels.drop = TRUE)) +
  tm_scalebar(bg.color = "white", position = c("left", "bottom")) +
  tm_compass(bg.color = "white", position = c("right", "top")) +
  tm_layout(legend.position = c("left", "top"), legend.bg.color = "white")

inset = tm_shape(utah) +
  tm_polygons() +
  tm_text("UTAH", size = 3) +
  #tm_shape(zion) +
  #tm_polygons(col = "red") +
  tm_shape(zion_region) +
  tm_borders(col = "red") +
  tm_layout(frame = FALSE)

library(grid)
norm_dim = function(obj){
    bbox = st_bbox(obj)
    width = bbox[["xmax"]] - bbox[["xmin"]]
    height = bbox[["ymax"]] - bbox[["ymin"]]
    w = width / max(width, height)
    h = height / max(width, height)
    return(unit(c(w, h), "snpc"))
}
main_dim = norm_dim(zion)
ins_dim = norm_dim(utah)

main_vp = viewport(width = main_dim[1], height = main_dim[2])
ins_vp = viewport(width = ins_dim[1] * 0.4, height = ins_dim[2] * 0.4,
                  x = unit(1, "npc") - unit(0.5, "cm"), y = unit(0.5, "cm"),
                  just = c("right", "bottom"))

grid.newpage()
print(main, vp = main_vp)
pushViewport(main_vp)
print(inset, vp = ins_vp)
```

E5. Eastern Africa の国々のファセットマップを作成しなさい。

- 1 つのファセットは HDI を表し、もう 1 つのファセットは人口増加を表す (ヒント: それぞれ変数`HDI`と`pop_growth`を使用)
- 国ごとに「小さな倍数」を設定

```{r}
ea = subset(africa, subregion == "Eastern Africa")
#1
tm_shape(ea) +
  tm_polygons(c("HDI", "pop_growth"))
#2
tm_shape(ea) +
  tm_polygons() +
  tm_facets_wrap("name")
```

E6. これまでのファセット地図の例に基づいて、East Africa の地図アニメーションを作成しなさい。

- 各国を順番に表示
- HDI を示す凡例とともに各国を順番に表示

```{r, eval=FALSE}
tma1 = tm_shape(ea) +
  tm_polygons() +
  tm_facets(by = "name", nrow = 1, ncol = 1)
tmap_animation(tma1, filename = "tma2.gif", width = 1000, height = 1000)
browseURL("tma1.gif")

tma2 = tm_shape(africa) +
  tm_polygons(fill = "lightgrey") +
  tm_shape(ea) +
  tm_polygons(fill = "darkgrey") +
  tm_shape(ea) +
  tm_polygons(fill = "HDI") +
  tm_facets(by = "name", nrow = 1, ncol = 1)
tmap_animation(tma2, filename = "tma2.gif", width = 1000, height = 1000)
browseURL("tma2.gif")
```

E7. Africa における HDI のインタラクティブ地図を作成しなさい。

- **tmap**
- **mapview**
- **leaflet**
- ボーナス: 各アプローチについて、凡例 (自動的に提供されない場合) とスケールバーを追加しなさい。

```{r, eval=FALSE}
# tmap
tmap_mode("view")
tm_shape(africa) + tm_polygons("HDI") + tm_scalebar()
# mapview
mapview::mapview(africa["HDI"])
# leaflet
africa4326 = st_transform(africa, "EPSG:4326")
library(leaflet)
pal = colorNumeric(palette = "YlGnBu", domain = africa4326$HDI)
leaflet(africa4326) |> 
  addTiles() |> 
  addPolygons(stroke = FALSE, smoothFactor = 0.2, fillOpacity = 1, color = ~pal(HDI)) |> 
  addLegend("bottomright", pal = pal, values = ~HDI, opacity = 1) |> 
  addScaleBar()
```

E8. 交通政策や土地利用政策をよりエビデンスに基づいたものにするために使用できるウェブ地図アプリのアイデアを紙にスケッチしなさい。

  - あなたが住んでいる都市で、1 日あたり数人のユーザー向け
  - あなたが住んでいる国で、1 日あたり数十人のユーザー向け
  - 世界中、1 日あたり数百人のユーザーと大規模なデータ配信が必要な場合

```{asis}
アイデアとしては、現在多くの人が車で短距離を移動しているルートの特定、公園へのアクセスを促す方法、長距離移動を減らすための新規開発の優先順位付けなどが考えられる。

都市レベルでは、ウェブ地図で十分である。

国レベルでは、例えば、shiny 地図アプリケーションが必要だろう。

世界レベルでは、データを提供するデータベースが必要になるだろう。そして、様々なフロントエンドがこれに接続できる。
```

E9. `coffeeApp/app.R` のコードを更新し、Brazil を中心に表示するのではなく、ユーザーがどの国を中心に表示するかを選択しなさい。

- `textInput()` を使いなさい
- `selectInput()` を使いなさい

```{asis}
The answer can be found in the `shinymod` branch of the geocompr repo: https://github.com/Robinlovelace/geocompr/pull/318/files
You create the new widget and then use it to set the center.
Note: the input data must be fed into the map earlier to prevent the polygons disappearing when you change the center this way.
```

E10. **ggplot2** パッケージを使用して、Figure 9.1 と Figure 9.7 をできるだけ忠実に再現しなさい。

```{r}
library(ggplot2)
ggplot() + 
  geom_sf(data = nz, color = NA) +
  coord_sf(crs = st_crs(nz), datum = NA) +
  theme_void()
ggplot() +
  geom_sf(data = nz, fill = NA) +
  coord_sf(crs = st_crs(nz), datum = NA) +
  theme_void()
ggplot() +
  geom_sf(data = nz) + 
  coord_sf(crs = st_crs(nz), datum = NA) +
  theme_void()
# fig 9.7
ggplot() +
  geom_sf(data = nz, aes(fill = Median_income)) +
  coord_sf(crs = st_crs(nz), datum = NA) +
  scale_fill_distiller(palette = "Blues", direction = 1) + 
  theme_void()
ggplot() + 
  geom_sf(data = nz, aes(fill = Island)) + 
  coord_sf(crs = st_crs(nz), datum = NA) +
  scale_fill_manual(values = c("#CC6677", "#332288")) + 
  theme_void()
```

E11. `us_states` と `us_states_df` を結合し、新しいデータセットを使って各州の貧困率を計算しなさい。
次に、総人口に基づいて連続的な範囲カートグラムを作成しなさい。
最後に、貧困率の 2 つの地図を作成し、比較しなさい：(1) 標準的なコロプレス地図と、(2) 作成したカートグラムの境界線を使った地図。
1 枚目と 2 枚目の地図から得られる情報は何か?
両者はどう違うのか?

```{r}
tmap_mode("plot")
library(cartogram)
# prepare the data
us = st_transform(us_states, "EPSG:9311")
us = left_join(us, us_states_df, by = c("NAME" = "state"))
# calculate a poverty rate
us$poverty_rate = us$poverty_level_15 / us$total_pop_15
# create a regular map
ecm1 = tm_shape(us) +
  tm_polygons("poverty_rate", fill.legend = tm_legend(title = "Poverty rate"))
# create a cartogram
us_carto = cartogram_cont(us, "total_pop_15")
ecm2 = tm_shape(us_carto) + 
  tm_polygons("poverty_rate", fill.legend = tm_legend(title = "Poverty rate"))
# combine two maps
tmap_arrange(ecm1, ecm2)
```

E12. Africa の人口増加を視覚化しなさい。
次に、**geogrid** パッケージを使って作成した六角形と正方形のグリッドの地図と比較しなさい。


```{r}
library(geogrid)

hex_cells = calculate_grid(africa, grid_type = "hexagonal", seed = 25, learning_rate = 0.03)
africa_hex = assign_polygons(africa, hex_cells)

reg_cells = calculate_grid(africa, grid_type = "regular", seed = 25, learning_rate = 0.03)
africa_reg = assign_polygons(africa, reg_cells)

tgg1 = tm_shape(africa) +
  tm_polygons("pop_growth", fill.legend = tm_legend(title = "Population's growth (annual %)"))
tgg2 = tm_shape(africa_hex) + 
  tm_polygons("pop_growth", fill.legend = tm_legend(title = "Population's growth (annual %)"))
tgg3 = tm_shape(africa_reg) + 
  tm_polygons("pop_growth", fill.legend = tm_legend(title = "Population's growth (annual %)"))

tmap_arrange(tgg1, tgg2, tgg3)
```