From da90ce9d9273fe71880b42aa787f5ca8d7be2c29 Mon Sep 17 00:00:00 2001 From: babayoshihiko Date: Mon, 26 Feb 2024 06:19:35 +0900 Subject: [PATCH] ch11-16 --- 11-algorithms-ja.Rmd | 89 ++++++++-------- 12-spatial-cv-ja.Rmd | 175 +++++++++++++++--------------- 13-transport-ja.Rmd | 145 +++++++++++++------------ 14-location-ja.Rmd | 11 +- 15-eco-ja.Rmd | 145 ++++++++++++------------- 16-synthesis-ja.Rmd | 16 +-- _11-ex-ja.Rmd | 103 ++++++++++++++++++ _12-ex-ja.Rmd | 237 +++++++++++++++++++++++++++++++++++++++++ _13-ex-ja.Rmd | 72 +++++++++++++ _14-ex-ja.Rmd | 173 ++++++++++++++++++++++++++++++ _15-ex-ja.Rmd | 247 +++++++++++++++++++++++++++++++++++++++++++ 11 files changed, 1116 insertions(+), 297 deletions(-) create mode 100644 _11-ex-ja.Rmd create mode 100644 _12-ex-ja.Rmd create mode 100644 _13-ex-ja.Rmd create mode 100644 _14-ex-ja.Rmd create mode 100644 _15-ex-ja.Rmd diff --git a/11-algorithms-ja.Rmd b/11-algorithms-ja.Rmd index 97881c5..ea11f3a 100644 --- a/11-algorithms-ja.Rmd +++ b/11-algorithms-ja.Rmd @@ -24,7 +24,7 @@ Chapter \@ref(intro) は、ジオコンピューティングは既存のツー しかし、プログラミングを学ぶ強い理由がある。 この章では、プログラミングそのものを教えるわけではない。プログラミングについては、@wickham_advanced_2019、@gillespie_efficient_2016、@xiao_gis_2016 を推奨する。これらの書籍は R や他の言語について教えてくれる。また、地理データに焦点を当て、プログラミング能力を伸ばすための基礎を作ることができる。 -本章は、再現性\index{reproducibility}の重要性について、例を示しながら強調していきたい。 +本章は、再現性\index{さいげんせい@再現性}の重要性について、例を示しながら強調していきたい。 再現性の利点は、他の人があなたの研究を複製することを可能にするだけではない。 再現性のあるコードは、一度だけ実行されるように書かれたコードよりも、計算効率、スケーラビリティ (より大きなデータセットに対して実行するコードの能力)、適応やメンテナンスのしやすさなど、あらゆる面で優れていることが多いのである。 @@ -32,8 +32,8 @@ Chapter \@ref(intro) は、ジオコンピューティングは既存のツー アルゴリズムは、Section \@ref(geometric-algorithms) で説明されているように、一連のステップを使用して入力を変更し、その結果、出力を得るためのレシピである。 共有と再現を容易にするために、アルゴリズムを関数に配置することができる。 それが、Section \@ref(functions) のトピックである。 -ポリゴンの重心\index{centroid}を求める例で、これらの概念を結びつけていこう。 -Chapter \@ref(geometry-operations) で、重心の関数\index{centroid} `st_centroid()` をすでに紹介したが、この例は、一見単純な操作が比較的複雑なコードの結果であることを強調し、次の観察を保証する [@wise_gis_2001] 。 +ポリゴンの重心\index{じゅうしん＠重心}を求める例で、これらの概念を結びつけていこう。 +Chapter \@ref(geometry-operations) で、重心の関数\index{じゅうしん＠重心} `st_centroid()` をすでに紹介したが、この例は、一見単純な操作が比較的複雑なコードの結果であることを強調し、次の観察を保証する [@wise_gis_2001] 。 > 空間データの問題で最も興味深いのは、人間にとっては些細なことに見えることが、コンピュータにとっては驚くほど難しいということである。 @@ -42,7 +42,7 @@ Chapter \@ref(geometry-operations) で、重心の関数\index{centroid} `st_cen ## スクリプト {#scripts} パッケージで配布される関数が R コードの構成要素だとすれば、スクリプトはそれらを論理的な順序でまとめる接着剤となる。 -スクリプトとは、再現可能なワークフローを作り出す目的で、手動または **targets** などの自動化ツールで保存・実行される [@landau_targets_2021]。 +スクリプトとは、再現可能なワークフローを作り出す目的で、手動または **targets** などの自動化ツールで保存・実行される [@landau_targets_2021]。\index{さいげんせい@再現性} プログラミングの初心者にとってスクリプトは敷居が高く聞こえるかもしれないが、単なるプレーンテキストファイルである。 スクリプトは、通常はその言語を表す拡張子で保存される。例えば、Python は `.py`、Rust は `.rs` である。 R スクリプトは一般に、`.R` 拡張子で保存され、実行内容を反映した名前が付けられる。 @@ -94,7 +94,7 @@ knitr::include_graphics("figures/codecheck.png") 詳細は reprex.tidyverse.org を参照。 ``` -\index{reproducibility} +\index{さいげんせい@再現性} このセクションの内容は、あらゆるタイプの R スクリプトに適用できる。 ジオコンピュテーションのためのスクリプトで特に考慮すべき点は、GDAL など外部ライブラリへの依存が多い点である。実際、Chapter \@ref(read-write) のデータ入出力では GDAL をたくさん使用した。 GIS ソフトウェアの依存は、Chaptger \@ref(gis) で解説したようにより多くの特別なジオアルゴリズムを実行するために必要になる。 @@ -138,32 +138,32 @@ if(curl::has_internet()) { ## 幾何学的アルゴリズム {#geometric-algorithms} -アルゴリズム\index{algorithm}は、コンピュータにおける料理のレシピに相当するものと理解することができる。 +アルゴリズム\index{あるごりずむ@アルゴリズム}は、コンピュータにおける料理のレシピに相当するものと理解することができる。 入力に対して実行されると、有用または美味しい出力が得られる指示の完全なセットである。 入力とは、料理においては小麦粉や砂糖などの食材で、アルゴリズムの場合はデータと入力パラメータとなる。 美味しいケーキはレシピの結果であるのと同様に、アルゴリズムの成功は環境や社会的に利点のある出力となりうる。 具体的なケーススタディに入る前に、アルゴリズムとスクリプト (Section \@ref(scripts) ) や関数 (Section \@ref(functions) で説明するように、アルゴリズムを一般化し、他でも使えるようにしたり簡単にする) の関係について簡単に説明する。 -「アルゴリズム」\index{algorithm}という言葉は、9 世紀のバグダッドで、初期の数学教科書である *Hisab al-jabr w'al-muqabala* の出版に端を発している。 +「アルゴリズム」\index{あるごりずむ@アルゴリズム}という言葉は、9 世紀のバグダッドで、初期の数学教科書である *Hisab al-jabr w'al-muqabala* の出版に端を発している。 この本はラテン語に翻訳され、人気を博しただけでなく、著者の名字である [al-Khwārizmī](https://en.wikipedia.org/wiki/Muhammad_ibn_Musa_al-Khwarizmi) が「科学用語として不滅の名声を得て、Alchoarismi、Algorismi、そして最終的には algorithm になった」 [@bellos_alex_2011] 。 -コンピューティング時代には、アルゴリズム\index{algorithm}は、問題を解決する一連のステップを指し、その結果、あらかじめ定義された出力が得られる。 +コンピュータ時代には、アルゴリズム\index{あるごりずむ@アルゴリズム}は、問題を解決する一連のステップを指し、その結果、あらかじめ定義された出力が得られる。 入力は、適切なデータ構造で正式に定義されなければならない [@wise_gis_2001]。 -アルゴリズムは、多くの場合、コードで実装される前に、処理の目的を示すフローチャートや疑似コード\index{pseudocode}として開始される。 +アルゴリズムは、多くの場合、コードで実装される前に、処理の目的を示すフローチャートや疑似コード\index{ぎじこーど@疑似コード}として開始される。 使い勝手をよくするために、一般的なアルゴリズムは関数内にパッケージ化されていることが多く、(関数のソースコードを見ない限り) 一部または全部の手順が隠されている場合がある (Section \@ref(functions) を参照)。 -Chapter \@ref(gis) で見たような ジオアルゴリズム\index{geoalgorithm}は、地理的なデータを取り込み、一般的には地理的な結果を返すアルゴリズムである (同じものを表す別の用語として、GIS アルゴリズム、幾何学的アルゴリズムがある)。 +Chapter \@ref(gis) で見たような ジオアルゴリズム\index{あるごりずむ@アルゴリズム}は、地理的なデータを取り込み、一般的には地理的な結果を返すアルゴリズムである (同じものを表す別の用語として、GIS アルゴリズム、幾何学的アルゴリズムがある)。 簡単そうに聞こえるだろうが、このテーマは奥が深く、計算幾何学という学問分野全体がその研究に専念している [@berg_computational_2008]。このテーマに関する書籍も多数出版されている。 例えば、@orourke_computational_1998 は、再現可能で自由に利用できる C コードを用いて、徐々に難しくなる幾何学的アルゴリズムの範囲を紹介している。 -幾何学的アルゴリズムの例としては、ポリゴンの重心\index{centroid}を求めるものがある。 -重心\index{centroid}の計算には多くのアプローチがあり、中には特定のタイプの[空間データ](https://en.wikipedia.org/wiki/Centroid)に対してのみ機能するものもある。 -本節では、視覚化しやすいアプローチを選択する。ポリゴンを多くの三角形に分割し、それぞれの重心\index{centroid}を求める。このアプローチは、他の重心アルゴリズムと並んで @kaiser_algorithms_1993 によって議論された [簡単な説明は @orourke_computational_1998]。 -コードを書く前に、このアプローチをさらに個別のタスクに分解するのに役立つ (以降、ステップ 1 からステップ 4 と呼ぶが、これらは模式図や疑似コード\index{pseudocode} として提示すこともできる)。 +幾何学的アルゴリズムの例としては、ポリゴンの重心\index{じゅうしん@重心}を求めるものがある。 +重心\index{じゅうしん@重心}の計算には多くのアプローチがあり、中には特定のタイプの[空間データ](https://en.wikipedia.org/wiki/Centroid)に対してのみ機能するものもある。 +本節では、視覚化しやすいアプローチを選択する。ポリゴンを多くの三角形に分割し、それぞれの重心\index{じゅうしん@重心}を求める。このアプローチは、他の重心アルゴリズムと並んで @kaiser_algorithms_1993 によって議論された [簡単な説明は @orourke_computational_1998]。 +コードを書く前に、このアプローチをさらに個別のタスクに分解するのに役立つ (以降、ステップ 1 からステップ 4 と呼ぶが、これらは模式図や疑似コード\index{ぎじこーど@疑似コード} として提示すこともできる)。 1. ポリゴンを連続した三角形に分割する -2. 各三角形の重心\index{centroid}を求める +2. 各三角形の重心\index{じゅうしん@重心}を求める 3. それぞれの三角形の面積を求める -4. 三角形の中心点の面積加重平均\index{centroid}を求める +4. 三角形の中心点の面積加重平均\index{じゅうしん@重心}を求める 一見、簡単そうに見えるが、言葉をコードに変換するには、入力に制約がある場合でも、試行錯誤を繰り返しながら作業を進める必要がある。 このアルゴリズムは、180°以上の内角を持たない凸ポリゴンに対してのみ動作し、星形は使用できない (パッケージの **decido** と **sfdct** は外部ライブラリを使用して非凸ポリゴンを三角測量できる。[geocompx.org](https://geocompx.org/) の [algorithm](https://geocompx.github.io/geocompkg/articles/algorithm.html) vignetteに示されている。)。 @@ -184,7 +184,7 @@ poly_mat = cbind(x_coords, y_coords) ``` これで、例となるデータセットができたので、上記のステップ 1 に着手する準備が整った。 -以下のコードでは、1つの三角形 (`T1`) を作成して、この方法を示している。また、 [数式](https://math.stackexchange.com/a/1702606) $1/3(a + b + c)$ ($a$ から $c$ は三角形の頂点を表す座標) に基づいて重心\index{centroid}を計算するステップ 2 も示している。 +以下のコードでは、1 つの三角形 (`T1`) を作成して、この方法を示している。また、 [数式](https://math.stackexchange.com/a/1702606) $1/3(a + b + c)$ ($a$ から $c$ は三角形の頂点を表す座標) に基づいて重心\index{じゅうしん@重心}を計算するステップ 2 も示している。 ```{r 10-algorithms-5} # 原点を作成: @@ -199,7 +199,6 @@ C1 = (T1[1,] + T1[2,] + T1[3,]) / 3 C1_alternative = (T1[1, , drop = FALSE] + T1[2, , drop = FALSE] + T1[3, , drop = FALSE]) / 3 ``` - ```{r polymat, echo=FALSE, fig.cap="ポリゴン重心計算問題の例", fig.height="100", warning=FALSE} # 最初のプロット。おそらく削除。 old_par = par(pty = "s") @@ -237,8 +236,8 @@ abs(T1[1, 1] * (T1[2, 2] - T1[3, 2]) + ステップ 4 では、ステップ 2 と 3 を 1 つの三角形だけでなく、すべての三角形に対して行う必要がある (上の例)。 このため、ポリゴンを表すすべての三角形を作成するためのイテレーション (繰り返し) が必要である。Figure \@ref(fig:polycent) に示す。 -`lapply()`\index{loop!lapply} と `vapply()`\index{loop!vapply} が各三角形の反復処理に使われているのは、base R で簡潔な解が得られるからである。^[ -`vapply()`\index{loop!vapply} のドキュメントについては `?lapply` を参照。イテレーションについては Chapter \@ref(location) を参照。 +`lapply()`\index{るーぷしょり@ループ処理!lapply} と `vapply()`\index{るーぷしょり@ループ処理!vapply} が各三角形の反復処理に使われているのは、Base R で簡潔な解が得られるからである。^[ +`vapply()`\index{るーぷしょり@ループ処理!vapply} のドキュメントについては `?lapply` を参照。イテレーションについては Chapter \@ref(location) を参照。 ] ```{r 10-algorithms-7} @@ -262,9 +261,9 @@ A = vapply(T_all, function(x) { source("code/11-polycent.R") ``` -これで、ステップ4を完了し、総面積を `sum(A)` で計算し、ポリゴンの重心\index{centroid}座標を `weighted.mean(C [, 1] , A)` と `weighted.mean(C [, 2] , A)` でポリゴンの重心座標を計算する (注意深い読者のための練習: これらのコマンドが動作することを確認してみよう)。 -アルゴリズム\index{algorithm}とスクリプトの関連性を示すために、このセクションの内容を凝縮して `11-centroid-alg.R` とした。 -Section \@ref(scripts) の最後で、このスクリプトが正方形の重心\index{centroid}を計算する方法を見た。 +これで、ステップ4を完了し、総面積を `sum(A)` で計算し、ポリゴンの重心\index{じゅうしん@重心}座標を `weighted.mean(C [, 1] , A)` と `weighted.mean(C [, 2] , A)` でポリゴンの重心座標を計算する (注意深い読者のための練習: これらのコマンドが動作することを確認してみよう)。 +アルゴリズム\index{あるごりずむ@アルゴリズム}とスクリプトの関連性を示すために、このセクションの内容を凝縮して `11-centroid-alg.R` とした。 +Section \@ref(scripts) の最後で、このスクリプトが正方形の重心\index{じゅうしん@重心}を計算する方法を見た。 アルゴリズムをスクリプト化することの素晴らしい点は、新しい `poly_mat` オブジェクト上で動作することである (`st_centroid()` を参照してこれらの結果を検証するには、以下の演習を参照)。 ```{r 10-algorithms-8} @@ -273,26 +272,26 @@ source("code/11-centroid-alg.R") 上記の例では、低レベルの地理的な操作は、base R で第一原理から開発することができることが示されている。 また、すでに試行錯誤したソリューションが存在する場合、車輪の再発明をする価値はないことも示している。 -もし、ポリゴンの重心\index{centroid}を求めるだけなら、`poly_mat` を **sf** オブジェクトとして表現し、代わりに既存の `sf::st_centroid()` 関数を使用する方が早かっただろう。 +もし、ポリゴンの重心\index{じゅうしん@重心}を求めるだけなら、`poly_mat` を **sf** オブジェクトとして表現し、代わりに既存の `sf::st_centroid()` 関数を使用する方が早かっただろう。 しかし、第一原理でアルゴリズムを書くことの大きな利点は、プロセスのすべての段階を理解できることであり、他の人のコードを使うときには保証されないことである。 さらに考慮すべきは性能である。R は C++\index{C++} のような低レベルの言語と比較すると数値計算が遅いことが多く、最適化が困難である (Section \@ref(software-for-geocomputation) 参照)。 新しい手法の開発を目的とするのであれば、計算効率を優先させるべきではない。 これは、「早すぎる最適化はプログラミングにおける諸悪の根源 (あるいは少なくともそのほとんど)」という言葉に集約される [@knuth_computer_1974]。 -アルゴリズム\index{algorithm}開発は大変な作業である。 -このことは、base R\index{R} を使った重心\index{centroid}アルゴリズム\index{algorithm}の開発に費やした作業量から明らかである。このアルゴリズムは、実世界での応用が限られている問題に対する一つの、むしろ非効率的なアプローチに過ぎない。というのも、実際には凸ポリゴンは珍しいからである。 +アルゴリズム\index{あるごりずむ@アルゴリズム}開発は大変な作業である。 +このことは、Base R\index{R} を使った重心\index{じゅうしん@重心}アルゴリズム\index{あるごりずむ@アルゴリズム}の開発に費やした作業量から明らかである。このアルゴリズムは、実世界での応用が限られている問題に対する一つの、むしろ非効率的なアプローチに過ぎない。というのも、実際には凸ポリゴンは珍しいからである。 この経験は、GEOS\index{GEOS} や CGAL\index{CGAL} (計算幾何学アルゴリズムライブラリ, Computational Geometry Algorithms Library) など、高速に動作し、かつ幅広い入力ジオメトリタイプに対応する低レベル地理ライブラリの理解につながるはずである。 -このようなライブラリのオープンソース化の大きな利点は、そのソースコード\index{source code}が、研究、理解、(技術と自信があれば) 改変のために容易に利用できることである。^[ +このようなライブラリのオープンソース化の大きな利点は、そのソースコード\index{そーすこーど@ソースコード}が、研究、理解、(技術と自信があれば) 改変のために容易に利用できることである。^[ CGAL\index{CGAL} の関数 `CGAL::centroid()` は、https://doc.cgal.org/latest/Kernel_23/group__centroid__grp.html で説明しているように、実際には 7 つのサブ関数で構成されており、幅広い入力データ型に対応できるようになっているが、私たちが作成したソリューションは特定の入力データ型にのみ対応する。 GEOS\index{GEOS} の関数 `Centroid::getCentroid()` の基礎となるソースコードは、https://github.com/libgeos/geos/search?q=getCentroid で見ることができる。 ] ## 関数 {#functions} -アルゴリズムと同様に \index{algorithm}、関数は入力を受け取り、出力を返す。 -しかし、関数 \index{function} は、「レシピ」そのものではなく、特定のプログラミング言語での実装を指している。 -R では、関数 \index{function} はそれ自体がオブジェクトであり、モジュール方式で作成したり結合したりすることができる。 -例えば、重心 \index{centroid} 生成アルゴリズム \index{algorithm} のステップ 2 を引き受ける関数を以下のように作成することができる。 +アルゴリズムと同様に \index{あるごりずむ@アルゴリズム}、関数は入力を受け取り、出力を返す。 +しかし、関数\index{かんすう@関数}は、「レシピ」そのものではなく、特定のプログラミング言語での実装を指している。 +R では、関数\index{かんすう@関数}はそれ自体がオブジェクトであり、モジュール方式で作成したり結合したりすることができる。 +例えば、重心\index{じゅうしん＠重心}生成アルゴリズム \index{\index{かんすう@関数}のステップ 2 を引き受ける関数を以下のように作成することができる。 ```{r 10-algorithms-9} t_centroid = function(x) { @@ -302,7 +301,7 @@ t_centroid = function(x) { 上記の例では、[関数](https://adv-r.hadley.nz/functions.html)の2つの重要な構成要素を示している。 1) 関数の *body* 、中括弧内のコードで、関数が入力に対して何をするかを定義する。 2) *formal* 、関数が扱う引数のリスト。この場合は `x` である (3番目の重要なコンポーネント、環境はこのセクションの範囲外)。 -デフォルトでは、関数は最後に計算したオブジェクトを返す (`t_centroid()` の場合は重心\index{centroid}の座標)。^[ +デフォルトでは、関数は最後に計算したオブジェクトを返す (`t_centroid()` の場合は重心\index{じゅうしん@重心}の座標)。^[ また、関数の本文に `return(output)` を追加することで、関数の出力を明示的に設定することができる。 `output` は返すべき結果である。 ] @@ -318,7 +317,7 @@ environment(t_centroid) t_centroid(T1) ``` -また、三角形の面積を計算する関数\index{function}を作成することができる。ここでは、`t_area()` と名付ける。 +また、三角形の面積を計算する関数\index{かんすう@関数}を作成することができる。ここでは、`t_area()` と名付ける。 ```{r 10-algorithms-12} t_area = function(x) { @@ -332,13 +331,13 @@ t_area = function(x) { なお、この関数を作成した後は、1行のコードで三角形の面積を計算できるようになり、冗長なコードの重複を避けることができる。 関数は、コードを一般化するためのメカニズムである。 -新たに作成した関数\index{function} `t_area()` は、これまで使ってきた「三角行列」データ構造と同じ寸法を持つと仮定した任意のオブジェクト `x` を受け取り、その面積を返すもので、`T1` で図示すと次のようになる。 +新たに作成した関数\index{かんすう@関数} `t_area()` は、これまで使ってきた「三角行列」データ構造と同じ寸法を持つと仮定した任意のオブジェクト `x` を受け取り、その面積を返すもので、`T1` で図示すと次のようになる。 ```{r 10-algorithms-13} t_area(T1) ``` -関数\index{function}を使って、高さ1、底辺3 の新しい三角行列の面積を求めることで、その一般化可能性を検証することができる。 +関数\index{かんすう@関数}を使って、高さ1、底辺3 の新しい三角行列の面積を求めることで、その一般化可能性を検証することができる。 ```{r 10-algorithms-14} t_new = cbind(x = c(0, 3, 3, 0), @@ -348,9 +347,9 @@ t_area(t_new) 関数の便利な点は、モジュール化されていることである。 出力が何であるかが分かっていれば、ある関数を別の関数の構成要素として利用することができる。 -したがって、関数 `t_centroid()` と `t_area()` は、スクリプト `11-centroid-alg.R`というより大きな関数\index{function}のサブコンポーネントとして使うことができる。このスクリプトは、任意の凸ポリゴンの面積を計算する。 +したがって、関数 `t_centroid()` と `t_area()` は、スクリプト `11-centroid-alg.R`というより大きな関数\index{かんすう@関数}のサブコンポーネントとして使うことができる。このスクリプトは、任意の凸ポリゴンの面積を計算する。 以下のコードでは、凸ポリゴンに対する `sf::st_centroid()` の動作を模倣する関数 `poly_centroid()` を作成している。^[ -なお、作成した関数は、`lapply()`\index{loop!lapply} と `vapply()`\index{loop!vapply} の関数呼び出しで繰り返し呼び出されている。 +なお、作成した関数は、`lapply()`\index{るーぷしょり@ループ処理!lapply} と `vapply()`\index{るーぷしょり@ループ処理!vapply} の関数呼び出しで繰り返し呼び出されている。 ] ```{r 10-algorithms-15} @@ -389,7 +388,7 @@ poly_centroid = function(poly_mat, output = "matrix") { poly_centroid(poly_mat) ``` -`poly_centroid()` などの関数\index{function}はさらに拡張して、さまざまなタイプの出力を提供することができる。 +`poly_centroid()` などの関数\index{かんすう@関数}はさらに拡張して、さまざまなタイプの出力を提供することができる。 例えば、結果をクラス `sfg` のオブジェクトとして返すには、結果を返す前に、「ラッパー」関数を用いて `poly_centroid()` の出力を変更することができる。 ```{r 10-algorithms-18} @@ -408,12 +407,12 @@ identical(poly_centroid_sfg(poly_mat), sf::st_centroid(poly_sfc)) ## プログラミング {#programming} -この章では、スクリプトからアルゴリズム\index{algorithm}という厄介なトピックを経由して関数へと移った。 +この章では、スクリプトからアルゴリズム\index{あるごりずむ@アルゴリズム}という厄介なトピックを経由して関数へと移った。 抽象的な議論だけでなく、具体的な問題を解決するために、それぞれの実用例を作成した。 -- スクリプト `11-centroid-alg.R` が紹介され、「ポリゴンマトリックス」での実演が行われることとした。 -- このスクリプトを動作させるための個々のステップは、アルゴリズム \index{algorithm}、計算レシピとして記述された。 -- アルゴリズムを一般化するために、 \index{algorithm} をモジュール関数に変換し、最終的にそれらを組み合わせて、前節の関数 `poly_centroid()` を作成した。 +- スクリプト `11-centroid-alg.R` を導入し、「ポリゴンマトリックス」で実行した。 +- このスクリプトを動作させるための個々のステップは、アルゴリズム\index{あるごりずむ@アルゴリズム}、つまり計算レシピとして記述した。 +- アルゴリズムを一般化するために、このアルゴリズムをモジュール関数に変換し、最終的にそれらを組み合わせて、前節の関数 `poly_centroid()` を作成した。\index{あるごりずむ@アルゴリズム} これらのステップは、簡単なことである。 しかし、プログラミングの技術とは、スクリプト、アルゴリズム、関数を組み合わせて、性能の良いシステムにすることなのである。 @@ -421,8 +420,8 @@ identical(poly_centroid_sfg(poly_mat), sf::st_centroid(poly_sfc)) この本を読んでいるほとんどの人がそうであるように、プログラミングの初心者であれば、前のセクションの結果を追って再現できることは、大きな達成感を得ることができるはずである。 プログラミングができるようになるまでには、何時間もかけて熱心に勉強し、練習する必要がある。 -新しいアルゴリズム\index{algorithm}を効率的に実装する際の課題は、実運用で使用することを意図していない単純な関数を作成するために費やした作業量を考慮すると、はるかに遠い。現在の状態では、`poly_centroid()` はほとんどの (非凸) ポリゴンで失敗する! -ここから生じる疑問は、関数\index{function}をどのように一般化するかということである。 +新しいアルゴリズム\index{あるごりずむ@アルゴリズム}を効率的に実装する際の課題は、実運用で使用することを意図していない単純な関数を作成するために費やした作業量を考慮すると、はるかに遠い。現在の状態では、`poly_centroid()` はほとんどの (非凸) ポリゴンで失敗する! +ここから生じる疑問は、関数\index{かんすう@関数}をどのように一般化するかということである。 (1) 非凸ポリゴンを三角測量する方法を探す (この章をサポートするオンライン記事 [algorithm](https://geocompx.github.io/geocompkg/articles/algorithm.html) で扱っている話題) と (2) 三角メッシュに依存しない他の重心のアルゴリズムを調べるという2つの選択肢がある。 もっと大きな疑問は、高性能なアルゴリズムがすでに実装され、`st_centroid()` のような関数にパッケージされているのに、ソリューションをプログラミングする価値があるのだろうか、ということである。 @@ -453,6 +452,6 @@ identical(poly_centroid_sfg(poly_mat), sf::st_centroid(poly_sfc)) ## 演習 {#ex-algorithms} ```{r, echo=FALSE, results='asis'} -res = knitr::knit_child('_11-ex.Rmd', quiet = TRUE, options = list(include = FALSE, eval = FALSE)) +res = knitr::knit_child('_11-ex-ja.Rmd', quiet = TRUE, options = list(include = FALSE, eval = FALSE)) cat(res, sep = '\n') ``` diff --git a/12-spatial-cv-ja.Rmd b/12-spatial-cv-ja.Rmd index f7bced9..4e6085d 100644 --- a/12-spatial-cv-ja.Rmd +++ b/12-spatial-cv-ja.Rmd @@ -10,8 +10,8 @@ knitr::opts_chunk$set(warning = FALSE) ## 必須パッケージ {- #prerequisites-12} -本章では、Chapter \@ref(spatial-class) から Chapter \@ref(reproj-geo-data) までの内容を学習し、演習を行うなどして、地理データ解析 \index{geographic data analysis} に習熟していることを前提としている。 -一般化線形モデル (GLM) \index{GLM}と機械学習\index{machine learning} に精通していることを強く推奨する [例えば @zuur_mixed_2009;@james_introduction_2013]。 +本章では、Chapter \@ref(spatial-class) から Chapter \@ref(reproj-geo-data) までの内容を学習し、演習を行うなどして、地理データ解析\index{ちりでーたかいせき@地理データ解析} に習熟していることを前提としている。 +一般化線形モデル (GLM)\index{GLM}と機械学習\index{きかいがくしゅう@機械学習} に精通していることを強く推奨する [例えば @zuur_mixed_2009;@james_introduction_2013]。 この章では、以下のパッケージを使用する。^[ パッケージ **GGally**、**lgr**、**kernlab**、**mlr3measures**、**paradox**、**pROC**、**progressr**、**spDataLarge** もインストールする必要があるが、ロードしておく必要はない。 @@ -21,49 +21,53 @@ knitr::opts_chunk$set(warning = FALSE) library(sf) library(terra) library(dplyr) -library(future) -library(lgr) -library(mlr3) -library(mlr3learners) -library(mlr3extralearners) -library(mlr3spatiotempcv) -library(mlr3tuning) -library(mlr3viz) -library(progressr) +library(future) # 並列処理 +library(lgr) # logging framework for R +library(mlr3) # 機械学習アルゴリズムへの統一インターフェイス +library(mlr3learners) # 最重要の機械学習アルゴリズム +library(mlr3extralearners) # その他の機械学習アルゴリズム +library(mlr3spatiotempcv) # 時空間リサンプリング +library(mlr3tuning) # ハイパーパラメータのチューニング +library(mlr3viz) # mlr3 オブジェクトのプロット関数 +library(progressr) # 進捗状況を報告 ``` データは必要に応じて読み込む。 ## イントロダクション {#intro-cv1} -統計的学習\index{statistical learning}は、データのパターンを特定し、そのパターンから予測するための統計的・計算的モデルの使用に関するものである。 -その起源から、統計的学習\index{statistical learning}は R\index{R} の の大きな強みの一つである ( Section \@ref(software-for-geocomputation) 参照)。^[ +統計的学習\index{とうけいてきがくしゅう@統計的学習}は、データのパターンを特定し、そのパターンから予測するための統計的・計算的モデルの使用に関するものである。 +その起源から、統計的学習\index{とうけいてきがくしゅう@統計的学習}は R\index{R} の の大きな強みの一つである ( Section \@ref(software-for-geocomputation) 参照)。^[ 地理データに統計的手法を適用することは、地理統計学、空間統計学、ポイントパターン解析の分野において、何十年にもわたって活発な研究テーマとなっている [@diggle_modelbased_2007; @gelfand_handbook_2010; @baddeley_spatial_2015]。 ] -統計的学習\index{statistical learning}とは、統計学\index{statistics}と機械学習\index{machine learning}の手法を組み合わせたもので、教師あり手法と教師なし手法に分類される。 +統計的学習\index{とうけいてきがくしゅう@統計的学習}とは、統計学\index{とうけいがく@統計学}と機械学習\index{きかいがくしゅう@機械学習}の手法を組み合わせたもので、教師あり手法と教師なし手法に分類される。 どちらも物理学、生物学、生態学から地理学、経済学に至るまで、ますます多くの分野で利用されるようになっている [@james_introduction_2013]。 この章では、クラスタリング\index{clustering}のような教師なし技術とは対照的に、学習データセットが存在する教師あり技術に焦点を当てていきたい。 応答変数は、二値 (地滑りの発生など)、カテゴリ (土地利用)、整数 (種の豊富さ数)、数値 (土壌酸性度の測定された pH 値) のいずれでもよい。 教師あり技術は、オブザベーションのサンプルについて既知の応答と、1つまたは複数の予測変数の間の関係をモデル化する。 -多くの機械学習\index{machine learning}において、研究の主な目的は優れた予測を行うことである。 -機械学習が「ビッグデータ」\index{big data}の時代に繁栄しているのは、その手法が入力変数に関する仮定をほとんど必要とせず、巨大なデータセットを扱えるからである。 +多くの機械学習\index{きかいがくしゅう@機械学習}において、研究の主な目的は優れた予測を行うことである。 +機械学習が「ビッグデータ」\index{びっぐでーた@ビッグデータ}の時代に繁栄しているのは、その手法が入力変数に関する仮定をほとんど必要とせず、巨大なデータセットを扱えるからである。 機械学習は、将来の顧客行動予測、推薦サービス (音楽、映画、次に買うもの)、顔認識、自律走行、テキスト分類、予知保全 (インフラ、産業) などのタスクに資するものである。 -この章では、地すべりの (空間) 予測という事例をもとに説明する。 -この応用例は、Chapter \@ref(intro) で定義されているジオコンピューティングの応用的な性質とリンクしており、機械学習\index{machine learning}が、予測を唯一の目的とする場合に統計学の分野\index{statistics}から借用する方法を示している。 -そこで、この章では、まず、一般化線形モデル\index{GLM} [@zuur_mixed_2009] の助けを借りて、モデリングとクロスバリデーション\index{cross-validation}の概念を紹介する。 -これを踏まえて、この章では、より典型的な機械学習\index{machine learning} アルゴリズム \index{algorithm}、すなわちサポートベクタマシン (Support Vector Machine, SVM) \index{SVM}を実装している。 -モデルの**予測性能**は、地理データが特殊であることを考慮した空間クロスバリデーション (Cross Validation, CV) \index{cross-validation!spatial CV}を用いて評価していこう。 +この章では、地すべりの発生モデルという事例をもとに説明する。 +この応用例は、Chapter \@ref(intro) で定義されているジオコンピューティングの応用的な性質とリンクしており、機械学習\index{きかいがくしゅう@機械学習}が、予測を唯一の目的とする場合に統計学の分野\index{とうけいがく@統計学}から借用する方法を示している。 +そこで、この章では、まず、一般化線形モデル\index{GLM} [@zuur_mixed_2009] の助けを借りて、モデリングとクロスバリデーション\index{くろすばりでーしょん@クロスバリデーション}の概念を紹介する。 +これを踏まえて、この章では、より典型的な機械学習\index{きかいがくしゅう@機械学習} アルゴリズム \index{あるごりずむ@アルゴリズム}、すなわちサポートベクタマシン (Support Vector Machine, SVM) \index{SVM}を実装している。 +モデルの**予測性能**は、地理データが特殊であることを考慮した空間クロスバリデーション (Cross Validation, CV) \index{くろすばりでーしょん@クロスバリデーション!くうかん CV@空間 CV}を用いて評価していこう。 -CV\index{cross-validation} データセットをトレーニングセットとテストセットに (繰り返し) 分割することで、モデルが新しいデータに対して汎化する能力を決定する。 +CV\index{くろすばりでーしょん@クロスバリデーション} データセットをトレーニングセットとテストセットに (繰り返し) 分割することで、モデルが新しいデータに対して汎化する能力を決定する。 学習データを使ってモデルを適合させ、テストデータに対して予測したときの性能をチェックする。 -CV は過適合\index{overfitting}を検出するのに役立つ。なぜなら、学習データをあまりに忠実に予測するモデル (ノイズ) は、テストデータでのパフォーマンスが低くなる傾向があるからである。 +CV は過適合\index{かてきごう@過適合}を検出するのに役立つ。なぜなら、学習データをあまりに忠実に予測するモデル (ノイズ) は、テストデータでのパフォーマンスが低くなる傾向があるからである。 空間データをランダムに分割することで、テスト点と空間的に隣接する学習点を得ることができる。 -空間的に自己相関\index{autocorrelation!spatial}していると、このシナリオではテストとトレーニングのデータセットが独立しておらず、結果として CV \index{cross-validation} は過適合\index{overfitting}の可能性を検出できなくなる。 -空間 CV\index{cross-validation!spatial CV} は、本章の**中心**テーマであり、この問題を軽減する。。 +空間的に自己相関\index{autocorrelation!spatial}していると、このシナリオではテストとトレーニングのデータセットが独立しておらず、結果として CV \index{くろすばりでーしょん@クロスバリデーション} は過適合\index{かてきごう@過適合}の可能性を検出できなくなる。 +空間 CV\index{くろすばりでーしょん@クロスバリデーション!くうかん CV@空間 CV} は、本章の**中心**テーマであり、この問題を軽減する。。 + +繰り返しになるが、この章ではモデルの**予測性能**に焦点を当てる。 +予測地図は**扱わない**。 +これは、Chapter \@ref(eco) で扱う。 ## ケーススタディ：地すべりの発生しやすさ {#case-landslide} @@ -75,12 +79,12 @@ data("lsl", "study_mask", package = "spDataLarge") ta = terra::rast(system.file("raster/ta.tif", package = "spDataLarge")) ``` -上記のコードでは、`lsl`、`sf` という名前の `data.frame`、`study_mask` という名前の という `sf` オブジェクト、そして `ta` (Section \@ref(raster-classes) を参照) という名前の地形属性ラスタ `SpatRaster` という3つのオブジェクトをロードしている。 +上記のコードでは、`lsl`、`sf` という名前の `data.frame`、`study_mask` という名前の `sf` オブジェクト、そして `ta` (Section \@ref(raster-classes) を参照) という名前の地形属性ラスタ `SpatRaster` という 3 つのオブジェクトをロードしている。 `lsl` は要因列 `lslpts` を含み、`TRUE`は観測された地すべり「開始点」に対応し、座標は列 `x` と `y` に格納されている。 ^[ 地すべり開始点は、地すべりポリゴンの崖に位置する。詳細は @muenchow_geomorphic_2012 を参照。 ] -`summary(lsl$lslpts)` に示すように、地すべり地点が175箇所、非地すべり地点が175箇所ある。 -非地すべり点175点は、地すべりポリゴン周辺の小さな緩衝地帯の外に位置しなければならないという制約のもと、調査地域からランダムにサンプル化されたものである。 +`summary(lsl$lslpts)` に示すように、地すべり地点が 175 箇所、非地すべり地点が 175 箇所ある。 +非地すべり点 175 点は、地すべりポリゴン周辺の小さな緩衝地帯の外に位置しなければならないという制約のもと、調査地域からランダムにサンプル化されたものである。 ```{r lsl-map, echo=FALSE, out.width="70%", fig.cap="エクアドル南部における地すべり発生地点 (赤) と地すべりの影響を受けていない地点 (青)。", fig.scap="Landslide initiation points."} # library(tmap) @@ -132,7 +136,7 @@ R-GIS ブリッジ (Chapter \@ref(gis) 参照) を用いて地形属性を計算 ## R による従来のモデリング手法 {#conventional-model} 何十もの学習アルゴリズムへの統一的なインタフェースを提供するアンブレラパッケージである **mlr3**\index{mlr3 (package)} パッケージを紹介する (Section \@ref(spatial-cv-with-mlr3)) が、その前に R\index{R} の従来のモデリングインタフェースについて見ておく価値がある。 -この教師あり統計学習\index{statistical learning}の入門は、空間 CV\index{cross-validation!spatial CV} を行うための基礎となり、この後に紹介する **mlr3**\index{mlr3 (package)} のアプローチの把握に貢献する。 +この教師あり統計学習\index{とうけいてきがくしゅう@統計的学習}の入門は、空間 CV\index{くろすばりでーしょん@クロスバリデーション!くうかん CV@空間 CV} を行うための基礎となり、この後に紹介する **mlr3**\index{mlr3 (package)} のアプローチの把握に貢献する。 教師あり学習では、予測変数の関数として応答変数を予測する (Section \@ref(intro-cv))。 R\index{R} では、モデリング関数は通常、数式を使って指定する (R の数式の詳細については、`?formula` を参照)。 @@ -231,11 +235,11 @@ pROC::auc(pROC::roc(lsl$lslpts, fitted(fit))) AUROC\index{AUROC} の値 0.82 は良好な適合性を示している。 しかし、これは完全なデータセットに対して計算したものであるため、楽観的すぎる推定値である。 -偏りを抑えた評価を導き出すためには、クロスバリデーション \index{cross-validation} を用いる必要があり、空間データの場合は空間 CV\index{cross-validation!spatial CV} を利用する必要がある。 +偏りを抑えた評価を導き出すためには、クロスバリデーション \index{くろすばりでーしょん@クロスバリデーション} を用いる必要があり、空間データの場合は空間 CV\index{くろすばりでーしょん@クロスバリデーション!くうかん CV@空間 CV} を利用する必要がある。 ## (空間) クロスバリデーションの紹介 {#intro-cv} -クロスバリデーション\index{cross-validation} は、リサンプリング法\index{resampling}のファミリーに属する [@james_introduction_2013]。 +クロスバリデーション\index{くろすばりでーしょん@クロスバリデーション} は、リサンプリング法\index{resampling}のファミリーに属する [@james_introduction_2013]。 基本的な考え方は、データセットをトレーニングセットとテストセットに分割 (繰り返し) し、トレーニングデータを使ってモデルを適合させ、それをテストセットに適用することである。 予測値とテストセットの既知の応答値を比較することにより (二項式の場合は AUROC\index{AUROC} のような性能指標を使用)、学習した関係を独立したデータに一般化するモデルの能力について、バイアスを低減した評価を得ることができる。 例えば、5倍クロスバリデーションを100回繰り返すとは、データをランダムに5分割 (フォールド) し、各フォールドをテストセットとして1回使用することを意味する (Figure \@ref(fig:partitioning) の上段を参照)。 @@ -246,12 +250,12 @@ AUROC\index{AUROC} の値 0.82 は良好な適合性を示している。 しかし、地理的なデータは特別である。 Chapter \@ref(transport) で見るように、地理学の「第一法則」は、互いに近い地点は、一般に、遠い地点よりも似ているとするものである [@miller_tobler_2004]。 -つまり、従来の CV\index{cross-validation} では学習点とテスト点が近すぎることが多いため、点が統計的に独立していないことになる (Figure \@ref(fig:partitioning) の最初の行を参照)。 +つまり、従来の CV\index{くろすばりでーしょん@クロスバリデーション} では学習点とテスト点が近すぎることが多いため、点が統計的に独立していないことになる (Figure \@ref(fig:partitioning) の最初の行を参照)。 「テスト」観測の近くにある「トレーニング」観測は、一種の「カンニング」を提供することができる。 すなわち、学習データセットでは利用できないはずの情報である。 この問題を軽減するために、オブザベーションを空間的に不連続なサブセットに分割する「空間分割」が使用される (*k*-means クラスタリング\index{clustering!kmeans}でオブザベーションの座標を使用; @brenning_spatial_2012; Figure \@ref(fig:partitioning) の2行目)。 この分割戦略が、従来のCVとの**唯一**の違いである。 -その結果、空間的 CV はモデルの予測性能のバイアスを低減させ、過適合\index{overfitting}を回避するのに役立つ。 +その結果、空間的 CV はモデルの予測性能のバイアスを低減させ、過適合\index{かてきごう@過適合}を回避するのに役立つ。 ```{r partitioning, fig.cap="1 回の繰り返しのクロスバリデーションで選択されたテストおよびトレーニングのオブザベーションの空間的な可視化。ランダム (上段) および空間分割 (下段)。", echo=FALSE, fig.scap="Spatial visualization of selected test and training observations."} knitr::include_graphics("figures/12_partitioning.png") @@ -260,11 +264,11 @@ knitr::include_graphics("figures/12_partitioning.png") ## **mlr3** を用いた空間 CV {#spatial-cv-with-mlr3} \index{mlr3 (package)} -統計的学習\index{statistical learning}のためのパッケージは何十種類もある。例えば [CRAN machine learning task view](https://CRAN.R-project.org/view=MachineLearning) で説明されている。 -クロスバリデーションやハイパーパラメータ ( \index{hyperparameter} ) のチューニング方法など、各パッケージに精通することは時間のかかる作業である。 +統計的学習\index{とうけいてきがくしゅう@統計的学習}のためのパッケージは何十種類もある。例えば [CRAN machine learning task view](https://CRAN.R-project.org/view=MachineLearning) で説明されている。 +クロスバリデーションやハイパーパラメータ\index{はいぱーぱらめーた@ハイパーパラメータ}のチューニング方法など、各パッケージに精通することは時間のかかる作業である。 異なるパッケージのモデル結果を比較するのは、さらに手間がかかる。 これらの問題を解決するために開発されたのが、**mlr3** パッケージとエコシステムである。 -これは「メタパッケージ」として機能し、分類、回帰 \index{regression}、生存時間分析、クラスタリング\index{clustering}など、一般的な教師あり・教師なしの統計学習技術への統一的なインタフェースを提供する [@lang_mlr3_2019; @becker_mlr3_2022] 。 +これは「メタパッケージ」として機能し、分類、回帰 \index{かいき@回帰}、生存時間分析、クラスタリング\index{clustering}など、一般的な教師あり・教師なしの統計学習技術への統一的なインタフェースを提供する [@lang_mlr3_2019; @becker_mlr3_2022] 。 標準化された **mlr3** インターフェースは、8 つの「ビルディングブロック」に基づいている。 Figure \@ref(fig:building-blocks) に示すように、これらは明確な順序を持っている。 @@ -275,15 +279,15 @@ knitr::include_graphics("figures/12_ml_abstraction_crop.png") ``` **mlr3** のモデリングプロセスは、主に3つのステージで構成されている。 -まず、**task** で、データ (応答変数と予測変数を含む) とモデルの種類 (回帰 \index{regression} や分類 \index{classification} など) を指定する。 +まず、**task** で、データ (応答変数と予測変数を含む) とモデルの種類 (回帰 \index{かいき@回帰} や分類 \index{ぶんるい@分類} など) を指定する。 次に、**learner**は、作成されたタスクに適用される特定の学習アルゴリズムを定義する。 第三に、**リサンプリング**アプローチでは、モデルの予測性能、すなわち新しいデータへの汎化能力を評価する (Section \@ref(intro-cv) も参照)。 ### 一般化線形モデル {#glm} GLM\index{GLM} を **mlr3**\index{mlr3 (package)} で使うためには、地すべりデータを含む **task** を作成する必要がある。 -応答は二値 (2カテゴリの変数) で、空間次元を持つので、**mlr3spatiotempcv** パッケージの `TaskClassifST$new()` を使用し、分類 \index{classification}タスクを作成する [@schratz_mlr3spatiotempcv_2021 、非空間 task には `mlr3::TaskClassif$new()`、回帰\index{regression} には `mlr3::TaskRegr$new()` を使用。他の task の詳細は、`?Task` を参照。] 。^[**mlr3** エコシステムは **data.table** と **R6** というクラスを多用する。これらのクラスの仕様を知らなくても **mlr3** を使うことはできるが、知っておくと有利である。 **data.table**, の詳細は https://rdatatable.gitlab.io/data.table/。**R6** について学びたい方は、Advanced R book (日本語版は共立出版「R 言語徹底解説」だが、旧版を元にしている) [@wickham_advanced_2019] の [Chapter 14](https://adv-r.hadley.nz/fp.html) を参照。] -`Task*$new()` 関数の最初の必須引数は、`backend` である。 +応答は二値 (2 カテゴリの変数) で、空間次元を持つので、**mlr3spatiotempcv** パッケージの `as_task_classif_st()` を使用し、分類\index{ぶんるい@分類}タスクを作成する [@schratz_mlr3spatiotempcv_2021 、非空間 task には `mlr3::as_task_classif()`、回帰\index{かいき@回帰} には `as_task_regr()` を使用。他の task の詳細は、`?Task` を参照。] 。^[**mlr3** エコシステムは **data.table** と **R6** というクラスを多用する。これらのクラスの仕様を知らなくても **mlr3** を使うことはできるが、知っておくと有利である。 **data.table**, の詳細は https://rdatatable.gitlab.io/data.table/。**R6** について学びたい方は、Advanced R book (日本語版は共立出版「R 言語徹底解説」だが、旧版を元にしている) [@wickham_advanced_2019] の [Chapter 14](https://adv-r.hadley.nz/fp.html) を参照。] +`as_task_` 関数の最初の必須引数は、`backend` である。 `backend` は、入力データが応答変数と予測変数を含んでいることを想定している。 `target` の引数は応答変数の名前を示し (ここでは `lslpts` )、`positive` は応答変数の2つの因子レベルのうちどちらが地滑り開始点を示すかを決定する (ここでは `TRUE`)。 `lsl` データセットの他のすべての変数が予測因子として機能する。 @@ -293,7 +297,7 @@ GLM\index{GLM} を **mlr3**\index{mlr3 (package)} で使うためには、地す ```{r 12-spatial-cv-11, eval=FALSE} # 1. task を作成 -task = mlr3spatiotempcv::TaskClassifST$new( +task = mlr3spatiotempcv::as_task_classif_st( id = "ecuador_lsl", backend = mlr3::as_data_backend(lsl), target = "lslpts", @@ -306,7 +310,7 @@ task = mlr3spatiotempcv::TaskClassifST$new( なお、`mlr3spatiotempcv::as_task_classif_st()` は、`backend` パラメータの入力として `sf`-オブジェクトも受け付ける。 この場合、引数 `coords_as_features` のみを追加して指定するとよいだろう。 -`lsl` を `sf` -オブジェクトに変換しなかったのは、`TaskClassifST$new()` がバックグラウンドで非空間的な `data.table` オブジェクトに戻してしまうだけだからである。 +`lsl` を `sf` -オブジェクトに変換しなかったのは、`as_task_classif_st()` がバックグラウンドで非空間的な `data.table` オブジェクトに戻してしまうだけだからである。 短時間のデータ探索では、**mlr3viz** パッケージの `autoplot()` 関数は、すべての予測因子に対する応答とすべての予測因子に対する応答をプロットするので便利だろう (図示していない)。 @@ -317,8 +321,8 @@ mlr3viz::autoplot(task, type = "duo") mlr3viz::autoplot(task, type = "pairs") ``` -task を作成したら、使用する統計的学習\index{statistical learning}方式を決定する **学習器** (learner) を選択する必要がある。 -分類\index{classification}**学習器** は `classif.` で、回帰\index{regression} 学習器は `regr.` で始まる (詳しくは `?Learner` を参照)。 +task を作成したら、使用する統計的学習\index{とうけいてきがくしゅう@統計的学習}方式を決定する **学習器** (learner) を選択する必要がある。 +分類\index{ぶんるい@分類}の**学習器** は `classif.` で始まり、回帰\index{かいき@回帰}の学習器は `regr.` で始まる (詳しくは `?Learner` を参照)。 `mlr3extralearners::list_mlr3learners()` は、利用可能なすべての学習器と、どのパッケージから **mlr3** がそれらをインポートしているかをリストアップする (Table \@ref(tab:lrns))。 二値応答変数をモデル化できる学習器について調べるには、次のように実行する。 @@ -348,9 +352,9 @@ knitr::kable(lrns_df, caption.short = "Sample of available learners.", booktabs = TRUE) ``` -これにより、すべての学習器が2クラス問題 (地滑りの有無) をモデル化することができるようになった。 -Section \@ref(conventional-model) で使用され、**mlr3learners** では `classif.log_reg` として実装されている二項分類 \index{classification} 方式を選択することにする。 -さらに、予測の種類を決める `predict.type` を指定する必要がある。 `prob` は、地滑り発生の予測確率を 0 から 1 の間で決定する (これは `type = response` の `predict.glm()` に対応する)。 +これにより、すべての学習器が 2 クラス問題 (地滑りの有無) をモデル化することができるようになった。 +Section \@ref(conventional-model) で使用され、**mlr3learners** では `classif.log_reg` として実装されている二項分類 \index{ぶんるい@分類} 方式を選択することにする。 +さらに、予測の種類を決める `predict.type` を指定する必要がある。`prob` は、地滑り発生の予測確率を 0 から 1 の間で決定する (これは `type = response` の `predict.glm()` に対応する)。 ```{r 12-spatial-cv-13, eval=FALSE} # 2. 学習器を指定 @@ -386,9 +390,9 @@ identical(fit$coefficients, learner$model$coefficients) **mlr3**\index{mlr3 (package)} でモデリングするためのセットアップ手順は、面倒に思えるだろう。 しかし、この一つのインターフェースで、`mlr3extralearners::list_mlr3learners()` が示す 130 種類以上の学習器にアクセスできることを思い出してほしい。各学習器のインターフェースを学ぶことはもっと退屈である。 -さらに、リサンプリング技術の簡単な並列化と、機械学習のハイパーパラメータ\index{hyperparameter}を調整できることも利点である ( Section \@ref(svm) を参照)。 +さらに、リサンプリング技術の簡単な並列化と、機械学習のハイパーパラメータ\index{はいぱーぱらめーた@ハイパーパラメータ}を調整できることも利点である ( Section \@ref(svm) を参照)。 最も重要なことは、**mlr3spatiotempcv** [@schratz_mlr3spatiotempcv_2021] の (空間) リサンプリングは簡単で、リサンプリング法の指定と実行という2つのステップを追加するだけでよいということである。 -100 回繰り返される5回空間 CV\index{cross-validation!spatial CV} : `task` で提供された座標に基づいて5つのパーティションが選ばれ、パーティショニングは100回繰り返される。 [^13] +100 回繰り返される5回空間 CV\index{くろすばりでーしょん@クロスバリデーション!くうかん CV@空間 CV} : `task` で提供された座標に基づいて5つのパーティションが選ばれ、パーティショニングは100回繰り返される。 [^13] [^13]: @@ -438,7 +442,7 @@ mean(score_spcv_glm$classif.auc) |> round(2) ``` -これらの結果を整理するために、100 回繰り返した 5 回非空間クロスバリデーションの AUROC\index{AUROC} 値と比較してみよう (Figure \@ref(fig:boxplot-cv) ; 非空間CV\index{cross-validation}のコードはここでは示さないが、演習セクションで検討する)。 +これらの結果を整理するために、100 回繰り返した 5 回非空間クロスバリデーションの AUROC\index{AUROC} 値と比較してみよう (Figure \@ref(fig:boxplot-cv) ; 非空間CV\index{くろすばりでーしょん@クロスバリデーション}のコードはここでは示さないが、演習セクションで検討する)。 予想通り (Section \@ref(intro-cv) 参照)、空間交差検証の結果は、従来の交差検証アプローチよりも平均して低い AUROC 値をもたらし、空間自己相関\index{autocorrelation!spatial}のため、後者の空間自己相関\index{autocorrelation!spatial}による楽観的な予測性能が強調された。 ```{r boxplot-cv, echo=FALSE, message=FALSE, out.width="75%", fig.cap="空間CV と従来の 100 回繰り返し 5 回 CV におけるGLM AUROC 値の差を示す箱ひげ図。", fig.scap="Boxplot showing AUROC values."} @@ -458,7 +462,7 @@ ggplot2::ggplot(data = score[learner_id == "classif.log_reg"], ### 機械学習のハイパーパラメータの空間的チューニング {#svm} -Section \@ref(intro-cv) では、統計的学習\index{statistical learning}の一環として、機械学習\index{machine learning}を導入した 。 +Section \@ref(intro-cv) では、統計的学習\index{とうけいてきがくしゅう@統計的学習}の一環として、機械学習\index{きかいがくしゅう@機械学習}を導入した 。 もう一度確認しよう。[Jason Brownlee](https://machinelearningmastery.com/linear-regression-for-machine-learning/) による機械学習の以下の定義に従う。 > 機械学習、より具体的には予測モデリングの分野では、説明可能性を犠牲にして、モデルの誤差を最小化すること、あるいは可能な限り正確な予測を行うことに主眼が置かれている。 @@ -466,21 +470,21 @@ Section \@ref(intro-cv) では、統計的学習\index{statistical learning}の Section \@ref(glm) では、GLM を用いて地滑りしやすさを予測した。 ここでは、同じ目的のためにサポートベクタマシン (SVM) \index{SVM}を紹介する。 -ランダムフォレスト\index{random forest}モデルは SVM よりも人気があるだろう。しかし、ハイパーパラメータ\index{hyperparameter}のチューニングがモデル性能に与えるプラスの効果は、SVM の場合の方が顕著である [@probst_hyperparameters_2018]。 +ランダムフォレスト\index{random forest}モデルは SVM よりも人気があるだろう。しかし、ハイパーパラメータ\index{はいぱーぱらめーた@ハイパーパラメータ}のチューニングがモデル性能に与えるプラスの効果は、SVM の場合の方が顕著である [@probst_hyperparameters_2018]。 本節では、 (空間) ハイパーパラメータのチューニングが主な目的であるため、SVM を用いることにする。 ランダムフォレストモデルを適用したい方は、この章を読んでから Chapter \@ref(eco) に進むことを勧める。この章では、現在取り上げられている概念と技術を応用して、ランダムフォレストモデルに基づく空間分布図を作成する方法を説明する。 -SVM\index{SVM} クラスを分離するための最適な「超平面」を探索し (分類 \index{classification} の場合)、特定のハイパーパラメータで「カーネル」を推定して、クラス間の非線形境界を作成する [@james_introduction_2013]。 -機械学習には、ハイパーパラメータ\index{hyperparameter}とパラメータがある。 -パラメータはデータから推定できるが、ハイパーパラメータ\index{hyperparameter}は学習開始前に設定しなければならない (mlr3 本の[machine mastery blog](https://machinelearningmastery.com/difference-between-a-parameter-and-a-hyperparameter/)と[hyperparameter optimization chapter](https://mlr3book.mlr-org.com/chapters/chapter4/hyperparameter_optimization.html) も参照)。 -最適なハイパーパラメータ\index{hyperparameter}は、通常、クロスバリデーション法を用いて定義された範囲内で決定する。 +SVM\index{SVM} クラスを分離するための最適な「超平面」を探索し (分類 \index{ぶんるい@分類} の場合)、特定のハイパーパラメータで「カーネル」を推定して、クラス間の非線形境界を作成する [@james_introduction_2013]。 +機械学習には、ハイパーパラメータ\index{はいぱーぱらめーた@ハイパーパラメータ}とパラメータがある。 +パラメータはデータから推定できるが、ハイパーパラメータ\index{はいぱーぱらめーた@ハイパーパラメータ}は学習開始前に設定しなければならない (mlr3 本の[machine mastery blog](https://machinelearningmastery.com/difference-between-a-parameter-and-a-hyperparameter/)と[hyperparameter optimization chapter](https://mlr3book.mlr-org.com/chapters/chapter4/hyperparameter_optimization.html) も参照)。 +最適なハイパーパラメータ\index{はいぱーぱらめーた@ハイパーパラメータ}は、通常、クロスバリデーション法を用いて定義された範囲内で決定する。 これをハイパーパラメータチューニングという。 **kernlab** が提供 SVM 実装の中には、ハイパーパラメータを自動的に、通常はランダムなサンプルに基づいて調整することができるものもある (Figure \@ref(fig:partitioning) の上段を参照)。 これは非空間データでは有効だが、空間データではあまり意味がなく、「空間チューニング」を行う必要がある。 空間チューニングを定義する前に、Section \@ref(glm) で紹介した **mlr3**\index{mlr3 (package)} ビルディングブロックを SVM 用に設定することにする。 -分類\index{classification}のタスクは変わらないので、Section \@ref(glm) で作成した `task` オブジェクトを再利用すればよい。 +分類\index{ぶんるい@分類}のタスクは変わらないので、Section \@ref(glm) で作成した `task` オブジェクトを再利用すればよい。 SVM を実装している学習器は、**mlr3extralearners** パッケージの `listLearners()` を用いて検索することができる。 ```{r 12-spatial-cv-23, eval=FALSE, echo=FALSE} @@ -502,10 +506,7 @@ lrn_ksvm$fallback = lrn("classif.featureless", predict_type = "prob") ``` 次の段階は、リサンプリング戦略を指定することである。 -ここでも 100 回繰り返しの 5 回空間 CV\index{cross-validation!spatial CV} を使用する。 - - +ここでも 100 回繰り返しの 5 回空間 CV\index{くろすばりでーしょん@クロスバリデーション!くうかん CV@空間 CV} を使用する。 ```{r 12-spatial-cv-25} # パフォーマンス推定レベル @@ -515,16 +516,16 @@ perf_level = mlr3::rsmp("repeated_spcv_coords", folds = 5, repeats = 100) これは Section \@ref(glm) の GLM \index{GLM} のリサンプリングに使われたコードと全く同じであることに注意しておこう。 ここまでは、Section \@ref(glm) で説明したものと同じである。 -しかし、次のステップは新しいもので、ハイパーパラメータ \index{hyperparameter} を調整することである。 +しかし、次のステップは新しいもので、ハイパーパラメータ \index{はいぱーぱらめーた@ハイパーパラメータ} を調整することである。 性能評価とチューニングに同じデータを使用すると、楽観的すぎる結果になる可能性がある [@cawley_overfitting_2010]。 -これは、ネストされた空間的な CV \index{cross-validation!spatial CV} を用いることで回避することができる。 +これは、ネストされた空間的な CV \index{くろすばりでーしょん@クロスバリデーション!くうかん CV@空間 CV} を用いることで回避することができる。 ```{r inner-outer, echo=FALSE, fig.cap="CV におけるハイパーパラメータのチューニングと性能推定レベルの模式図 (図は Schratz et al. (2019) から引用した。快く再利用の許可をいただいた)。", fig.scap="Schematic of hyperparameter tuning."} knitr::include_graphics("figures/12_cv.png") ``` -これは、各フォールドを空間的に不連続な 5 つのサブフォールドに再び分割し、最適なハイパーパラメータ \index{hyperparameter} ( `tune_level` 以下のコードチャンクのオブジェクト。視覚的表現については Figure \@ref(fig:inner-outer) を参照) を決定するために使用することを意味する。 -さらに、値 C と Sigma のランダムな選択は、あらかじめ定義された調整空間 ( `search_space` オブジェクト) に制限されている。 +これは、各フォールドを空間的に不連続な 5 つのサブフォールドに再び分割し、最適なハイパーパラメータ \index{はいぱーぱらめーた@ハイパーパラメータ} ( `tune_level` 以下のコードチャンクのオブジェクト。視覚的表現については Figure \@ref(fig:inner-outer) を参照) を決定するために使用することを意味する。 +さらに、値 C と Sigma のランダムな選択は、あらかじめ定義された調整空間 (`search_space` オブジェクト) に制限されている。 同調空間の範囲は、文献で推奨されている値で選択した [@schratz_hyperparameter_2019]。 最適なハイパーパラメータの組み合わせを見つけるために、これらのサブフォルダそれぞれにおいて、ハイパーパラメータ C とシグマにランダムな値を選択して 50 のモデル (以下のコードチャンクの `terminator` オブジェクト) を適合させた。 @@ -541,10 +542,10 @@ terminator = mlr3tuning::trm("evals", n_evals = 50) tuner = mlr3tuning::tnr("random_search") ``` -次の段階は、`AutoTuner$new()` を用いてハイパーパラメータチューニングを定義するすべての特性に従って学習器 `lrn_ksvm` を修正することである。 +次の段階は、`auto_tuner()` を用いてハイパーパラメータチューニングを定義するすべての特性に従って学習器 `lrn_ksvm` を修正することである。 ```{r 12-spatial-cv-27, eval=FALSE} -at_ksvm = mlr3tuning::AutoTuner$new( +at_ksvm = mlr3tuning::auto_tuner( learner = lrn_ksvm, resampling = tune_level, measure = mlr3::msr("classif.auc"), @@ -554,7 +555,7 @@ at_ksvm = mlr3tuning::AutoTuner$new( ) ``` -このチューニングは、1つのフォールドに対して最適なハイパーパラメータを決定するために、250 のモデルを適合させるように設定されている。 +このチューニングは、1 つのフォールドに対して最適なハイパーパラメータを決定するために、250 のモデルを適合させるように設定されている。 これを 1 回ずつ繰り返すと、1,250 個 (250 \* 5) のモデルができあがる。 100 回繰り返すということは、合計 125,000 個のモデルを適合して最適なハイパーパラメータ( Figure \@ref(fig:partitioning) )を特定することになる。 これを性能推定に使用し、さらに 500 個のモデル (5 folds \* 100 repetitions; Figure \@ref(fig:partitioning) 参照) の適合が必要である。 @@ -562,13 +563,13 @@ at_ksvm = mlr3tuning::AutoTuner$new( 1. パフォーマンスレベル ( Figure \@ref(fig:inner-outer) の左上部分) - データセットを5つの空間的に不連続な (外側の) サブフォールドに分割する。 1. チューニング・レベル ( Figure \@ref(fig:inner-outer) の左下部分) - パフォーマンス・レベルの最初のフォールドを使用し、ハイパーパラメータのチューニングのために、それを再び5つの (内側の) サブフォールドに空間的に分割する。 -これらの内部サブフォールドのそれぞれで、ランダムに選択された50個のハイパーパラメータ\index{hyperparameter} を使用する、つまり、250個のモデルを適合させる。 +これらの内部サブフォールドのそれぞれで、ランダムに選択された 50 個のハイパーパラメータ\index{はいぱーぱらめーた@ハイパーパラメータ} を使用する、つまり、250 個のモデルを適合させる。 1. 性能推定 - 前のステップ (チューニング・レベル) から最適なハイパーパラメータの組み合わせを使用し、性能レベルの最初の外側のフォールドに適用して性能を推定する (AUROC \index{AUROC} )。 -1. 残りの4つの外側の折り目について、手順2と3を繰り返す -1. 手順2～4を100回繰り返す +1. 残りの 4 つの外側の折り目について、手順 2 と 3 を繰り返す +1. 手順 2～4 を 100 回繰り返す ハイパーパラメータのチューニングと性能推定のプロセスには、計算量が必要である。 -モデルの実行時間を短縮するために、**mlr3** では、**future** パッケージの助けを借りて、並列化\index{parallelization}を使用する可能性を提供している。 +モデルの実行時間を短縮するために、**mlr3** では、**future** パッケージの助けを借りて、並列化\index{へいれつか@並列化}を使用する可能性を提供している。 これからネストした CV を実行するので、内側ループと外側ループのどちらを並列化するか決めることができる (Figure \@ref(fig:inner-outer) の左下部分を参照)。 前者は 125,000 個のモデルを実行するのに対し、後者は 500 個しか実行しないので、内側のループを並列化するのは当然である。 内側のループの並列化を設定するために、実行する。 @@ -584,11 +585,11 @@ future::plan(list("sequential", "multisession"), これで、ネストされた空間 CV を計算するための準備ができた。 `resample()` パラメータの指定は、GLM\index{GLM} を使用したときと全く同じ手順で行う。唯一の違いは、`store_models` と `encapsulate` の引数である。 -前者を `TRUE` に設定すると、ハイパーパラメータ\index{hyperparameter}のチューニング結果を抽出できる。これは、チューニングに関するフォローアップ分析を計画する場合に重要である。 +前者を `TRUE` に設定すると、ハイパーパラメータ\index{はいぱーぱらめーた@ハイパーパラメータ}のチューニング結果を抽出できる。これは、チューニングに関するフォローアップ分析を計画する場合に重要である。 後者は、モデルの 1 つがエラーを投げても処理が継続されるようにするものである。 これにより、1 つのモデルが失敗しただけで処理が停止することを避けることができ、大規模なモデルの実行には望ましい。 処理が完了すると、故障したモデルを見ることができる。 -処理終了後、`future::ClusterRegistry("stop")` で明示的に並列化\index{parallelization}を停止するのがよいだろう。 +処理終了後、`future::ClusterRegistry("stop")` で明示的に並列化\index{へいれつか@並列化}を停止するのがよいだろう。 最後に、出力オブジェクト (`result`) を、別の R セッションで使用する場合に備えてディスクに保存する。 125,500 個のモデルで空間クロスバリデーションを行うため、時間がかかることをご了承の上、実行してみよう。 現代のコンピュータで、たった半日で実行できる。 @@ -628,39 +629,39 @@ score_spcv_svm = score[learner_id == "classif.ksvm.tuned" & round(mean(score_spcv_svm$classif.auc), 2) ``` -GLM\index{GLM} (集計された AUROC \index{AUROC} は `r score [score $resampling_id == "repeated_spcv_coords" & score$ learner_id == "classif.log_reg", round(mean(score$classif.auc), 2)] `) は、この特定のケースでは、SVM\index{SVM} よりもわずかに優れているようである。 +GLM\index{GLM} (集計された AUROC\index{AUROC} は `r score [score $resampling_id == "repeated_spcv_coords" & score$ learner_id == "classif.log_reg", round(mean(score$classif.auc), 2)] `) は、この特定のケースでは、SVM\index{SVM} よりもわずかに優れているようである。 絶対的に公平な比較を保証するためには、2 つのモデルが全く同じパーティションを使用していることを確認する必要がある。ここでは示していないが、バックグラウンドで黙々と使用しているものである (詳しくは本書の GitHub リポジトリにある `code/12_cv.R` を参照)。 そのために、**mlr3** は関数 `benchmark_grid()` と `benchmark()` を提供している [https://mlr3book.mlr-org.com/chapters/chapter3/evaluation_and_benchmarking.html#sec-benchmarking, @becker_mlr3_2022 参照] 。 これらの機能については、「演習」でより詳しく解説する。 -また、SVM のランダムな探索に 50 回以上の反復を使用すると、おそらくより良い AUROC [@schratz_hyperparameter_2019] を持つモデルになるハイパーパラメータ\index{hyperparameter}が得られるであろうことに注意しておこう。 +また、SVM のランダムな探索に 50 回以上の反復を使用すると、おそらくより良い AUROC [@schratz_hyperparameter_2019] を持つモデルになるハイパーパラメータ\index{はいぱーぱらめーた@ハイパーパラメータ}が得られるであろうことに注意しておこう。 一方、ランダムサーチの反復回数を増やすと、モデルの総数も増え、その分実行時間も長くなる。 -これまで、空間 CV \index{cross-validation!spatial CV} は、学習アルゴリズムが未知のデータに対して汎化する能力を評価するために利用されていた。 -空間分布図作成では、完全なデータセットでハイパーパラメータ\index{hyperparameter}を調整する。 +これまで、空間 CV \index{くろすばりでーしょん@クロスバリデーション!くうかん CV@空間 CV} は、学習アルゴリズムが未知のデータに対して汎化する能力を評価するために利用されていた。 +空間分布図作成では、完全なデータセットでハイパーパラメータ\index{はいぱーぱらめーた@ハイパーパラメータ}を調整する。 これについては、Chapter \@ref(eco) で説明する。 ## 結論 {#conclusions} -リサンプリング手法は、データサイエンティストのツールボックスの重要な一部である [@james_introduction_2013]。 -この章では、CV\index{cross-validation} を用いて、様々なモデルの予測性能を評価した。 -Section \@ref(intro-cv) で述べたように、空間座標を持つオブザベーションは、空間自己相関\index{autocorrelation!spatial}のために統計的に独立でない場合があり、クロスバリデーションの基本的な仮定に違反する。 -空間 CV \index{cross-validation!spatial CV} この問題は、空間的自己相関\index{autocorrelation!spatial}によってもたらされるバイアスを低減することで解決される。 +リサンプリング手法は、データサイエンティストのツールボックスの重要なものの一つである [@james_introduction_2013]。 +この章では、CV\index{くろすばりでーしょん@クロスバリデーション} を用いて、様々なモデルの予測性能を評価した。 +Section \@ref(intro-cv) で述べたように、空間座標を持つオブザベーションは、空間自己相関\index{じこそうかん@自己相関!くうかんてき@空間的}のために統計的に独立でない場合があり、クロスバリデーションの基本的な仮定に違反する。 +空間 CV\index{くろすばりでーしょん@クロスバリデーション!くうかん CV@空間 CV} この問題は、空間的自己相関\index{じこそうかん@自己相関!くうかんてき@空間的}によってもたらされるバイアスを低減することで解決される。 -**mlr3**\index{mlr3 (package)} パッケージは、線形回帰\index{regression!linear}、一般化加法モデル\index{generalized additive model}などのセミパラメトリックモデルなどの統計学習\index{statistical learning}、あるいはランダムフォレスト\index{random forest}、SVM \index{SVM}、ブースト回帰木 [@bischl_mlr:_2016;@schratz_hyperparameter_2019] などの機械学習\index{machine learning} 技術と組み合わせることで、 (空間) リサンプリング技法\index{resampling}を容易にしている。 -機械学習アルゴリズムは、ハイパーパラメータ\index{hyperparameter}の入力を必要とすることがある。その最適な「チューニング」には、大規模な計算資源を必要とする数千回のモデル実行が必要で、多くの時間、RAM、コアを消費することがある。 -**mlr3** は、並列化\index{parallelization}を可能にすることでこの問題に取り組んでいる。 +**mlr3**\index{mlr3 (package)} パッケージは、線形回帰\index{かいき@回帰!せんけい@線形}、一般化加法モデル\index{いっぱんかかほうもでる@一般化加法モデル}などのセミパラメトリックモデルなどの統計学習\index{とうけいてきがくしゅう@統計的学習}、あるいはランダムフォレスト\index{らんだむふぉれすと@ランダムフォレスト}、SVM \index{SVM}、ブースト回帰木 [@bischl_mlr:_2016;@schratz_hyperparameter_2019] などの機械学習\index{きかいがくしゅう@機械学習} 技術と組み合わせることで、 (空間) リサンプリング技法\index{resampling}を容易にしている。 +機械学習アルゴリズムは、ハイパーパラメータ\index{はいぱーぱらめーた@ハイパーパラメータ}の入力を必要とすることがある。その最適な「チューニング」には、大規模な計算資源を必要とする数千回のモデル実行が必要で、多くの時間、RAM、コアを消費することがある。 +**mlr3** は、並列化\index{へいれつか@並列化}を可能にすることでこの問題に取り組んでいる。 機械学習全体、そして空間データを理解するための機械学習は大きな分野であり、この章では基本的なことを説明したが、まだまだ学ぶべきことはある。 このような方向性で、以下の資料を勧める。 - **mlr3 book** (@becker_mlr3_2022; https://mlr3book.mlr-org.com/) と、特に [chapter on the handling of spatio-temporal data](https://mlr3book.mlr-org.com/chapters/chapter13/beyond_regression_and_classification.html#sec-spatiotemporal) を参考。 -- ハイパーパラメータ\index{hyperparameter}チューニングに関する学術論文 [@schratz_hyperparameter_2019] +- ハイパーパラメータ\index{はいぱーぱらめーた@ハイパーパラメータ}チューニングに関する学術論文 [@schratz_hyperparameter_2019] - **mlr3spatiotempcv** の使用方法に関する学術論文 [@schratz_mlr3spatiotempcv_2021] -- 時空間データの場合、空間的\index{autocorrelation!spatial}と時間的\index{autocorrelation!temporal}の自己相関を考慮した上で CV\index{cross-validation} を行う必要がある [@meyer_improving_2018]。 +- 時空間データの場合、空間的\index{autocorrelation!spatial}と時間的\index{autocorrelation!temporal}の自己相関を考慮した上で CV\index{くろすばりでーしょん@クロスバリデーション} を行う必要がある [@meyer_improving_2018]。 ## 演習 ```{r, echo=FALSE, results='asis'} -res = knitr::knit_child('_12-ex.Rmd', quiet = TRUE, options = list(include = FALSE, eval = FALSE)) +res = knitr::knit_child('_12-ex-ja.Rmd', quiet = TRUE, options = list(include = FALSE, eval = FALSE)) cat(res, sep = '\n') ``` diff --git a/13-transport-ja.Rmd b/13-transport-ja.Rmd index 5bb8f27..03d6760 100644 --- a/13-transport-ja.Rmd +++ b/13-transport-ja.Rmd @@ -33,7 +33,7 @@ library(sfnetworks) # 空間ネットワークのクラスと関数 > Everything is related to everything else, but near things are more related than distant things. -この「法則」は、空間自己相関 \index{autocorrelation!spatial} など、地理学の重要な概念の基礎となるものである。 +この「法則」は、空間自己相関\index{じこそうかん@自己相関!くうかんてき@空間的}など、地理学の重要な概念の基礎となるものである。 それは、友情のネットワークや生態系の多様性といった多様な現象に適用され、「距離の摩擦」を構成する時間、エネルギー、金銭といった輸送コストによって説明することができる。 この観点からすると、交通テクノロジーは破壊的であり、移動する人間やモノを含む地理的な実体間の関係を変化させる。「交通の目的は空間を克服することである」 [@rodrigue_geography_2013]。 @@ -42,20 +42,20 @@ library(sfnetworks) # 空間ネットワークのクラスと関数 本章では、交通システムの地理的分析について、さまざまな地理的レベルで紹介する。 -- **実際の単位**：交通パターンは、主な移動手段 (車、自転車、徒歩など) や、特定のゾーンに住む人々の平均移動距離など、ゾーンごとの集計を参照して理解することができる ( Section \@ref(transport-zones) で解説する)。 -- **希望線** (desire line)\index{desire lines}: 地理空間における場所 (点またはゾーン) 間を何人が移動したか (または移動できたか) を記録した「起点-終点」データを表す直線で、Section \@ref(desire-lines) の話題である。 -- **ノード** (node)\index{node} : ノードとは、共通の起点と終点を表すことができる交通システムの点であり、バス停や鉄道駅などの公共交通機関の駅、Section \@ref(nodes) のトピックである。 -- **ルート**：希望線に沿ったルート網とノード間のルートを表す線である。 +- **面積単位**\index{めんせきたんい@面積単位}: 交通パターンは、主な移動手段 (車、自転車、徒歩など) や、特定のゾーンに住む人々の平均移動距離など、ゾーンごとの集計を参照して理解することができる (Section \@ref(transport-zones) で解説する)。 +- **希望線** (desire line)\index{きぼうせん@希望線}: 地理空間における場所 (点またはゾーン) 間を何人が移動したか (または移動できたか) を記録した「起点-終点」データを表す直線で、Section \@ref(desire-lines) の話題である。 +- **ノード** (node)\index{のーど@ノード} : ノードとは、共通の起点と終点を表すことができる交通システムの点であり、バス停や鉄道駅などの公共交通機関の駅、Section \@ref(nodes) のトピックである。 +- **ルート**\index{るーと@ルート}: 希望線に沿ったルート網とノード間のルートを表す線である。 ルート (１つのラインストリングまたは複数のセグメントの場合がある) と、ルートを生成するルートエンジンは、Section \@ref(routes) で詳しく述べる。 -- **ルートネットワーク** (network)\index{network} : これらは、ある地域の道路、小道、その他の線形特徴のシステムを表し、Section \@ref(route-networks) で取り上げている。 +- **ルートネットワーク** (network)\index{ねっとわーく@ネットワーク}: これらは、ある地域の道路、小道、その他の線形特徴のシステムを表し、Section \@ref(route-networks) で取り上げている。 これらは地理的な特徴 (短いセグメントでネットワーク全体を作り上げる) として表現することもできるし、相互接続されたグラフとして構造化することもできる。異なるセグメント上の交通のレベルは、交通モデルでは「フロー」と呼ばれる [@hollander_transport_2016]。 もうひとつの重要なレベルは、私やあなたのような移動する存在を表す**エージェント**である。 -エージェントは、 [MATSim](http://www.matsim.org/) のようなソフトウェアのおかげで、計算によって表現することができる。これは、エージェントベースモデリング (ABM)\index{agent-based modeling} のアプローチを用いて、高い空間および時間分解能で輸送システムのダイナミクスを捉えるものである [@horni_multi-agent_2016]。 +エージェントは、 [MATSim](http://www.matsim.org/) のようなソフトウェアのおかげで、計算によって表現することができる。これは、エージェントベースモデリング (ABM)\index{えーじぇんとべーすもでりんぐ@エージェントベースモデリング}のアプローチを用いて、高い空間および時間分解能で輸送システムのダイナミクスを捉えるものである [@horni_multi-agent_2016]。 ABM は、R の空間クラスと統合できる可能性が高い交通研究の強力なアプローチであるが [@thiele_r_2014; @lovelace_spatial_2016]、本章の範囲外である。 地理的レベルやエージェントの次に来るのは、ほとんどの交通モデルの分析の基本単位である**トリップ**であり、出発地「A」から目的地「B」までの単一目的の旅である [@hollander_transport_2016]。 -トリップは、異なるレベルの交通システムを結合し、単純化すると、ゾーンの重心\index{centroid} (ノード) を結ぶ地理的な希望線として、あるいは交通ルートネットワークに沿った経路として表現することができる。 -この文脈では、エージェント\index{agent-based modeling}は通常、交通ネットワーク内を移動する点である。 +トリップは、異なるレベルの交通システムを結合し、単純化すると、ゾーンの重心\index{じゅうしん@重心} (ノード) を結ぶ地理的な希望線として、あるいは交通ルートネットワークに沿った経路として表現することができる。 +この文脈では、エージェント\index{えーじぇんとべーすもでりんぐ@エージェントベースモデリング}は通常、交通ネットワーク内を移動する点である。 交通システムは動的である [@xie_evolving_2011]。 この章では、交通システムの地理的な分析に焦点を当てるが、変化のシナリオをシミュレートするために、この手法をどのように利用できるかについて、Section \@ref(prioritizing-new-infrastructure) にて詳しくみていこう。 @@ -65,7 +65,7 @@ ABM は、R の空間クラスと統合できる可能性が高い交通研究 一般的に、モデルは特定の問題を解決するために設計される。 そのため、本章では、次章で紹介する政策シナリオを軸に、問いかけを行っている。ブリストル市内で自転車を増やす方法は？ Chapter \@ref(location) では、ジオコンピュテーションの応用として、新しい自転車店の立地の優先順位付けを行う。 -新しいサイクリング・インフラが効果的に配置されると、人々がサイクリングするようになり、サイクリング・ショップの需要や地域の経済活動を高めることができる。 +新しい自転車利用・インフラが効果的に配置されると、人々が自転車利用するようになり、自転車利用・ショップの需要や地域の経済活動を高めることができる。 これは、交通システムの重要な特徴を強調している。交通システムは、より広範な現象や土地利用パターンと密接に関連している。 ## ブリストルのケーススタディ {#bris-case} @@ -93,13 +93,13 @@ tm_shape(region_all[1, ], bbox = region_all) + tm_basemap(server = leaflet::providers$Esri.WorldTopoMap) ``` -```{r bristol, echo=FALSE, fig.cap="ブリストルの交通網は、アクティブ (緑)、公共 (鉄道、黒)、自家用車 (赤) の各移動手段を色分けして表現している。黒い境界線は、市街地の境界線 (黄色でハイライト) と、より広い TTWA (Travel To Work Area) を表している。", fig.scap="Bristol's transport network."} +```{r bristol, echo=FALSE, fig.cap="ブリストルの交通網は、アクティブ (緑)、公共 (鉄道、青)、自家用車 (赤) の各移動手段を色分けして表現している。黒い境界線は、市街地の境界線 (黄色でハイライト) と、より広い TTWA (Travel To Work Area) を表している。", fig.scap="Bristol's transport network."} knitr::include_graphics("figures/13_bristol.png") ``` -ブリストルはイングランドで 10 番目に大きい市で、人口は 50 万人であるが、その旅行キャッチメントエリア\index{catchment area}はもっと大きい (Section \@ref(transport-zones) 参照)。 +ブリストルはイングランドで 10 番目に大きい市で、人口は 50 万人であるが、その旅行キャッチメントエリア\index{きゃっちめんとえりあ@キャッチメントエリア}はもっと大きい (Section \@ref(transport-zones) 参照)。 市内には航空宇宙、メディア、金融サービス、観光などの企業が集まり、2 つの主要な大学とともに、活気ある経済が展開されている。 -ブリストルは一人当たりの平均所得が高いが、深刻な貧困地域も含まれている [@bristol_city_council_deprivation_2015]。 +Bristol は一人当たりの平均所得が高いが、深刻な貧困地域も含まれている [@bristol_city_council_deprivation_2015]。 交通の面では、ブリストルは鉄道や道路の便がよく、アクティブ・トラベルのレベルも比較的高い。 [Active People Survey](https://www.gov.uk/government/statistical-data-sets/how-often-and-time-spent-walking-and-cycling-at-local-authority-level-cw010#table-cw0103) によると、国民の19％が月に1回以上自転車を利用し、88％が歩いている (全国平均はそれぞれ15％、81％)。 @@ -119,15 +119,15 @@ View(cw0103) 多くの都市と同様にブリストルも渋滞、大気質、運動不足などの大きな問題を抱えている。 自転車は、これらの問題すべてに効率的に取り組むことができる。典型的な[速度](https://en.wikipedia.org/wiki/Bicycle_performance)は、徒歩の時速4〜6kmに対して時速15〜20kmと、徒歩よりも自動車による移動を置き換える可能性が大きい。 -このため、ブリストルの[交通戦略](https://www.bristol.gov.uk/council-and-mayor/policies-plans-and-strategies/bristol-transport-strategy)では、サイクリングについて野心的な計画を立てている。 +このため、ブリストルの[交通戦略](https://www.bristol.gov.uk/council-and-mayor/policies-plans-and-strategies/bristol-transport-strategy)では、自転車利用について野心的な計画を立てている。 この章では、交通研究における政策的配慮の重要性を強調するため、人々を車から解放し、より持続可能な手段、特に徒歩と自転車に乗せることを任務とする人々 (交通プランナー、政治家、その他の利害関係者) にエビデンスを提供する目的で書かれている。 より広い目的として、ジオコンピューティングがどのように証拠に基づく交通計画をサポートできるかを示すことである。 この章では、次のことを学ぶ。 - 都市における交通行動の地理的パターンを説明する -- マルチ・モデル・トリップをサポートする主要な公共交通機関のノードを特定する。 -- 移動の「希望線」\index{desire lines}を分析、多くの人が短距離をドライブする場所を見つける +- マルチモーダル・トリップをサポートする主要な公共交通機関のノードを特定する。 +- 移動の「希望線」\index{きぼうせん@希望線}を分析、多くの人が短距離をドライブする場所を見つける - 自動車を減らし自転車を増やすような自転車ルートの位置を特定する 本章の実用的な面から始めるため、次節で、移動パターンに関するゾーンデータをロードする。 @@ -135,7 +135,7 @@ View(cw0103) ## 交通ゾーン {#transport-zones} -交通システムは主に線形のフィーチャやノード\index{node} (例えば、ルートや駅) に基づいているが、連続した空間を具体的な単位に分割するために、面的なデータから始めることがしばしば意味を持つ [@hollander_transport_2016]。 +交通システムは主に線形のフィーチャやノード\index{のーど@ノード} (例えば、ルートや駅) に基づいているが、連続した空間を具体的な単位に分割するために、面的なデータから始めることがしばしば意味を持つ [@hollander_transport_2016]。 調査地域 (ここではブリストル) を定義する境界線に加えて、交通研究者にとって特に関心の高い2つのゾーンタイプ、すなわち発着地ゾーンがある。 多くの場合、発着地には同じ地理的単位が使用される。 しかし、学校や商店などの「トリップアトラクター」が多い地域では、「 [Workplace Zones](https://www.data.gov.uk/dataset/022e455a-b5cc-4247-bfc5-1aaf9da78379/workplace-zones-a-new-geography-for-workplace-statistics)」のような異なるゾーニングシステムが、トリップ先の密度上昇を表すのに適切だろう [@office_for_national_statistics_workplace_2014]。 @@ -177,7 +177,7 @@ TTWA はまず、人口の75％が通勤で利用する連続したゾーンと names(bristol_zones) ``` -旅行データを追加するために、Section \@ref(vector-attribute-joining) で説明されている一般的なタスクである *attribute join* \index{attribute!join} を実行する。 +旅行データを追加するために、Section \@ref(vector-attribute-joining) で説明されている一般的なタスクである *attribute join*\index{ぞくせい@属性!けつごう@結合} を実行する。 ここでは、英国の2011年国勢調査の出勤時間に関する質問から、 [ons.gov.uk](https://www.ons.gov.uk/help/localstatistics) データポータルから提供された、`bristol_od` に保存されている旅行データを使用する。 `bristol_od` は、英国の2011年国勢調査によるゾーン間の通勤に関する起点 (Origin)ー終点 (Desitination) (OD) データセットである (Section \@ref(desire-lines) 参照)。 Section \@ref(desire-lines) 1列目は出発地のゾーンID、2列目は目的地のゾーンである。 @@ -248,9 +248,6 @@ zones_od = inner_join(zones_joined, zones_destinations, by = "geo_code") Figure \@ref(fig:zones) の簡易版は、以下のコードで作成する (図を再現するには、本書の GitHub リポジトリの [`code`](https://github.com/Robinlovelace/geocompr/tree/main/code) フォルダ の `13-zones.R` を参照。**tmap**\index{tmap (package)} によるファセット・マップの詳細は Section \@ref(faceted-maps) を参照)。 - - - ```{r 13-transport-11, eval=FALSE} qtm(zones_od, c("all", "all_dest")) + tm_layout(panel.labels = c("Origin", "Destination")) @@ -310,7 +307,7 @@ od_inter = filter(bristol_od, o != d) 次のステップは、ゾーン間 OD ペアを、**stplanr**\index{stplanr (package)} 関数 `od2line()` を用いて地図上にプロットできる希望線を表す `sf` オブジェクトに変換することである。[^13-transport-5] [^13-transport-5]:`od2line()` は、 オブジェクトの最初の 2 列の ID を、地理的な オブジェクトの ID 列にマッチングさせることで動作する。 - `bristol_od` `zones_od` `zone_code` なお、この操作は警告を発するが、これは `od2line()` が各起点と終点のペアをその起点と終点のゾーンの重心\index{centroid}に割り当てることで動作するためである。 + `bristol_od` `zones_od` `zone_code` なお、この操作は警告を発するが、これは `od2line()` が各起点と終点のペアをその起点と終点のゾーンの重心\index{じゅうしん@重心}に割り当てることで動作するためである。 実際の使用では、投影データから生成された重心を使用するか、できれば人口重み付けされた重心を使用することになるだろう [@lovelace_propensity_2017]。 ```{r 13-transport-15, warning=FALSE} @@ -334,24 +331,24 @@ source("code/13-desire.R", print.eval = TRUE) ## ノード {#nodes} -地理的な交通データにおけるノード\index{node}は、ネットワーク\index{network}を構成する主に一次元のフィーチャ (線) のうち、ゼロ次元のフィーチャ (点) である。 +地理的な交通データにおけるノード\index{のーど@ノード}は、ネットワーク\index{ねっとわーく@ネットワーク}を構成する主に一次元のフィーチャ (線) のうち、ゼロ次元のフィーチャ (点) である。 交通ノードには2種類ある。 -1. ネットワーク\index{network}上に直接存在しないノード\index{node}、例えばゾーン重心\index{centroid} (次のセクションで取り上げる) あるいは家や職場などの個人の発着地。 +1. ネットワーク\index{ねっとわーく@ネットワーク}上に直接存在しないノード\index{のーど@ノード}、例えばゾーン重心\index{じゅうしん@重心} (次のセクションで取り上げる) あるいは家や職場などの個人の発着地。 2. 交通網の一部であるノード。 技術的には、ノードは交通ネットワーク上のどの点にも位置することができるが、実際には、ルートの交差点 (ジャンクション) やバス停や駅など交通ネットワークに出入りする点など、特殊な頂点である場合が多い。[^13-transport-6] [^13-transport-6]: **sfnetworks** パッケージの関数 [`st_network_blend()`](https://luukvdmeer.github.io/sfnetworks/reference/st_network_blend.html) により、ネットワーク上またはネットワーク外の任意の点への近接性に基づいて、ネットワーク上に新しいノードを作成することができる。 -交通ネットワークは、グラフ\index{graph}として表すことができ、その中で各セグメントは、ネットワーク内の1つ以上の他のエッジ\index{edge}に (地理的な線を表すエッジを介して) 接続されている。 +交通ネットワークは、グラフ\index{ぐらふ@グラフ}として表すことができ、その中で各セグメントは、ネットワーク内の 1 つ以上の他のエッジ\index{えっじ@エッジ}に (地理的な線を表すエッジを介して) 接続されている。 ネットワーク外のノードは「重心コネクタ」で追加可能。重心コネクタとは、ネットワーク上の近くのノードへの新しいルートセグメントである [@hollander_transport_2016] 。[^13-transport-7] -ネットワークの各ノード\index{node}は、ネットワーク上の個々のセグメントを表す 1 つ以上の「エッジ」によって接続されている。 +ネットワークの各ノード\index{のーど@ノード}は、ネットワーク上の個々のセグメントを表す 1 つ以上の「エッジ」によって接続されている。 交通ネットワークがグラフとして表現できることを Section \@ref(route-networks) で確認する。 [^13-transport-7]: これらのコネクターは、その近辺の交通量を過大に見積もる可能性があるため、その位置は慎重に選ぶ必要がある [@jafari_investigation_2015]。 -公共交通機関の駅や停留所は特に重要なノード\index{node}である。道路の一部であるバス停、または線路から数百メートルの歩行者入口ポイントによって表される大規模な鉄道駅のいずれかのタイプのノードとして表現されることができる。 -ここでは、ブリストルにおける自転車の増加という研究課題に関連して、公共交通機関のノード\index{node}を説明するために鉄道駅を使用しよう。 +公共交通機関の駅や停留所は特に重要なノード\index{のーど@ノード}である。道路の一部であるバス停、または線路から数百メートルの歩行者入口ポイントによって表される大規模な鉄道駅のいずれかのタイプのノードとして表現されることができる。 +ここでは、ブリストルにおける自転車の増加という研究課題に関連して、公共交通機関のノード\index{のーど@ノード}を説明するために鉄道駅を使用しよう。 鉄道駅のデータは、`bristol_stations` の **spDataLarge** で提供されている。 通勤において自動車から他の交通への転換を阻む共通の障壁は、自宅から職場までの距離が遠すぎると徒歩や自転車では無理だということである。 @@ -363,16 +360,16 @@ source("code/13-desire.R", print.eval = TRUE) - 降車駅から目的地までの行程。 Section \@ref(desire-lines) で行った分析に基づき、公共交通機関のノードを使って、バスと (この例では) 鉄道を利用できる旅行のための3分割の希望線を構築することができる。 -最初の段階は、公共交通機関の利用が多い希望線を特定することである。ここでは、先に作成したデータセット `desire_lines` にすでに電車での移動回数を表す変数が含まれているので簡単である (公共交通機関の利用可能性は、[OpenTripPlanner](http://www.opentripplanner.org/) などの公共交通ルート探索\index{routing}サービスを使って推定することも可能)。 -アプローチを簡単にするために、レールの使用量の上位3つの希望線\index{desire lines}のみを選択することにする。 +最初の段階は、公共交通機関の利用が多い希望線を特定することである。ここでは、先に作成したデータセット `desire_lines` にすでに電車での移動回数を表す変数が含まれているので簡単である (公共交通機関の利用可能性は、[OpenTripPlanner](http://www.opentripplanner.org/) などの公共交通ルート探索\index{るーとけんさく@ルート検索}サービスを使って推定することも可能)。 +アプローチを簡単にするために、レールの使用量の上位 3 つの希望線\index{きぼうせん@希望線}のみを選択することにする。 ```{r 13-transport-20} desire_rail = top_n(desire_lines, n = 3, wt = train) ``` そこで、これらの線を 3 つに分解し、公共交通機関の乗り換えを表現することにする。 -これは、希望線を、旅行の出発地、公共交通機関、目的地の3つの線ジオメトリからなるマルチラインストリングに変換することで実現できる。 -この作業は、行列の作成 (起点、終点、鉄道駅を表す「経由点」)、最近傍\index{nearest neighbor}の特定、マルチラインストリング\index{multilinestrings}への変換の3段階に分けることができる。 +これは、希望線を、旅行の出発地、公共交通機関、目的地の 3 つの線ジオメトリからなる複合線に変換することで実現できる。 +この作業は、行列の作成 (起点、終点、鉄道駅を表す「経由点」)、最近傍\index{さいきんぼう@最近傍}の特定、複合線\index{ふくごうせん@複合線}への変換の 3 段階に分けることができる。 この一連の処理は、`line_via()` 関数が行う。 この **stplanr**\index{stplanr (package)} 関数は入力された線と点を受け取り、希望線のコピーを返す (この動作の詳細については geocompr.github.io ウェブサイトの [Desire Lines Extended](https://geocompr.github.io/geocompkg/articles/linevia.html) vignette と `?line_via` を参照)。 出力は入力と同じだが、公共交通機関のノードを使った旅を表す新しいジオメトリの列がある (以下に示す)。 @@ -383,7 +380,7 @@ desire_rail = line_via(desire_rail, bristol_stations) ncol(desire_rail) ``` -Figure \@ref(fig:stations) に示すように、最初の `desire_rail` の行に、自宅から出発駅まで、そこから目的地まで、そして最後に目的地から目的地までの移動を表す 3 つのジオメトリリスト列\index{list column}が追加されている。 +Figure \@ref(fig:stations) に示すように、最初の `desire_rail` の行に、自宅から出発駅まで、そこから目的地まで、そして最後に目的地から目的地までの移動を表す 3 つのジオメトリリスト列\index{りすとれつ@リスト列}が追加されている。 この場合、目的地までの距離は非常に短いが (歩行距離)、出発地までの距離は十分であるため、往路の駅まで自転車での移動を促すためのインフラ投資を正当化することができる。Figure \@ref(fig:stations) にある 3 つの出発地の駅周辺の住宅地では、人々が通勤するために自転車を利用することができる。 ```{r stations, echo=FALSE, message=FALSE, warning=FALSE, fig.cap="鉄道利用率の高い直線的な希望線 (黒) を、公共交通機関 (グレー) を経由して出発駅 (赤) へ、そして目的地 (ごく短い青線) へという 3 レグに変換する中間点として使用される駅ノード (赤い点)。", fig.scap="Station nodes."} @@ -427,7 +424,8 @@ tm_shape(zones_od) + ## ルート {#routes} -地理学者から見れば、ルートとは直線でなくなった希望線\index{desire lines}である。出発地と目的地は同じだが、A から B へのルートはより複雑である。 +\index{るーと@ルート} +地理学者から見れば、ルートとは直線でなくなった希望線\index{きぼうせん@希望線}である。出発地と目的地は同じだが、A から B へのルートはより複雑である。 ルートは、ローカルまたはリモートで実行されるルート探索サービスを使用して、希望線 (より一般的には起点と終点のペア) から生成される。 希望線が2つの頂点 (始点と終点) しか持たないのに対し、ルートは任意の数の頂点を持つことができ、A-B間の点を直線で結んだものと定義される。これは線ジオメトリ (linestring) の定義である。 @@ -458,9 +456,9 @@ R5 のようなマルチモーダルなルーティングエンジンは、OSRM ルーティングエンジンは、ルート、レッグ、セグメントという 3 つの地理的なレベルで出力を生成することができる。 -- **Route** レベルの出力には、出発地と目的地のペアごとに 1 つのフィーチャー (通常はデータフレーム表現におけるマルチライン文字列と関連する行) が含まれ、トリップごとに 1 つのデータ行があることを意味する -- **Leg** レベルの出力には、各起点と終点のペアに 1 つのフィーチャーと関連する属性が含まれまる。1 つのモードを含むだけのトリップの場合 (たとえば、自宅から職場まで車で行き、車まで の短い徒歩は無視)、レッグはルートと同じで、車の旅になる。公共交通機関を利用するトリップの場合、レッグは重要な情報を提供する。**r5r** の関数 `detailed_itineraries()` はレッグを返すが、紛らわしいことに、これは「セグメント」と呼ばれることもある -- Segment レベルの出力は、交通ネットワークの各小セクションのレコードを持つ、ルートに関する最も詳細な情報を提供する。一般的にセグメントは、OpenStreetMap の道と同じか、同じ長さになっている。**cyclestreets** 関数 `journey()` はセグメントレベルのデータを返し、これは **stplanr** の `route()` 関数が返す出発地と目的地レベルのデータでグループ化することによって集約できる +- **ルート** レベルの出力には、出発地と目的地のペアごとに 1 つのフィーチャ (通常はデータフレーム表現におけるマルチライン文字列と関連する行) が含まれ、トリップごとに 1 つのデータ行があることを意味する +- **レッグ** レベルの出力には、各起点と終点のペアに 1 つのフィーチャと関連する属性が含まれまる。1 つのモードを含むだけのトリップの場合 (たとえば、自宅から職場まで車で行き、車まで の短い徒歩は無視)、レッグはルートと同じで、車の旅になる。公共交通機関を利用するトリップの場合、レッグは重要な情報を提供する。**r5r** の関数 `detailed_itineraries()` はレッグを返すが、紛らわしいことに、これは「セグメント」と呼ばれることもある +- セグメントレベルの出力は、交通ネットワークの各小セクションのレコードを持つ、ルートに関する最も詳細な情報を提供する。一般的にセグメントは、OpenStreetMap の道と同じか、同じ長さになっている。**cyclestreets** 関数 `journey()` はセグメントレベルのデータを返し、これは **stplanr** の `route()` 関数が返す出発地と目的地レベルのデータでグループ化することによって集約できる ほとんどのルーティングエンジンは、デフォルトでルートレベルを返すが、マルチモーダルエンジンは一般的にレッグレベル (単一の輸送モードによる連続した移動ごとに1つの機能) の出力を提供する。 セグメントレベルの出力は、より詳細な情報を提供するという利点がある。 @@ -483,14 +481,14 @@ R ネイティブなルーティングは高速で柔軟な反面、現実的な ### ローカル型専用ルーティングエンジン {#localengine} **ローカル型**ルーティングエンジンには、OpenTripPlanner、[Valhalla](https://github.com/valhalla/valhalla)、OpenStreetMap Routing Machine (OSRM) (これは「ユニモーダル」のみ対応)、R5 などがある。 -R からこうしたサービスにアクセスするには、**opentripplanner**、[**valhalla**](https://github.com/chris31415926535/valhallr)、**r5r**、[**osrm**](https://github.com/riatelab/osrm) などのパッケージがある [@morgan_opentripplanner_2019; @pereira_r5r_2021]。 +R からこうしたサービスにアクセスするには、**opentripplanner**、[**valhallr**](https://github.com/chris31415926535/valhallr)、**r5r**、[**osrm**](https://github.com/riatelab/osrm) などのパッケージがある [@morgan_opentripplanner_2019; @pereira_r5r_2021]。 ローカルにホストされたルーティングエンジンは、ユーザのコンピュータ上で実行されるが、R とは別のプロセスで実行される。 利点としては、実行速度が速く、異なる輸送手段に対する重み付けプロファイルを制御できるという点がある。 反対に欠点は、複雑なネットワークをローカルに表現することが難しく、時間的ダイナミクス (主にトラフィックによる)、特殊なソフトウェアが必要となるなどの点が挙げられる。 ### リモート型専用ルーティングエンジン {#remoteengine} -**リモート型**\index{routing}ルーティングエンジンは、Web API\index{API} を使用して、起点と終点に関するクエリーを送信し、専用のソフトウェアが動作する強力なサーバーで生成された結果を返す。 +**リモート型**\index{るーとけんさく@ルート検索}ルーティングエンジンは、Web API\index{API} を使用して、起点と終点に関するクエリーを送信し、専用のソフトウェアが動作する強力なサーバーで生成された結果を返す。 OSRM の一般公開されているサービスのような、オープンソースルーティングエンジンに基づくルーティングサービスは、R から呼び出された場合、ローカルにホストされたインスタンスと全く同じように動作し、単に更新される「ベース URL」を指定する引数を必要とするだけである。 しかし、外部のルーティングサービスは専用のマシンでホストされているため (通常、正確なルートを生成するインセンティブを持つ営利企業が資金を提供している)、以下のような利点がある。 @@ -519,13 +517,13 @@ R では [**cyclestreets**](https://rpackage.cyclestreets.net/) パッケージ ### ルーティング: 実例紹介 -前のセクションで生成された希望線をすべてルート探索\index{routing}する代わりに、政策上関心のある希望線\index{desire lines}に焦点を当てることにする。 +前のセクションで生成された希望線をすべてルート探索\index{るーとけんさく@ルート検索}する代わりに、政策上関心のある希望線\index{きぼうせん@希望線}に焦点を当てることにする。 データ全体を処理する前に、一部のデータに対して計算量の多い処理を実行することは、賢明な方法である。ルーティングも然り。 ルーティングは、ジオメトリが詳細でかつローとオブジェクトの属性が多いと、時間とメモリを消費し、オブジェクトが巨大になる。 このため、この節では、ルートを計算する前に希望線をフィルタする。 自転車旅行の利点は、自動車旅行を置き換えるときに最も大きく、比較的短い旅行は長い旅行よりも自転車に乗る可能性が高いという観察に基づいて、自転車旅行の可能性を推定することに重点を置いて希望線のサブセットをフィルタリングする [@lovelace_propensity_2017]。 -短い距離 (5 km 程度、時速 20 km で 15 分程度で走れる距離) は、比較的自転車で移動する確率が高く、電動自転車で移動すると最大距離が延びる [@lovelace_propensity_2017]。 +短いトリップ (5 km 程度、時速 20 km で 15 分程度で走れる距離) は、比較的自転車で移動する確率が高く、電動自転車で移動すると最大距離が延びる [@lovelace_propensity_2017]。 これらの考察に基づき、希望線をフィルタし、多くの (100 以上の) 短い (ユークリッド距離 2.5〜5 km) トリップを駆動する OD ペアを表すオブジェクト `desire_lines_short` を返す次のコードチャンクに反映しよう。 ```{r 13-transport-17, message=FALSE} @@ -535,7 +533,7 @@ desire_lines_short = desire_lines |> ``` 上記のコードで、`st_length()` は Section \@ref(distance-relations) にあるように、各希望線の長さを計算している。 -Section \@ref(vector-attribute-subsetting) で述べるように、**dplyr** の `filter()` 関数で、上記の条件に基づいて `desire_lines` データセットにフィルタをかけている\index{filter operation|see{attribute!subsetting}},。 +Section \@ref(vector-attribute-subsetting) で述べるように、**dplyr** の `filter()` 関数で、上記の条件に基づいて `desire_lines` データセットにフィルタをかけている\index{ふぃるたそうさ@フィルタ操作|{ぞくせい@属性!ぶぶんしゅうごう＠部分集合} 参照}。 次に、この希望線をルートに変換する。 これは一般に公開されている OSRM サービスを用いて、以下のコードにある **stplanr** 関数 `route()` と `route_osrm()`\index{stplanr (package)} で行われる。 @@ -544,12 +542,12 @@ routes_short = route(l = desire_lines_short, route_fun = route_osrm, osrm.profile = "bike") ``` -出力は `routes_short` で、(少なくとも OSRM ルーティングエンジンによれば) サイクリングに適したトランスポートネットワーク\index{network} 上のルートを表す `sf` オブジェクトで、各欲求行に対して一つずつ出力される。 +出力は `routes_short` で、(少なくとも OSRM ルーティングエンジンによれば) 自転車利用に適したトランスポートネットワーク\index{ねっとわーく@ネットワーク} 上のルートを表す `sf` オブジェクトで、各欲求行に対して一つずつ出力される。 注意：上記のコマンドのような外部のルーティングエンジンの呼び出しは、インターネット接続 (そして、今回は必要ないが、環境変数に保存された API キーも必要な場合がある) でのみ動作する。 `desire_lines` オブジェクトに含まれるカラムに加えて、新しいルートデータセットには `distance` (今回はルートの距離を参照) と `duration` (秒単位) のカラムが含まれ、それぞれのルートについて有用な追加情報を提供する可能性がある。 -サイクリングルートと一緒に、多くの短い車での移動が行われる希望線をプロットすることにする。 +自転車利用ルートと一緒に、多くの短い車での移動が行われる希望線をプロットすることにする。 ルートの幅を置き換えられる可能性のある車の旅の数に比例させることは、道路網への介入に優先順位をつける効果的な方法を提供する [@lovelace_propensity_2017]。 -以下のコードチャンクは、希望線とルートをプロットし、その結果、人々が短距離を運転するルートを示す Figure \@ref(fig:routes) になる。[^13-transport-8] +Figure \@ref(fig:routes) は、希望線とルートをプロットし、人々が短距離を運転するルートを示している (ソースコードは github.com/geocompx を参照)。[^13-transport-8] [^13-transport-8]: なお、赤いルートと黒い希望線は、まったく同じ地点からスタートしているわけではない。 これは、ゾーンの重心がルートネットワーク上に存在することはほとんどなく、重心に最も近い輸送ネットワークノードからルートが発信されるためである。 @@ -581,28 +579,29 @@ tm_shape(zones_od) + tm_scalebar() ``` -例えば `mapview::mapview(desire_carshort$geom_car)` でインタラクティブな地図にプロットしてみると、ブリストル中心部から約 10 km 北のブラッドリー・ストーク周辺で多くの短距離自動車旅行が行われていることがわかる。 -[Wikipedia](https://en.wikipedia.org/wiki/Bradley_Stoke) によると、ブラッドリー・ストークは「民間投資によって建設されたヨーロッパ最大のニュータウン」であり、公共交通機関の整備が限定的であることを示唆している。 -さらに、この町は、「M4高速道路とM5高速道路など、大規模な (サイクリングに不利な) 道路構造に囲まれている」 [@tallon_bristol_2007]。 +インタラクティブな地図で可視化してみると、Bristol 中心部から約 10 km 北の Bradley Stoke 周辺で多くの短距離自動車旅行が行われていることがわかる。 +[Wikipedia](https://en.wikipedia.org/wiki/Bradley_Stoke) によると、Bradley Stoke は「民間投資によって建設されたヨーロッパ最大のニュータウン」であり、公共交通機関の整備が限定的であることを示唆している。 +さらに、この町は、「M4 高速道路と M5 高速道路など、大規模な (自転車利用に不利な) 道路構造に囲まれている」 [@tallon_bristol_2007]。 -旅行希望線\index{desire lines}を旅行ルートに変換することは、政策の観点から多くの利点がある。 +旅行希望線\index{きぼうせん@希望線}を旅行ルートに変換することは、政策の観点から多くの利点がある。 ルーティングエンジンによって計算された正確なルートをたどるトリップがどれだけあるか (あったとしても) 確認することはできないことを覚えておくことが重要である。 しかし、ルートや道路・区間レベルの結果は、政策に大きく関連する可能性がある。 ルートセグメントの結果は、利用可能なデータ に従って、最も必要な場所に投資を優先させることを可能にすることができる [@lovelace_propensity_2017]。 ## ルートネットワーク {#route-networks} +\index{ねっとわーく@ネットワーク} 一般に路線は、希望線と同じレベルのデータ (または重複する可能性のあるセグメントのレベル) を含むが、路線網データセットは、交通網をほぼ完全に表現するものである。 路線網の各セグメント (交差点間の連続した道路区間にほぼ相当) は、一度だけ存在する。しかし、セグメントの平均長はデータソースによって異なる (このセクションで使用した OSM 由来の `bristol_ways` データセットのセグメントの平均長は 200 m 強で、標準偏差は 500 m 近い)。 セグメント長にばらつきがあるのは、地方では交差点が遠く、密集した都市部では数メートルおきに交差点があるなど、セグメントの切れ目があるためと思われる。 -路線網は、インプットである場合もあれば、アウトプットである場合もあり、両方である場合もある。 +路線網は、インプットの場合もあれば、アウトプットの場合もあり、両方の場合もある。 ルート計算を行う交通研究は、内部または外部のルーティングエンジンのルートネットワークデータを必要とする (後者の場合、ルート網データは必ずしも R にインポートされるわけではない)。 しかし、路線ネットワークは、多くの交通研究プロジェクトにおいて重要なアウトプットでもある。特定の区間で発生しうるトリップ数などのデータをまとめ、路線ネットワークとして表現することで、最も必要なところに優先的に投資することができる。 -\index{network} +\index{ねっとわーく@ネットワーク} -ルートレベルのデータから得られる出力としてルートネットワークを作成する方法を示すために、モーダルシフトの簡単なシナリオを想像してください。 -ルート距離 0～3 km の車移動の 50％ が自転車に置き換えられ、その割合はルート距離が 1 km 増えるごとに 10 ポイントずつ下がり、6 km の車移動の 20％ が自転車に置き換えられ、8 km 以上の車移動が自転車に置き換えられないと想像してみよう。 +ルートレベルのデータから得られる出力としてルートネットワークを作成する方法を示すために、モーダルシフトの簡単なシナリオを想像してみたい。 +ルート距離 0～3 km の車移動の 50% が自転車に置き換えられ、その割合はルート距離が 1 km 増えるごとに 10 ポイントずつ下がり、6 km の車移動の 20％ が自転車に置き換えられ、8 km 以上の車移動が自転車に置き換えられないと想像しよう。 これはもちろん非現実的なシナリオ [@lovelace_propensity_2017] ではあるものの、出発点としては有用であろう。 この場合、自動車から自転車へのモーダルシフトを次のようにモデル化することができる。 @@ -621,9 +620,9 @@ routes_short_scenario = routes_short |> sum(routes_short_scenario$bicycle) - sum(routes_short$bicycle) ``` -約 4000 のトリップが車から自転車に切り替わるというシナリオを作成したので、この更新されたモデル化されたサイクリング活動がどこで行われるかをモデル化することができるようになった。 +約 4000 のトリップが車から自転車に切り替わるというシナリオを作成したので、この更新されたモデル化された自転車利用活動がどこで行われるかをモデル化することができるようになった。 これには、**stplanr** パッケージの関数 `overline()` を使用する。 -この関数は、ルートと要約する属性の名前を含むオブジェクトを第1と第2の引数として受け取り、分岐点 (2 つ以上の線の形状が出会うところ) で線を切断し、それぞれのユニークなルートセグメント [@morgan_travel_2020] について集約した統計量を計算するものである。 +この関数は、ルートと要約する属性の名前を含むオブジェクトを第 1 と第 2 の引数として受け取り、分岐点 (2 つ以上の線の形状が出会うところ) で線を切断し、それぞれのユニークなルートセグメント [@morgan_travel_2020] について集約した統計量を計算するものである。 ```{r rnet1} route_network_scenario = overline(routes_short_scenario, attrib = "bicycle") @@ -655,8 +654,8 @@ OSM\index{OpenStreetMap}のダウンロードと準備の時間を節約する summary(bristol_ways) ``` -出力は、`bristol_ways` が輸送ネットワーク\index{network}上の6,000以上のセグメントを表していることを示している。 -このネットワークや他の地理的なネットワークは、数学的なグラフとして表現することができ、ネットワーク上のノード\index{node}とエッジ\index{edge}で接続されている。 +出力は、`bristol_ways` が輸送ネットワーク\index{ねっとわーく@ネットワーク}上の6,000以上のセグメントを表していることを示している。 +このネットワークや他の地理的なネットワークは、数学的なグラフとして表現することができ、ネットワーク上のノード\index{のーど@ノード}とエッジ\index{えっじ@エッジ}で接続されている。 このようなグラフを扱うために、多くのRパッケージが開発されており、特に **igraph**\index{igraph (package)} が有名である。 ルートネットワークを手動で `igraph` オブジェクトに変換することはできるが、地理的な属性は失われる。 この **igraph** の制限を克服するために、路線ネットワークをグラフと地理的な線として同時に表現する **sfnetworks**\index{sfnetworks (package)} パッケージ [@R-sfnetworks] が開発された。 @@ -682,9 +681,9 @@ ways_sfn 前のコードチャンクの出力(スペースの関係上、最終的な出力は最も重要な 8 行のみを含むように短縮されている)は、 `ways_sfn` がグラフ形式と空間形式の両方のノードとエッジを含む複合オブジェクトであることを表している。 `ways_sfn` は `sfnetwork` クラスで、 **igraph** パッケージの `igraph` クラスをベースにしている。 -下の例では、各エッジを通る最短ルート\index{shortest route}の数を意味する 'edge betweenness' \index{edge} が計算されている (詳しくは `?igraph::betweenness` を参照)。 -このデータセットに対して、エッジの間 隔を計算した結果を図に示す。図には、比較のために `overline()` 関数で計算した自転車ルートネットワークデータセットをオーバーレイで表示している。 -この結果は、グラフの各エッジがセグメントを表していることを示している。道路ネットワークの中心に近いセグメントは最も高い betweenness 値を持ち、一方、ブリストル中心部に近いセグメントはより高いサイクリングの可能性を持っていることが、これらの単純化されたデータセットに基づいて示されている。 +下の例では、各エッジを通る最短ルート\index{さいたんるーと@最短ルート}の数を意味する 'edge betweenness' \index{えっじ@エッジ} が計算されている (詳しくは `?igraph::betweenness` を参照)。 +このデータセットに対して、エッジの間隔を計算した結果を図に示す。図には、比較のために `overline()` 関数で計算した自転車ルートネットワークデータセットをオーバーレイで表示している。 +この結果は、グラフの各エッジがセグメントを表していることを示している。道路ネットワークの中心に近いセグメントは最も高い betweenness 値を持ち、一方、Bristol 中心部に近いセグメントはより高い自転車利用の可能性を持っていることが、これらの単純化されたデータセットに基づいて示されている。 ```{r wayssln-gen} ways_centrality = ways_sfn |> @@ -692,7 +691,7 @@ ways_centrality = ways_sfn |> mutate(betweenness = tidygraph::centrality_edge_betweenness(lengths)) ``` -```{r wayssln, fig.cap="路線ネットワークデータセットの図。グレーの線は簡略化された道路網を表し、セグメントの太さは betweenness に比例する。緑色の線は、上記のコードで計算された潜在的なサイクリングフロー (片道) である。", fig.scap="Illustration of a small route network.", dev="ragg_png"} +```{r wayssln, fig.cap="路線ネットワークデータセットの図。グレーの線は簡略化された道路網を表し、セグメントの太さは betweenness に比例する。緑色の線は、上記のコードで計算された潜在的な自転車利用フロー (片道) である。", fig.scap="Illustration of a small route network.", dev="ragg_png"} bb_wayssln = tmaptools::bb(route_network_scenario, xlim = c(0.1, 0.9), ylim = c(0.1, 0.6), relative = TRUE) tm_shape(zones_od) + tm_fill(fill_alpha = 0.2, lwd = 0.1) + @@ -716,7 +715,7 @@ tm_shape(zones_od) + # gdf = igraph::as_long_data_frame(ways_sfn@g) ``` -また、**sfnetworks** パッケージを用いると、このルートネットワークのグラフ表現を使って、出発地と目的地の間の最短ルート\index{shortest route}を求めることができる。 +また、**sfnetworks** パッケージを用いると、このルートネットワークのグラフ表現を使って、出発地と目的地の間の最短ルート\index{さいたんるーと@最短ルート}を求めることができる。 このセクションで紹介した方法は比較的単純で、実際にはもっと可能なことはある。 **sfnetworks** が提供するグラフと空間の二つの機能により、多くの新しい強力な技法が可能になるが、このセクションで完全にカバーすることはできない。 @@ -731,8 +730,8 @@ tm_shape(zones_od) + 本章で紹介するデータ駆動型アプローチの利点は、モジュール化されていることである。 この段階に至るまでには、(希望線から生成された)短いが車に依存する通勤経路を特定するための手順があり、(ルート・ネットワーク)セクションで **sfnetworks** パッケージを使用してルート・ネットワークの特性を分析することが含まれている。 -本章の最後のコードチャンクは、サイクリングインフラから短い距離のエリアを表す新しいデータセットに、前のセクションのサイクリング潜在力の推定値を重ねることによって、これらの一連の分析を結合する。 -この新しいデータセットは、以下のコードで作成される。1) 交通ネットワークを表す`bristol_ways`オブジェクトからサイクルウェイのエンティティをフィルタリングする。2) サイクルウェイの個々の LINESTRING エンティティを単一の multilinestring オブジェクトに「統合」する (バッファ作成が速くなるため)。 +本章の最後のコードチャンクは、自転車利用インフラから短い距離のエリアを表す新しいデータセットに、前のセクションの自転車利用潜在力の推定値を重ねることによって、これらの一連の分析を結合する。 +この新しいデータセットは、以下のコードで作成される。1) 交通ネットワークを表す `bristol_ways` オブジェクトから自転車道のエンティティをフィルタリングする。2) 自転車道の個々の線エンティティを単一の複合線オブジェクトに「統合」する (バッファ作成が速くなるため)。 ```{r 13-transport-25} existing_cycleways_buffer = bristol_ways |> @@ -741,7 +740,7 @@ existing_cycleways_buffer = bristol_ways |> st_buffer(dist = 100) # 3) create buffer ``` -次の段階は、ネットワークの中で、サイクリングの可能性が高いが、サイクリングのための設備がほとんどないポイントを表すデータセットを作成することである。 +次の段階は、ネットワークの中で、自転車利用の可能性が高いが、自転車利用のための設備がほとんどない点を表すデータセットを作成することである。 ```{r, echo=FALSE, eval=FALSE} waldo::compare( @@ -779,7 +778,7 @@ qtm(route_network_no_infra, basemaps = leaflet::providers$Esri.WorldTopoMap, -```{r cycleways, echo=FALSE, message=FALSE, fig.cap="ブリストルにおける自動車依存度を下げるために、自転車インフラを優先的に整備するルートの候補。静的マップは、既存のインフラと自動車と自転車の乗り換えの可能性が高いルートとの間のオーバーレイの概要を提供する (左)。`qtm()` 関数から生成されたインタラクティブマップのスクリーンショットは、新しいサイクルウェイから利益を得ることができる場所としてWhiteladies Roadを強調している (右)。", out.width="50%", fig.show='hold', fig.scap="Routes along which to prioritise cycle infrastructure.", fig.height=9, dev="ragg_png"} +```{r cycleways, echo=FALSE, message=FALSE, fig.cap="ブリストルにおける自動車依存度を下げるために、自転車インフラを優先的に整備するルートの候補。静的マップは、既存のインフラと自動車と自転車の乗り換えの可能性が高いルートとの間のオーバーレイの概要を提供する (左)。`qtm()` 関数から生成されたインタラクティブマップのスクリーンショットは、新しい自転車道から利益を得ることができる場所として Whiteladies Road を強調している (右)。", out.width="50%", fig.show='hold', fig.scap="Routes along which to prioritise cycle infrastructure.", fig.height=9, dev="ragg_png"} # Previous verson: # source("https://github.com/Robinlovelace/geocompr/raw/main/code/13-cycleways.R") tm_shape(existing_cycleways_buffer, bbox = bristol_region) + @@ -792,16 +791,16 @@ tm_shape(existing_cycleways_buffer, bbox = bristol_region) + knitr::include_graphics("figures/bristol_cycleways_zoomed.png") ``` -この方法には限界がある。現実には、人々はゾーン重心に移動したり、特定のモードの最短ルート\index{shortest route}アルゴリズムを常に使用するわけではない。 +この方法には限界がある。現実には、人々はゾーン重心に移動したり、特定のモードの最短ルート\index{さいたんるーと@最短ルート}アルゴリズムを常に使用するわけではない。 しかし、この結果は、自動車依存と公共交通の観点から、自転車専用道路を優先的に整備できるルートを示している。 この分析は、現実の交通計画の立案に使うためには、より大きなデータを使うなど、大幅に拡大して行う必要がある。 ## 今後の方向性 {#future-directions-of-travel} -この章では、交通研究にジオコンピューティングを利用する可能性を紹介し、オープンデータと再現可能なコードを用いて、都市の交通システムを構成するいくつかの重要な地理的要素を調査した。 +この章では、交通研究にジオコンピュテーションを利用する可能性を紹介し、オープンデータと再現可能なコードを用いて、都市の交通システムを構成するいくつかの重要な地理的要素を調査した。 この結果は、どこに投資が必要かを計画するのに役立つだろう。 -交通システムは複数の相互作用レベルで機能しているため、ジオコンピューティングの手法は交通システムの仕組みや異なる介入の効果を理解する上で大きな可能性を秘めている。 +交通システムは複数の相互作用レベルで機能しているため、ジオコンピュテーションの手法は交通システムの仕組みや異なる介入の効果を理解する上で大きな可能性を秘めている。 この分野でできることはまだまだたくさんある。この章で紹介した基礎の上に、さまざまな方向性を構築することが可能であろう。 交通は、多くの国で温室効果ガスの排出源として最も急速に増加しており、「特に先進国では最大の温室効果ガス排出部門」になると言われている ([EURACTIV.com](https://www.euractiv.com/section/agriculture-food/opinion/transport-needs-to-do-a-lot-more-to-fight-climate-change/) を参照)。 交通機関の排出量は社会全体で非常に不平等に配分されており、交通機関は (食料や暖房とは異なり) 幸福に不可欠なものではない。 @@ -812,16 +811,16 @@ knitr::include_graphics("figures/bristol_cycleways_zoomed.png") こうした「交通の未来」を地域や国レベルでさらに探求することで、新たな知見を得られるだろう。 方法論的には、本章で示した基礎は、より多くの変数を分析に含めることで拡張することが可能である。 -速度制限、交通量の多さ、自転車や歩行者用保護道の設置などのルートの特徴は、「モーダルスプリット」 (modal split, 異なる交通手段による旅行の割合) に関連付けることができる。 +速度制限、交通量の多さ、自転車や歩行者用保護道の設置などのルートの特徴は、「モーダルスプリット」 (modal split、異なる交通手段による旅行の割合) に関連付けることができる。 例えば、OpenStreetMap\index{OpenStreetMap} のデータをバッファや、Chapter \@ref(attr) と Chapter \@ref(spatial-operations) で紹介した地理データ手法で集約すれば、交通ルートの近くに緑地があるかどうかを検出することも可能である。 -R\index{R} の統計モデリング機能を使えば、例えば、現在と将来のサイクリングのレベルを予測することができるだろう。 +R\index{R} の統計モデリング機能を使えば、例えば、現在と将来の自転車利用のレベルを予測することができるだろう。 -この種の分析のベースには、Propensity to Cycle Tool (PCT)がある。 これは、R\index{R} で開発された一般にアクセス可能な ( [www.pct.bike](https://www.pct.bike/) 参照) マッピングツールで、イングランド全域のサイクリングへの投資を優先させるために使用されている [@lovelace_propensity_2017]。 +この種の分析のベースには、Propensity to Cycle Tool (PCT)がある。 これは、R\index{R} で開発された一般にアクセス可能な ( [www.pct.bike](https://www.pct.bike/) 参照) マッピングツールで、イングランド全域の自転車利用への投資を優先させるために使用されている [@lovelace_propensity_2017]。 同様のツールは、世界中の大気汚染や公共交通機関へのアクセスなど、他のテーマに関連したエビデンスに基づく交通政策を奨励するためにも使用できる。 ## 演習 {#ex-transport} ```{r, echo=FALSE, results='asis'} -res = knitr::knit_child('_13-ex.Rmd', quiet = TRUE, options = list(include = FALSE, eval = FALSE)) +res = knitr::knit_child('_13-ex-ja.Rmd', quiet = TRUE, options = list(include = FALSE, eval = FALSE)) cat(res, sep = '\n') ``` \ No newline at end of file diff --git a/14-location-ja.Rmd b/14-location-ja.Rmd index 96ce263..07096cc 100644 --- a/14-location-ja.Rmd +++ b/14-location-ja.Rmd @@ -212,12 +212,11 @@ reclass # 出力は一部省略 #> max values : 8000, 3, 3, 3 ``` - ## 大都市圏を定義 {#define-metropolitan-areas} -大都市圏とは、50万人以上が住む 20 km^2^ のピクセルと定義している。 +大都市圏とは、50 万人以上が住む 20 km^2^ のピクセルと定義している。 この粗い解像度のピクセルは、Section \@ref(aggregation-and-disaggregation) で紹介したように、`aggregate()`\index{aggregation}、速やかに作成することができる。 -以下のコマンドは、引数 `fact = 20`、結果の解像度を20倍にしている (元のラスタの解像度が1km^2^であったことを思い出してみよう)。 +以下のコマンドは、引数 `fact = 20`、結果の解像度を 20 倍にしている (元のラスタの解像度が 1 km^2^ であったことを思い出してみよう)。 ```{r 14-location-11, warning=FALSE, cache=TRUE, cache.lazy=FALSE} pop_agg = aggregate(reclass$pop, fact = 20, fun = sum, na.rm = TRUE) @@ -233,7 +232,7 @@ pop_agg = pop_agg[pop_agg > 500000, drop = FALSE] これをプロットすると、8 つの大都市圏 (Figure \@ref(fig:metro-areas)) が見えてくる。 各領域は、1つ以上のラスタセルで構成される。 1つのリージョンに属するすべてのセルを結合できればコマンドは、 -**terra**\index{terra} の `patches()` である。 +**terra**\index{terra (package)} の `patches()` である。 その後、`as.polygons()` でラスタオブジェクトを空間ポリゴンに変換し、`st_as_sf()` で `sf`. オブジェクトに変換する。 ```{r 14-location-13, warning=FALSE, message=FALSE} @@ -294,7 +293,7 @@ metro_names = metro_names$city |> - `map()`\index{loop!map} (`lapply()`\index{loop!lapply} の **tidyverse** 相当)。これは、OSM \index{OpenStreetMap} クエリー関数 `opq()` (Section \@ref(retrieving-data) 参照) のバウンディングボックス\index{bounding box} を定義する、8 つの大都市名すべてを繰り返し処理するものである - `add_osm_feature()` で、キー値が `shop` の OSM\index{OpenStreetMap} 要素を指定する (共通のキー:値のペアの一覧は [wiki.openstreetmap.org](https://wiki.openstreetmap.org/wiki/Map_Features) を参照) - `osmdata_sf()`、これは OSM\index{OpenStreetMap} データを空間オブジェクト (クラス `sf`) に変換するものである -- `while()`\index{loop!while}、1 回目のダウンロードに失敗すると繰り返し (今回は 3 回) ダウンロードを試みる。^[OSM-downloadは1回目で失敗することもあるようである +- `while()`\index{loop!while}、ダウンロードに失敗すると、さらに 2 回ダウンロードを試みる。^[OSM-download は 1 回目で失敗することもあるようである ] このコードを実行する前に: 約 2 GB のデータをダウンロードすることを考慮してみよう。 @@ -457,6 +456,6 @@ if (knitr::is_latex_output()) { ## 演習 ```{r, echo=FALSE, results='asis'} -res = knitr::knit_child('_14-ex.Rmd', quiet = TRUE, options = list(include = FALSE, eval = FALSE)) +res = knitr::knit_child('_14-ex-ja.Rmd', quiet = TRUE, options = list(include = FALSE, eval = FALSE)) cat(res, sep = '\n') ``` diff --git a/15-eco-ja.Rmd b/15-eco-ja.Rmd index b62a543..fbec855 100644 --- a/15-eco-ja.Rmd +++ b/15-eco-ja.Rmd @@ -40,27 +40,25 @@ library(vegan) # 生態学 残念ながら、霧のオアシスは、主に人間の活動により、大きく危機に瀕している。 残されたユニークな植生生態系を効果的に保護するためには、原生植物相の構成と空間分布に関するエビデンスが必要である [@muenchow_predictive_2013; @muenchow_soil_2013]。 -この章では、前の章で学んだ技術を生態学的に実際に応用していこう。 -その中で、ペルー中央北岸のカスマ近郊にある *lomas* のある山、Mongón 山の南斜面における維管束植物の構成と空間分布を分析する (Figure \@ref(fig:study-area-mongon))。 + +この章では、Peru の中央北岸に位置するカスマ近郊の**ロマ植生**山である Mongón 山の南斜面における維管束植物 (ここでは主に顕花植物) の組成と空間分布を分析する (Figure \@ref(fig:study-area-mongon))。 +モゴン山での野外調査において、2011年の冬に無作為にサンプリングした 4 x 4 m^2^ の 100 区画に生息する維管束植物をすべて記録した [@muenchow_predictive_2013]。 +このサンプリングは、米国海洋大気庁 ([NOAA](https://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php)) が発表したデータに示されているように、その年の強いラニーニャ現象が発生した時期と重なった。 +このため、沿岸部の砂漠では例年よりもさらに高いレベルの乾燥が発生し、Peru の**ロマ植生**山脈の南斜面では霧の活動が活発化した。 ```{r study-area-mongon, echo=FALSE, fig.cap="Mongón 山調査地、Muenchow, Schratz, and Brenning (2017) より。", out.width="60%", fig.scap="The Mt. Mongón study area."} knitr::include_graphics("figures/15_study_area_mongon.png") # knitr::include_graphics("https://user-images.githubusercontent.com/1825120/38989956-6eae7c9a-43d0-11e8-8f25-3dd3594f7e74.png") ``` -Mongón 山への野外調査において、2011年の冬、ランダムにサンプルした 4x4 m^2^ の 100 区画に生息するすべての維管束植物を記録した [@muenchow_predictive_2013]。 -National Oceanic and Atmospheric Administration ([NOASS](https://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php)を参照) の報告書にある通り、サンプル採取は、その年の強いラニーニャ現象に重なった。 -そのため、沿岸部の砂漠は通常よりもさらに高い乾燥度を示すようになり、ペルーの *lomas* のある山脈の南斜面では、霧の活動も活発になっている。 - - +また本章では、前章で扱ったテクニックを、重要な応用分野である生態学に応用する方法を示す。 +具体的には、 -Ordinations\index{ordination} は、 (ノイズの多い) データセットから主要な勾配を抽出するための次元削減技術であり、私たちの場合は南山斜面に沿って発達した植物学的勾配である (次のセクションを参照)。 -本章では、最初の序列軸である植物相の勾配を、標高、傾斜、集水域\index{catchment area}および NDVI\index{NDVI} といった環境予測因子の関数としてモデル化することにする。 -このために、非常に有名な機械学習\index{machine learning}アルゴリズムであるランダムフォレストモデル\index{random forest}を利用する [@breiman_random_2001]。 -このモデルを使えば、調査地のどこでも植物組成の空間分布図作成が可能になる。 -最適な予測を保証するために、空間クロスバリデーション \index{cross-validation!spatial CV} (Section \@ref(svm) 参照) の助けを借りて、ハイパーパラメータ\index{hyperparameter}を事前に調整することが推奨される。 +- 必要なデータを読み込み、環境予測因子を計算する (Section \@ref(data-and-data-preparation)) +- 次元削減技術 (座標付け; Section \@ref(nmds)) を用いて、種組成行列から主な植物学的勾配を抽出する\index{じげんさくげんぎじゅつ@次元削減技術} +- 最初の順序軸、すなわち植物相勾配を、標高、傾斜、集水域、NDVI などの環境予測変数の関数としてモデル化する (Section \@ref(modeling-the-floristic-gradient))\index{ざひょうづけ@座標付け}\index{しゅうすいいき@集水域}\index{NDVI} +そのために、ランダムフォレストモデルを利用する。機械学習アルゴリズム [@breiman_random_2001]の1つである。最適な予測を行うために、空間交差検証を用いて、ハイパーパラメータを事前に調整することが望ましい (see Section \@ref(svm))\index{らんだむふぉれすと@ランダムフォレスト}\index{きかいがくしゅう@機械学習}\index{はいぱーぱらめーた@ハイパーパラメータ}\index{くろすばりでーしょん@クロスバリデーション!くうかんてき@空間的 CV} +- 調査地内の任意の場所の植物組成の空間分布図を作成する (Section \@ref(predictive-mapping)) ## データとデータ準備 {#data-and-data-preparation} @@ -219,12 +217,12 @@ random_points = cbind(random_points, ep_rp) ## 次元性を低減 {#nmds} -順序付け (Ordination\index{ordination}) は、植生学において、0 で埋め尽くされた大規模な種間プロット行列から主要情報 (生態的勾配に相当することが多い) を抽出するための一般的なツールである。 +座標付け (Ordination\index{ざひょうづけ@座標付け}) は、植生学において、0 で埋め尽くされた大規模な種間プロット行列から主要情報 (生態的勾配に相当することが多い) を抽出するための一般的なツールである。 しかし、リモートセンシング\index{remote sensing}、土壌学、ジオマーケティング\index{geomarketing} などの分野でも利用されている。 -順序付け\index{ordination}テクニックに馴染みがない場合、または復習が必要な場合は、生態学で人気の順序付けテクニックを簡単に紹介した Michael W. Palmer の [web page](https://ordination.okstate.edu/overview.htm) と R でこれらのテクニックを適用する方法について深く調べた @borcard_numerical_2011 に目を通してみよう。 +座標付け\index{ざひょうづけ@座標付け}テクニックに馴染みがない場合、または復習が必要な場合は、生態学で人気の座標付けテクニックを簡単に紹介した Michael W. Palmer の [web page](https://ordination.okstate.edu/overview.htm) と R でこれらのテクニックを適用する方法について深く調べた @borcard_numerical_2011 に目を通してみよう。 **vegan**\index{vegan (package)} のパッケージのドキュメントも非常に有用なリソースである (`vignette(package = "vegan")`)。 -主成分分析 (principal component analysis, PCA\index{PCA}) は、おそらく最も有名な順序付け \index{ordination}の手法である。 +主成分分析 (principal component analysis, PCA\index{しゅせいぶんぶんせき@主成分分析}) は、おそらく最も有名な座標付け \index{ざひょうづけ@座標付け}の手法である。 変数間の線形関係が期待でき、2 つのプロット (オブザベーション) における変数の共同不在が類似性とみなせる場合、次元を削減するためのすばらしいツールである。 これは植生データではほとんどない。 @@ -235,15 +233,15 @@ random_points = cbind(random_points, ep_rp) というのも、この2つの全く異なる環境設定に共通するのは、 (稀少なユビキタス種を除いて) 種が存在しないという点だけである可能性が非常に高いからである。 Non-metric multidimensional scaling (NMDS\index{NMDS}) は、生態学でよく使われる次元削減手法の一つである [@vonwehrden_pluralism_2009]。 -NMDS\index{NMDS} は、元の行列のオブジェクト間の距離と、順序付けられたオブジェクト間の距離の間のランクベースの差異を低減する。 +NMDS\index{NMDS} は、元の行列のオブジェクト間の距離と、座標付けられたオブジェクト間の距離の間のランクベースの差異を低減する。 その差をストレスとして表現している。 -ストレス値が低いほど、順序付け、すなわち元の行列の低次元表現が良好であることを示す。 +ストレス値が低いほど、座標付け、すなわち元の行列の低次元表現が良好であることを示す。 ストレス値が 10 より小さいと適合性が高く、15程度でも良好、20より大きいと適合性が低いことを表している [@mccune_analysis_2002]。 R では、**vegan**\index{vegan (package)} パッケージの `metaMDS()` で NMDS を実行することができる。 入力として、サイトを行、種を列とする群集行列を期待する。 -しばしば、有-無データを用いた順序付け\index{ordination} は、 (説明される分散の点で) より良い結果をもたらするが、その代償として、もちろん、入力行列の情報量は少なくなる (演習も参照)。 -`decostand()` は、数値の観測結果を、1が種の発生、0が種の不在を示す有無に変換する。 -NMDS\index{NMDS} のようなオーダリング技術では、部位ごとに少なくとも1回の観測が必要である。 +しばしば、有-無データを用いた座標付け\index{ざひょうづけ@座標付け} は、 (説明される分散の点で) より良い結果をもたらするが、その代償として、もちろん、入力行列の情報量は少なくなる (演習も参照)。 +`decostand()` は、数値の観測結果を、1 が種の発生、0 が種の不在を示す有無に変換する。 +NMDS\index{NMDS} のようなオーダリング技術では、部位ごとに少なくとも 1 回のオブザベーションが必要である。 したがって、種が発見されなかったサイトはすべて除外する必要がある。 ```{r 15-eco-11} @@ -251,13 +249,13 @@ NMDS\index{NMDS} のようなオーダリング技術では、部位ごとに少 pa = vegan::decostand(comm, "pa") # 100 行 (箇所) , 69 列 (種) # 1種以上発見されたサイトのみを残す pa = pa[rowSums(pa) != 0, ] # 84 行, 69 列 (種) -``` +``` 得られた行列は、NMDS\index{NMDS} の入力として機能する。 `k` は出力軸数を指定し、ここでは 4 とする。 ^[ `k` を選択する一つの方法として、1 から 6 の間で `k` を試し、最も良いストレス値が得られる結果を使用する方法がある [@mccune_analysis_2002]。 ] -NMDS\index{NMDS} は、各ステップで順序付けされた空間をより入力行列に近づけようとする反復的な手順である。 +NMDS\index{NMDS} は、各ステップで座標付けされた空間をより入力行列に近づけようとする反復的な手順である。 アルゴリズムの収束を確認するために、`try` パラメータを使ってステップ数を 500 に設定した。 ```{r 15-eco-12, eval=FALSE, message=FALSE} @@ -283,12 +281,12 @@ saveRDS(nmds, "extdata/15-nmds.rds") nmds = readRDS("extdata/15-nmds.rds") ``` -ストレス値 9 は非常に良い結果を表し、縮小された順序付け空間が入力行列の分散の大部分を表していることを意味する。 -全体として、NMDS\index{NMDS} は、順序付け空間において、 (種の構成という点で) より類似したオブジェクトがより近くに配置される。 -しかし、他の多くの順序付け \index{ordination} の手法とは対照的に、軸は任意であり、必ずしも重要度 [@borcard_numerical_2011] によって順序付けされるわけではない。 +ストレス値 9 は非常に良い結果を表し、縮小された座標付け空間が入力行列の分散の大部分を表していることを意味する。 +全体として、NMDS\index{NMDS} は、座標付け空間において、(種の構成という点で) より類似したオブジェクトがより近くに配置される。 +しかし、他の多くの座標付け \index{ざひょうづけ@座標付け} の手法とは対照的に、軸は任意であり、必ずしも重要度 [@borcard_numerical_2011] によって座標付けされるわけではない。 しかし、調査地では湿度が主な勾配を表していることが既に分かっている [@muenchow_predictive_2013;@muenchow_rqgis:_2017]。 湿度は標高と高い相関があるため、標高に応じて NMDS\index{NMDS} 軸を回転させる (NMDS 軸の回転の詳細については `?MDSrotate` も参照)。 -結果をプロットしてみると、意図したとおり、第1軸は明らかに高度 (Figure \@ref(fig:xy-nmds)) と関連していることがわかる。 +結果をプロットしてみると、意図したとおり、第 1 軸は明らかに高度 (Figure \@ref(fig:xy-nmds)) と関連していることがわかる。 ```{r xy-nmds-code, fig.cap="NMDS の第1軸を高度に対してプロット。", fig.scap = "First NMDS axis against altitude plot.", fig.asp=1, out.width="60%", eval=FALSE} elev = dplyr::filter(random_points, id %in% rownames(pa)) |> @@ -363,15 +361,15 @@ plot(rotnmds, display = "sites") plot(fit_2) ``` -最初の NMDS\index{NMDS} 軸のスコアは、Mont. Mongón の斜面に沿って現れる異なる植生形態、すなわち植物学的勾配を表している。 -それらを空間的に可視化するために、NMDS\index{NMDS} のスコアを先に作成した予測因子 (Section \@ref(data-and-data-preparation)) でモデル化し、得られたモデルを予測マッピングに利用する (次節参照) ことができる。 +最初の NMDS\index{NMDS} 軸のスコアは、Mongón 山の斜面に沿って現れる異なる植生形態、すなわち植物学的勾配を表している。 +それらを空間的に可視化するために、NMDS\index{NMDS} のスコアを先に作成した予測因子 (Section \@ref(data-and-data-preparation)) でモデル化し、得られたモデルを予測地図に利用する (次節参照) ことができる。 ## 植物相の勾配をモデル化 {#modeling-the-floristic-gradient} -植物相の勾配を空間的に予測するために、ランダムフォレスト\index{random forest}モデルを使用する。 -ランダムフォレスト\index{random forest}モデルは、環境・生態系のモデリングに頻繁に適用され、予測性能の面で最良の結果をもたらすことが多い [@hengl_random_2018;@schratz_hyperparameter_2019]。 -ここで、これらはランダムフォレスト\index{random forest}の基礎を形成するものであるため、決定木とバギングについて簡単に紹介する。 -ランダムフォレスト\index{random forest}と関連する技術についてのより詳細な説明は @james_introduction_2013 を参照。 +植物相の勾配を空間的に予測するために、ランダムフォレスト\index{らんだむふぉれすと@ランダムフォレスト}モデルを使用する。 +ランダムフォレスト\index{らんだむふぉれすと@ランダムフォレスト}モデルは、環境・生態系のモデリングに頻繁に適用され、予測性能の面で最良の結果をもたらすことが多い [@hengl_random_2018;@schratz_hyperparameter_2019]。 +ここで、これらはランダムフォレスト\index{らんだむふぉれすと@ランダムフォレスト}の基礎を形成するものであるため、決定木とバギングについて簡単に紹介する。 +ランダムフォレスト\index{らんだむふぉれすと@ランダムフォレスト}と関連する技術についてのより詳細な説明は @james_introduction_2013 を参照。 決定木を例として紹介すると、まず、回転した NMDS\index{NMDS} スコアと現場観測値 (`random_points`) を結合して、応答予測行列を構築する。 また、得られたデータフレームは、後で **mlr3**\index{mlr3 (package)} のモデリングに使用する予定である。 @@ -405,10 +403,7 @@ knitr::include_graphics("figures/15_tree.png") 結果として得られる木は、3 つの内部ノードと 4 つの終端ノード (Figure \@ref(fig:tree)) で構成される。 木の一番上にある最初の内部ノードは、328.5 m 以下のすべてのオブザベーションを左に、それ以外のすべてのオブザベーションを右の枝に割り当てる。 - - 左の枝に入るオブザベーションは、平均 NMDS\index{NMDS} スコアが -1.198 である。 - 全体として、標高が高いほどNMDS\index{NMDS} のスコアが高くなる、というように解釈できる。 つまり、単純な決定木によって、すでに4つの異なる植物群像が明らかにされているのである。 詳しく学びたい方は Section \@ref(predictive-mapping) を参照。 @@ -417,8 +412,8 @@ knitr::include_graphics("figures/15_tree.png") アンサンブル技術は、複数のモデルの予測値を単純に結合するものである。 このように、バギングでは同じ入力データから繰り返しサンプルを取り、その予測値を平均化する。 これにより、分散と過適合\index{overfitting}を減らし、決定木と比較してはるかに優れた予測精度を実現する。 -最後に、ランダムフォレスト\index{random forest}は、相関の高い木の予測を平均化すると、相関の低い木の予測を平均化するよりも分散が大きく、信頼性が低くなるので、バギングを拡張して改良することが望ましい [@james_introduction_2013]。 -これを実現するために、ランダムフォレスト\index{random forest}はバギングを使用するが、従来のバギングでは各木が利用可能なすべての予測子を使用できるのとは対照的に、ランダムフォレストは利用可能なすべての予測子のランダムサンプルだけを使用する。 +最後に、ランダムフォレスト\index{らんだむふぉれすと@ランダムフォレスト}は、相関の高い木の予測を平均化すると、相関の低い木の予測を平均化するよりも分散が大きく、信頼性が低くなるので、バギングを拡張して改良することが望ましい [@james_introduction_2013]。 +これを実現するために、ランダムフォレスト\index{らんだむふぉれすと@ランダムフォレスト}はバギングを使用するが、従来のバギングでは各木が利用可能なすべての予測子を使用できるのとは対照的に、ランダムフォレストは利用可能なすべての予測子のランダムサンプルだけを使用する。 既に入力変数 (`rp`) を構築しているので、**mlr3**\index{mlr3 (package)} の構成要素 (タスク、学習器、リサンプリング) を指定するための準備は全て整っている。 @@ -460,7 +455,7 @@ task = mlr3spatiotempcv::as_task_regr_st(select(rp, -id, -spri), バックエンドとして `sf` オブジェクトを使用すると、後の空間分割に必要なジオメトリ情報が自動的に提供される。 さらに、`id` と `spri` の列は、これらの変数をモデリングにおける予測因子として使用すべきではないため、削除した。 -次に、**ranger** パッケージのランダムフォレスト\index{random forest} 学習器を構築する [@wright_ranger_2017]。 +次に、**ranger** パッケージのランダムフォレスト\index{らんだむふぉれすと@ランダムフォレスト} 学習器を構築する [@wright_ranger_2017]。 ```{r 15-eco-21} lrn_rf = lrn("regr.ranger", predict_type = "response") @@ -468,7 +463,7 @@ lrn_rf = lrn("regr.ranger", predict_type = "response") 例えばサポートベクタマシン\index{SVM} (Section \@ref(svm) 参照) とは対照的に、ランダムフォレストはハイパーパラメータのデフォルト値で使用した場合、既に良い性能を示すことが多い (これが人気の理由の一つだろう)。 それでも、チューニングによってモデル結果が適度に改善されることが多いので、努力する価値はある [@probst_hyperparameters_2018]。 -ランダムフォレスト\index{random forest}では、ハイパーパラメータ\index{hyperparameter} `mtry`、`min.node.size`、`sample.fraction` がランダム性の度合いを決定するので、これらを調整する必要がある [@probst_hyperparameters_2018]。 +ランダムフォレスト\index{らんだむふぉれすと@ランダムフォレスト}では、ハイパーパラメータ\index{はいぱーぱらめーた@ハイパーパラメータ} `mtry`、`min.node.size`、`sample.fraction` がランダム性の度合いを決定するので、これらを調整する必要がある [@probst_hyperparameters_2018]。 `mtry` は、各ツリーでいくつの予測変数を使用すべきかを示す。 すべての予測変数が使用される場合、これは事実上バギングに相当する (Section \@ref(modeling-the-floristic-gradient) の冒頭を参照)。 `sample.fraction` パラメータは、各ツリーで使用されるオブザベーションの割合を指定する。 @@ -476,7 +471,7 @@ lrn_rf = lrn("regr.ranger", predict_type = "response") `min.node.size` パラメータは、端末ノードが少なくとも持つべき観測値の数を示す (Figure \@ref(fig:tree) も参照)。 当然ながら、木や演算時間が大きくなればなるほど、`min.node.size` は小さくなる。 -ハイパーパラメータ\index{hyperparameter}の組み合わせはランダムに選択されるが、特定のチューニング限界 (`paradox::ps()` で作成) の範囲内に収まる必要がある。 +ハイパーパラメータ\index{はいぱーぱらめーた@ハイパーパラメータ}の組み合わせはランダムに選択されるが、特定のチューニング限界 (`paradox::ps()` で作成) の範囲内に収まる必要がある。 `mtry` は 1 から予測変数の数 `r ncol(task$data()) - 1`) までの範囲でなければならない。`sample.fraction` は 0.2 から 0.9 の間、`min.node.size` は 1 から 10 の間でなければならない [@probst_hyperparameters_2018]。 ```{r 15-eco-22, eval=FALSE} @@ -489,13 +484,13 @@ search_space = paradox::ps( ``` 探索空間を定義したことで、`AutoTuner()` 関数でチューニングを指定する準備が整いた。 -地理的なデータを扱うので、今回も空間クロスバリデーションを用いてハイパーパラメータ\index{hyperparameter}を調整する (Section \@ref(intro-cv) と Section \@ref(spatial-cv-with-mlr3) を参照)。 +地理的なデータを扱うので、今回も空間クロスバリデーションを用いてハイパーパラメータ\index{はいぱーぱらめーた@ハイパーパラメータ}を調整する (Section \@ref(intro-cv) と Section \@ref(spatial-cv-with-mlr3) を参照)。 具体的には、1回だけ繰り返す5分割の空間分割 (`rsmp()`) を使用することにする。 -これらの空間分割のそれぞれにおいて、あらかじめ定義された限界値 (`seach_space`) の範囲内でランダムに選択されたハイパーパラメータ構成 (`tnr()` ) を用いながら50個のモデル(`trm()`)を実行し、最適なハイパーパラメータ\index{hyperparameter}の組合せを見出す [Section \@ref(svm) と https://mlr3book.mlr-org.com/chapters/chapter4/hyperparameter_optimization.html#sec-autotuner, @becker_mlr3_2022 を参照]。 +これらの空間分割のそれぞれにおいて、あらかじめ定義された限界値 (`seach_space`) の範囲内でランダムに選択されたハイパーパラメータ構成 (`tnr()` ) を用いながら50個のモデル(`trm()`)を実行し、最適なハイパーパラメータ\index{はいぱーぱらめーた@ハイパーパラメータ}の組合せを見出す [Section \@ref(svm) と https://mlr3book.mlr-org.com/chapters/chapter4/hyperparameter_optimization.html#sec-autotuner, @becker_mlr3_2022 を参照]。 性能指標は二乗平均平方根誤差 (root mean squared error, RMSE\index{RMSE}) である。 ```{r 15-eco-23, eval=FALSE} -autotuner_rf = mlr3tuning::AutoTuner$new( +autotuner_rf = mlr3tuning::auto_tuner( learner = lrn_rf, resampling = mlr3::rsmp("spcv_coords", folds = 5), # 空間分割 measure = mlr3::msr("regr.rmse"), # パフォーマンス測定 @@ -505,7 +500,7 @@ autotuner_rf = mlr3tuning::AutoTuner$new( ) ``` -`AutoTuner` -オブジェクトの `train()` -メソッドを呼び出すと、最終的にハイパーパラメータ\index{hyperparameter}のチューニングが実行され、指定したパラメータに対して最適なハイパーパラメータ\index{hyperparameter}の組み合わせが見つかる。 +`AutoTuner` -オブジェクトの `train()` -メソッドを呼び出すと、最終的にハイパーパラメータ\index{はいぱーぱらめーた@ハイパーパラメータ}のチューニングが実行され、指定したパラメータに対して最適なハイパーパラメータ\index{はいぱーぱらめーた@ハイパーパラメータ}の組み合わせが見つかる。 ```{r 15-eco-24, eval=FALSE, cache=TRUE, cache.lazy=FALSE} # ハイパーパラメータのチューニング @@ -525,16 +520,10 @@ autotuner_rf = readRDS("extdata/15-tune.rds") autotuner_rf$tuning_result ``` +### 予測地図 {#predictive-mapping} - - - - - -### 予測マッピング {#predictive-mapping} - -調整されたハイパーパラメータ\index{hyperparameter}は、これで予測に使用することができる。 -そのためには、フィットした `AutoTuner` オブジェクトの `predict` メソッドを実行するだけでよいのである。 +調整されたハイパーパラメータ\index{はいぱーぱらめーた@ハイパーパラメータ}は、これで予測に使用することができる。 +そのためには、適合した `AutoTuner` オブジェクトの `predict` メソッドを実行するだけでよい。 ```{r 15-eco-27, cache=TRUE, cache.lazy=FALSE, warning=FALSE} # 最適なハイパーパラメータの組み合わせで予測 @@ -592,45 +581,45 @@ all(values(pred - pred_2) == 0) ``` 予測地図には、はっきりとした植生帯 (Figure \@ref(fig:rf-pred)) が描かれている。 -**lomas** 山の植生帯の詳細については、@muenchow_soil_2013 を参照。 +**ロマ植生**山の植生帯の詳細については、@muenchow_soil_2013 を参照。 青い色調は、いわゆる *Tillandsia*-帯を表している。 -*Tillandsia* は、特に **lomas** 山脈の砂地やかなり砂漠的な麓で大量に見られる、高度に適応した属植物である。 +*Tillandsia* は、特に**ロマ植生**山脈の砂地やかなり砂漠的な麓で大量に見られる、高度に適応した属植物である。 黄色は草本植生帯で、*Tillandsia*-植生帯に比べ植物被度が高いことを表している。 オレンジ色はブロメリア帯を表し、種の豊富さと植物被覆率が最も高いことを特徴としている。 霧による湿度が最も高い気温逆転地帯 (標高約 750 - 850 m) の真下で見られる。 逆転温度以上になると当然水分は減少し、再び砂漠化し、数種の多肉植物が見られるようになる (多肉植物帯；赤色)。 -興味深いのは、空間予測によってブロメリア帯が途切れていることが明らかになったことである。これは、予測マッピングなしでは発見できなかった非常に興味深い発見であった。 +興味深いのは、空間予測によってブロメリア帯が途切れていることが明らかになったことである。これは、予測地図なしでは発見できなかった非常に興味深い発見であった。 ## 結論 -本章では、NMDS\index{NMDS} (Section \@ref(nmds)) を用いて、**lomas** の Mongón 山の群集行列を順序付け\index{ordination}した。 +本章では、NMDS\index{NMDS} (Section \@ref(nmds)) を用いて、**ロマ植生**の Mongón 山の群集行列を座標付け\index{ざひょうづけ@座標付け}した。 最初の軸は、調査地域の主な植物相の勾配を表し、部分的に R-GIS\index{GIS} ブリッジ ( Section \@ref(data-and-data-preparation) ) を使って導き出した環境予測因子の関数としてモデル化された。 -**mlr3**\index{mlr3 (package)} パッケージは、ハイパーパラメータ\index{hyperparameter} `mtry`、`sample.fraction` および `min.node.size` (Section \@ref(mlr3-building-blocks) ) を空間的に調整するためのビルディングブロックを提供している。 -調整されたハイパーパラメータ\index{hyperparameter}は最終モデルの入力となり、このモデルを環境予測変数に適用して植物相の勾配を空間的に表現した (Section \@ref(predictive-mapping))。 +**mlr3**\index{mlr3 (package)} パッケージは、ハイパーパラメータ\index{はいぱーぱらめーた@ハイパーパラメータ} `mtry`、`sample.fraction` および `min.node.size` (Section \@ref(mlr3-building-blocks) ) を空間的に調整するためのビルディングブロックを提供している。 +調整されたハイパーパラメータ\index{はいぱーぱらめーた@ハイパーパラメータ}は最終モデルの入力となり、このモデルを環境予測変数に適用して植物相の勾配を空間的に表現した (Section \@ref(predictive-mapping))。 その結果、砂漠の真ん中にある驚異的な生物多様性を空間的に示すことができたのである。 -**ロマス**山は絶滅の危機に瀕しているため、予測地図は保護区域を定める際の判断材料となり、地域住民に身近にあるユニークな存在であることを認識させることができるのである。 +**ロマ植生**山は絶滅の危機に瀕しているため、予測地図は保護区域を定める際の判断材料となり、地域住民に身近にあるユニークな存在であることを認識させることができるのである。 方法論の面では、いくつかの追加的な指摘ができる。 -- 2つ目の軸\index{ordination}もモデル化し、2つの軸のモデル化されたスコアを1つの予測地図に統合して可視化する革新的な方法を見つけるのは興味深いことである。 +- 2 つ目の軸\index{ざひょうづけ@座標付け}もモデル化し、2つの軸のモデル化されたスコアを1つの予測地図に統合して可視化する革新的な方法を見つけるのは興味深いことである。 - もし、生態学的に意味のある方法でモデルを解釈することに興味があれば、おそらく (セミ) パラメトリックモデルを使うべきだろう [@muenchow_predictive_2013;@zuur_mixed_2009;@zuur_beginners_2017] -しかし、少なくともランダムフォレスト\index{random forest}のような機械学習モデルの解釈を助けるアプローチは存在する (例えば、[https://mlr-org.com/posts/2018-04-30-interpretable-machine-learning-iml-and-mlr/](https://mlr-org.com/posts/2018-04-30-interpretable-machine-learning-iml-and-mlr/) を参照) -- 本章で使用したハイパーパラメータ\index{hyperparameter}のランダム化最適化よりも、逐次モデルベース最適化 (sequential model-based optimization, SMBO) の方が望ましいかもしれない [@probst_hyperparameters_2018] +しかし、少なくともランダムフォレスト\index{らんだむふぉれすと@ランダムフォレスト}のような機械学習モデルの解釈を助けるアプローチは存在する (例えば、[https://mlr-org.com/posts/2018-04-30-interpretable-machine-learning-iml-and-mlr/](https://mlr-org.com/posts/2018-04-30-interpretable-machine-learning-iml-and-mlr/) を参照) +- 本章で使用したハイパーパラメータ\index{はいぱーぱらめーた@ハイパーパラメータ}のランダム化最適化よりも、逐次モデルベース最適化 (sequential model-based optimization, SMBO) の方が望ましいかもしれない [@probst_hyperparameters_2018] -最後に、ランダムフォレスト\index{random forest}と他の機械学習\index{machine learning}モデルは、多くのオブザベーションと多くの予測因子、この章で使われるよりもはるかに多く、どの変数と変数の相互作用が応答を説明するのに寄与するかが不明である設定で頻繁に使用されることに注意しておこう。 +最後に、ランダムフォレスト\index{らんだむふぉれすと@ランダムフォレスト}と他の機械学習\index{きかいがくしゅう@機械学習}モデルは、多くのオブザベーションと多くの予測因子、この章で使われるよりもはるかに多く、どの変数と変数の相互作用が応答を説明するのに寄与するかが不明である設定で頻繁に使用されることに注意しておこう。 さらに、その関係は高度に非線形である可能性もある。 私たちのユースケースでは、レスポンスと予測変数の関係はかなり明確で、非線型はわずかであり、オブザベーションと予測変数の数は少ない。 -したがって、線形モデル\index{regression!linear}を試してみる価値はあるかもしれない。 -線形モデルは、ランダムフォレスト\index{random forest}モデルよりも説明や理解がしやすいので好まれ (パーシモンの法則)、さらに計算負荷が少ない (演習を参照)。 -線形モデルがデータに存在する非線形性の程度に対処できない場合、一般化加法モデル\index{generalized additive model} (generalized additive model, GAM) を試してみることもできる。 +したがって、線形モデル\index{かいき@回帰!せんけい@線形}を試してみる価値はあるかもしれない。 +線形モデルは、ランダムフォレスト\index{らんだむふぉれすと@ランダムフォレスト}モデルよりも説明や理解がしやすいので好まれ (パーシモンの法則)、さらに計算負荷が少ない (演習を参照)。 +線形モデルがデータに存在する非線形性の程度に対処できない場合、一般化加法モデル\index{いっぱんかかほうもでる@一般化加法モデル} (generalized additive model, GAM) を試してみることもできる。 ここで重要なのは、データサイエンティストのツールボックスは複数のツールで構成されており、目の前のタスクや目的に最適なツールを選択するのはあなたの責任であるということである。 -ここでは、ランダムフォレスト\index{random forest}のモデリングと、それに対応した空間予測図への利用方法を紹介したいと思われる。 +ここでは、ランダムフォレスト\index{らんだむふぉれすと@ランダムフォレスト}のモデリングと、それに対応した空間予測図への利用方法を紹介したいと思われる。 この目的のためには、反応と予測因子の関係が既知の、よく研究されたデータセットが適切である。 -しかし、これはランダムフォレスト\index{random forest}モデルが予測性能の面で最良の結果を返したことを意味するものではない。 +しかし、これはランダムフォレスト\index{らんだむふぉれすと@ランダムフォレスト}モデルが予測性能の面で最良の結果を返したことを意味するものではない。 ## 演習 ```{r, echo=FALSE, results='asis'} -res = knitr::knit_child('_15-ex.Rmd', quiet = TRUE, options = list(include = FALSE, eval = FALSE)) +res = knitr::knit_child('_15-ex-ja.Rmd', quiet = TRUE, options = list(include = FALSE, eval = FALSE)) cat(res, sep = '\n') ``` diff --git a/16-synthesis-ja.Rmd b/16-synthesis-ja.Rmd index d82a1ab..fbcdf2d 100644 --- a/16-synthesis-ja.Rmd +++ b/16-synthesis-ja.Rmd @@ -46,7 +46,7 @@ identical(nz_u1, nz_u3$geom) Base R の `aggregate()` 関数を使うか、**dplyr** の `summarise()` を使うかは好みの問題であるが、後者の方が読みやすいだろう。 つまり、R で地理データを扱う場合、たとえ 1 つのパッケージであっても、複数の選択肢から選ぶことができる場合が多いということである。 -R のパッケージが増えれば、さらに選択肢は広がる。例えば、古いパッケージの **sp** を使っても同じ結果を得ることができる。 +R のパッケージが増えれば、さらに選択肢は広がる。例えば、古いパッケージの **sp** を使っても同じ結果を得ることができる。\index{sp (package)} しかしながら、良いアドバイスをしたいという本書のゴールに従うと、パフォーマンスもよく将来性のある **sf** パッケージを推奨する。 このことは、本書のすべてのパッケージに当てはまるが、他の手段の存在を知って自分の選択について正当化できることは、(邪魔にならない程度に) 役に立つことがある。 @@ -113,22 +113,22 @@ length(spdeps) # 463 ## ギャップとオーバーラップ {#gaps} -ジオコンピュテーションは巨大な分野であり、多くのギャップがあることは避けられない。 +ジオコンピュテーションは巨大な分野であり、多くのギャップがあることは避けられない。\index{じおこんぴゅてーしょん@ジオコンピュテーション} 特定にトピックやテクニック、パッケージを意図的に強調し、あるいは省略するなど、選択的に行っている。 地理データの操作、座標参照系の基本、データの読み書き、可視化技術など、実際のアプリケーションで最もよく必要とされるトピックを重視するよう努めた。 また、本書は、ジオコンピュテーションに必要なスキルを身につけ、さらに高度なトピックや特定のアプリケーションに進む方法を紹介することを目的としており、いくつかのトピックやテーマは何度も登場する。 また、他で深く取り上げられているトピックをあえて省略した。 -例えば、点パターン解析、空間補間 (クリギング)、空間回帰といった空間データの統計的モデリングは、機械学習の文脈で Chapter \@ref(spatial-cv) で触れているが、詳細には触れていない。 +例えば、点パターン解析\index{てんぱたーんかいせき@点パターン解析}、空間補間\index{くうかんほかん@空間補間} (例えばクリギング)、空間回帰\inde{くうかんかいき@空間回帰}といった空間データの統計的モデリングは、機械学習の文脈で Chapter \@ref(spatial-cv) で触れているが、詳細には触れていない。 これらの手法については、@pebesma_spatial_2023 の統計学的指向の章や、ポイントパターン分析に関する書籍 [@baddeley_spatial_2015]、空間データに適用するベイズ手法 [@gomez-rubio_bayesian_2020;@moraga_spatial_2023]、健康 [@moraga_geospatial_2019] や[山火事深刻度分析](https://bookdown.org/mcwimberly/gdswr-book/application---wildfire-severity-analysis.html) など特定のアプリケーションに焦点を当てた書籍といった優れた資料がすでに存在している [@wimberly_geographic_2023]。 その他の話題としては、リモートセンシングや、GIS 専用ソフトと並行しての R の利用 (ブリッジとしてではなく) などが挙げられるが、これらは限定的である。 これらの話題については、[R におけるリモートセンシングについての議論](https://github.com/r-spatial/discuss/issues/56)、@wegmann_remote_2016 や [Marburg University](https://geomoer.github.io/moer-info-page/courses.html) から入手できる GIS 関連教材など、多くの資料がある。 -Chapter \@ref(spatial-cv) と Chapter \@ref(eco) で空間統計推論よりも機械学習に焦点を当てたのは、このトピックに関する質の高いリソースが豊富にあるためである。 +Chapter \@ref(spatial-cv) と Chapter \@ref(eco) で空間統計推論よりも機械学習に焦点を当てたのは、このトピックに関する質の高いリソースが豊富にあるためである。\index{とうけいすいろん@統計推論} これらのリソースには、生態系のユースケースに焦点を当てた @zuur_mixed_2009、@zuur_beginners_2017、そして [css.cornell.edu/faculty/dgr2](https://css.cornell.edu/faculty/dgr2/teach/) でホストされている *Geostatistics & Open-source Statistical Computing* の自由に利用できる教材とコードがある。 [*R for Geographic Data Science*](https://sdesabbata.github.io/r-for-geographic-data-science/) では、地理データサイエンスとモデリングのための R の紹介をしている。 -また、「ビッグデータ」(ハイスペックなラップトップに収まらないデータセットという意味) に対するジオコンピュテーションはほとんど省略している。 +また、「ビッグデータ」(ハイスペックなラップトップに収まらないデータセットという意味) に対するジオコンピュテーションはほとんど省略している。\index{びっぐでーた@ビッグデータ} この決定は、一般的な研究や政策アプリケーションに必要な地理的データセットの大部分は、個人用ハードウェアに収まるという事実によって正当化される (Section \@ref(cloud) 参照)。 コンピュータの RAM を増やしたり、[GitHub Codespaces: 本書のコードを実行可能](https://github.com/codespaces/new?hide_repo_select=true&ref=main&repo=84222786&machine=basicLinux32gb&devcontainer_path=.devcontainer.json&location=WestEurope) のようなプラットフォームで利用できる計算能力を一時的に「借りる」ことは可能である。 さらに、小さなデータセットで問題を解くことを学ぶことは、巨大なデータセットで問題を解くための前提条件であり、本書で強調しているのは「始めること」であり、ここで学んだスキルは大きなデータセットに移行したときに役立つものである。 @@ -187,13 +187,13 @@ Chapter \@ref(gis) で紹介した空間データベースは、メモリで処 - 大型・汎用プログラミングQ&Aサイト [stackoverflow.com](https://stackoverflow.com/) - [RStudio Community](https://community.rstudio.com/)、[rOpenSci Discuss](https://discuss.ropensci.org/) ウェブフォーラム、 [Stan](https://discourse.mc-stan.org/) フォーラムなど、特定のソフトウェアツールに関連するフォーラムなど、特定のエンティティに関連するオンラインフォーラム。 - GitHub などのソフトウェア開発プラットフォームは、R-spatial パッケージの大半の課題トラッカーや、最近では **sfnetworks** パッケージに関する議論 (バグ報告だけでなく) を促すために作られた議論ページなどをホストしている ([luukvdmeer/sfnetworks/discussions](https://github.com/luukvdmeer/sfnetworks/discussions/)を参照)。 -- チャットルームやフォーラム [rOpenSci](https://ropensci.org/blog/2022/09/13/contributing-ropensci/) や [geocompx](https://geocompx.org) コミュニティ (これには質問をすることができる [Discord server](https://discord.com/invite/PMztXYgNxp) もある)。本書もここに関係している。 +- チャットルームやフォーラム [rOpenSci](https://ropensci.org/blog/2022/09/13/contributing-ropensci/) や [geocompx](https://geocompx.org) コミュニティ (これには質問をすることができる [Discord server](https://discord.com/invite/PMztXYgNxp) もある)。本書もここに関係している。\index{geocompx} - [OSGeoJapan のメーリングリスト](https://www.osgeo.jp/mailing_list) (日本語) - [r-wakalang Slack](https://github.com/tokyor/r-wakalang) (日本語) ### **reprex** による再現性の例 {#reprex} -良い質問とは、明確に述べられた質問で、さらにアクセスしやすい完全に再現可能な例があるとよい (https://r4ds.hadley.nz/workflow-help.html も参照)。 +良い質問とは、明確に述べられた質問で、さらにアクセスしやすい完全に再現可能な例があるとよい (https://r4ds.hadley.nz/workflow-help.html も参照)。\index{さいげんせい@再現性} また、ユーザーの視点から「うまくいかなかった」コードを示した後、何を見たいかを説明することも有効である。 再現可能な例を作成するための非常に便利なツールが、**reprex** パッケージである。 予期せぬ動作を強調するために、問題を示す完全に再現可能なコードを書き、`reprex()`関数を使って、フォーラムや他のオンラインスペースに貼り付けられるようなコードのコピーを作成することができる。 @@ -256,7 +256,7 @@ Chapter \@ref(intro) で紹介されているように、ジオコンピュテ R の強みは、科学的な再現性を重視し、学術研究において広く使用され、地理データの統計的モデリングを比類なくサポートすることから、特にジオコンピュテーションの定義に関連しています。 さらに、他の言語やフレームワークを学ぶのはコンテキストスイッチというコストが発生するため、ジオコンピュテーションのための 1 つの言語を深く学ぶことを推奨する。 ] -*Python*、*C++*、*JavaScript*、*Scala*\index{Scala}、*Rust*\index{Rust} を使っても、同じ深さでジオコンピュテーションを勉強することができる。 +*Python*\index{Python}、*C++*、*JavaScript*、*Scala*\index{Scala}、*Rust*\index{Rust} を使っても、同じ深さでジオコンピュテーションを勉強することができる。 それぞれが進化した地理空間能力を有している。 例えば、Pythonのパッケージ [**rasterio**](https://github.com/rasterio/rasterio) は、この本で使われている **terra** パッケージを補足/置換できるものである。 Python のジオコンピュテーションについては、[*Geocomputation with Python*](https://py.geocompx.org/) を参照。 diff --git a/_11-ex-ja.Rmd b/_11-ex-ja.Rmd new file mode 100644 index 0000000..6689997 --- /dev/null +++ b/_11-ex-ja.Rmd @@ -0,0 +1,103 @@ +```{asis 11-ex-asis1, message=FALSE} +回答するには、以下のパッケージをアタッチすることとする (他のパッケージも必要に応じてアタッチする)。 +``` + +```{r setup11, include=FALSE} +poly_centroid = function(poly_mat) { + Origin = poly_mat[1, ] # create a point representing the origin + i = 2:(nrow(poly_mat) - 2) + T_all = lapply(i, function(x) {rbind(Origin, poly_mat[x:(x + 1), ], Origin)}) + C_list = lapply(T_all, t_centroid) + C = do.call(rbind, C_list) + A = vapply(T_all, t_area, FUN.VALUE = double(1)) + c(weighted.mean(C[, 1], A), weighted.mean(C[, 2], A)) +} +``` + +```{r 11-ex-e0} +library(sf) +``` + +E1. 本書の GitHub リポジトリの [`11-centroid-alg.R`](https://github.com/geocompx/geocompr/blob/main/code/11-centroid-alg.R) スクリプトを読みなさい。 + + - ベストプラクティスのうちどれを使っているか？ + - RStudio\index{RStudio} などの IDE\index{IDE} を使い、自分のパソコンでスクリプトを作成しなさい (スクリプトを 1 行ずつ打ち込み、適宜コメントを入れると良い。コピペはしない。この作業でスクリプトの入力の仕方を学ぶことができる。)。正方形ポリゴン (`poly_mat = cbind(x = c(0, 9, 9, 0, 0), y = c(0, 0, 9, 9, 0))` で作成) の例を使い、スクリプトを 1 行ずつ実行しなさい。 + - 再現可能性を高めるためにはどのように変更したら良いか? + - ドキュメンテーションをより良くするためにはどうしたら良いか? + +```{asis 11-ex-e1, message=FALSE} +スクリプトは論理的な場所に、適切なファイル名で保存されている。 +スクリプトはコメント付きできちんと文書化されており、コードの書式も整っている。 +スクリプトが再現可能である。 + +例えば、キーボードショートカット `Ctrl + Shift + N` (Windows) または `Cmd + Shift + N` (Mac) を使用して、`File > New File > R Script` をクリックするか、`Source` ペインの左上にある `+` アイコンをクリックします。 +また、R コンソールから `file.create("11-centroid-alg.R")` というコマンドで新しい R スクリプトを作成することもできる。 + +このスクリプトはすでに再現可能であり、`poly_mat` というオブジェクトが必要であることを示すメッセージが表示され、もし存在しなければ、テスト用のサンプルデータセットが最初に作成される。 +R を初めて使う人のために、スクリプトを実行する前に R をインストールする必要があることを示すコメントを含めることもできる。 + +ドキュメンテーションは、本の関連セクションへのリンクを含む、アルゴリズムのより詳細な説明で改善できるだろう。 +さらに、無名関数を名前付き関数に置き換え、Roxygen2 コメントで文書化することもできる。 +``` + + +E2. 幾何アルゴリズムのセクションで、`poly_mat` のポリゴンの面積は 245 で、重心は 8.8, 9.2 であると計算した。 + + - このアルゴリズムのスクリプトである [`11-centroid-alg.R`](https://github.com/geocompx/geocompr/blob/main/code/11-centroid-alg.R) を参照し、自分のパソコンで結果を再現しなさい (ボーナス: コピペせずに自分で入力しなさい)。 + - 結果は正しいか? `poly_mat` を `st_polygon()` 関数で `sfc` オブジェクトに変換し (`poly_sfc` という名前)、`st_area()` 関数と `st_centroid()` 関数を用いて検証しなさい (ヒント: この関数は、クラス `list()` を引数に取る)。 + +```{r 11-ex-e} +# We can verify the answer by converting `poly_mat` into a simple feature collection +# as follows, which shows the calculations match: +x_coords = c(10, 20, 12, 0, 0, 10) +y_coords = c(0, 15, 20, 10, 0, 0) +poly_mat = cbind(x_coords, y_coords) +poly_sfc = sf::st_polygon(list(poly_mat)) +sf::st_area(poly_sfc) +sf::st_centroid(poly_sfc) +# By calling the script: +# source("https://github.com/geocompx/geocompr/raw/main/code/11-centroid-alg.R") +``` + +E3. 我々が作成したアルゴリズム\index{あるごりずむ@アルゴリズム}は**凸包**に対してのみ動作すると記載されている。凸包を定義し (ジオメトリ操作の章を参照)、凸包でないポリゴンでアルゴリズムをテストしなさい。\index{とつほう@凸包} + +```{r 11-ex-e3} +x_coords = c(10, 20, 12, 0, 0, 5, 10) +y_coords = c(0, 15, 20, 10, 0, 5, 0) +plot(x_coords, y_coords, type = "l") +poly_mat = cbind(x_coords, y_coords) +# source("https://github.com/geocompx/geocompr/raw/main/code/11-centroid-alg.R") +# Area from our script: 270 +poly_sfc = sf::st_polygon(list(poly_mat)) +sf::st_area(poly_sfc) # Actual area: 220 +``` + + - ボーナス 1: なぜこの方法が凸の外皮に対してのみ機能するのかを考え、他の種類の多角形に対して機能させるためにアルゴリズムに加える必要がある変更点に注意する。 + - ボーナス 2: `11-centroid-alg.R` の内容を基に、行列形式で表現された線分の全長を求めることができる、Base R 関数のみを使ったアルゴリズムを書きなさい。 + + + +```{asis 11-ex-e3-bonus1} +このアルゴリズムは、正の面積値だけでなく負の面積値も持つことができる必要がある。 + +ボーナス 2 は読者のための練習として残しておく。 +``` + +E4. 関数のセクションでは、`sfg` クラスの出力 (`poly_centroid_sfg()`) と `matrix` 型の出力 (`poly_centroid_type_stable()`) を生成する `poly_centroid()` 関数の異なるバージョンを作成した。 +さらに、型が安定で (`sf` クラスの入力しか受け付けない) `sf` オブジェクトを返すバージョン (例えば `poly_centroid_sf()`) を作成し、関数を拡張しなさい (ヒント: `sf::st_coordinates(x)` コマンドでオブジェクト `x` を行列に変換する必要があるかもしれない)。 + + - `poly_centroid_sf(sf::st_sf(sf::st_sfc(poly_sfc)))` を実行し、動作するか検証しなさい + - `poly_centroid_sf(poly_mat)` を実行しようとした時、どのようなエラーメッセージが表示されたか? + +```{r 11-ex-e4} +poly_centroid_sf = function(x) { + stopifnot(is(x, "sf")) + xcoords = sf::st_coordinates(x) + centroid_coords = poly_centroid(xcoords) + centroid_sf = sf::st_sf(geometry = sf::st_sfc(sf::st_point(centroid_coords))) + centroid_sf +} +poly_centroid_sf(sf::st_sf(sf::st_sfc(poly_sfc))) +poly_centroid_sf(poly_sfc) +poly_centroid_sf(poly_mat) +``` diff --git a/_12-ex-ja.Rmd b/_12-ex-ja.Rmd new file mode 100644 index 0000000..0df3d6b --- /dev/null +++ b/_12-ex-ja.Rmd @@ -0,0 +1,237 @@ +```{asis 12-ex-asis1, message=FALSE} +回答するには、以下のパッケージをアタッチすることとする (他のパッケージも必要に応じてアタッチする)。 +``` + +```{r 12-ex-e0, message=FALSE, warning=FALSE, eval=FALSE} +library(dplyr) +# library(kernlab) +library(mlr3) +library(mlr3learners) +library(mlr3extralearners) +library(mlr3spatiotempcv) +library(mlr3tuning) +library(qgisprocess) +library(terra) +library(sf) +library(tmap) +``` + +E1. `terra::rast(system.file("raster/ta.tif", package = "spDataLarge"))$elev` で読み込んだ `elev` データセットから、R-GIS ブリッジ (GIS ソフトウェアへのブリッジの章を参照) を用いて以下の地形属性を計算しなさい。 + + - 傾斜角度 + - 平面曲率 + - プロファイル曲率 + - 集水域 + +```{r 12-ex-e1-1, eval=FALSE} +# データをアタッチ +dem = terra::rast(system.file("raster/ta.tif", package = "spDataLarge"))$elev + +algs = qgisprocess::qgis_algorithms() +qgis_search_algorithms("curvature") +alg = "sagang:slopeaspectcurvature" +qgisprocess::qgis_show_help(alg) +qgisprocess::qgis_get_argument_specs(alg) +# terrain attributes (ta) +out_nms = paste0(tempdir(), "/", c("slope", "cplan", "cprof"), + ".sdat") +args = rlang::set_names(out_nms, c("SLOPE", "C_PLAN", "C_PROF")) +out = qgis_run_algorithm(alg, ELEVATION = dem, METHOD = 6, + UNIT_SLOPE = "[1] degree", + !!!args, + .quiet = TRUE + ) +ta = out[names(args)] |> unlist() |> terra::rast() +names(ta) = c("slope", "cplan", "cprof") +# catchment area +qgis_search_algorithms("[Cc]atchment") +alg = "sagang:catchmentarea" +qgis_show_help(alg) +qgis_get_argument_specs(alg) +carea = qgis_run_algorithm(alg, + ELEVATION = dem, + METHOD = 4, + FLOW = file.path(tempdir(), "carea.sdat")) +# transform carea +carea = terra::rast(carea$FLOW[1]) +log10_carea = log10(carea) +names(log10_carea) = "log10_carea" +# add log_carea and dem to the terrain attributes +ta = c(ta, dem, log10_carea) +``` + +E2. `slope`、`cplan`、`cprof`、`elev`、`log_carea` という新しい変数を追加し、対応する出力ラスタから `lsl` データフレーム (`data("lsl", package = "spDataLarge"`)) に値を抽出しなさい。 + +```{r 12-ex-e2, eval=FALSE} +# +# terrain attribute raster stack をアタッチ (前の演習で行っていない場合) +data("lsl", package = "spDataLarge") +ta = terra::rast(system.file("raster/ta.tif", package = "spDataLarge")) +lsl = select(lsl, x, y, lslpts) +# 値を点に抽出、predictor を作成 +lsl[, names(ta)] = terra::extract(ta, lsl[, c("x", "y")]) |> + select(-ID) +``` + +E3. 導き出された地形属性ラスタを GLM と組み合わせて、Figure 12.2に示すような空間予測マップを作成しなさい。 +`data("study_mask", "package="spDataLarge")` を実行すると、調査地域のマスクが添付される。 + +```{r 12-ex-e3, eval=FALSE} +# データをアタッチ (E1 と E2 で行っていない場合) +# landslide points with terrain attributes and terrain attribute raster stack +data("lsl", "study_mask", package = "spDataLarge") +ta = terra::rast(system.file("raster/ta.tif", package = "spDataLarge")) + +# モデルに適合 +fit = glm(lslpts ~ slope + cplan + cprof + elev + log10_carea, + data = lsl, family = binomial()) + +# 予測を作成 +pred = terra::predict(object = ta, model = fit, type = "response") + +# 地図を作成 +lsl_sf = sf::st_as_sf(lsl, coords = c("x", "y"), crs = 32717) +lsl_sf = sf::st_as_sf(lsl, coords = c("x", "y"), crs = 32717) +hs = terra::shade(ta$slope * pi / 180, + terra::terrain(ta$elev, v = "aspect", unit = "radians")) +rect = tmaptools::bb_poly(raster::raster(hs)) +bbx = tmaptools::bb(raster::raster(hs), xlim = c(-0.00001, 1), + ylim = c(-0.00001, 1), relative = TRUE) + +tm_shape(terra::mask(hs, study_mask), bbox = bbx) + + tm_grid(col = "black", n.x = 1, n.y = 1, labels.inside.frame = FALSE, + labels.rot = c(0, 90), lines = FALSE) + + tm_raster(col.scale = tm_scale(values = gray(0:100 / 100), n = 100), + col.legend = tm_legend_hide()) + + # 予測ラスタ + tm_shape(terra::mask(pred, study_mask)) + + tm_raster(col_alpha = 0.5, + col.scale = tm_scale(values = "Reds", n = 6), + col.legend = tm_legend(title = "Susceptibility")) + + # 矩形と外側マージン + tm_shape(rect) + + tm_borders() + + tm_layout(legend.position = c("left", "bottom"), + legend.title.size = 0.9) +``` + +E4. GLM 学習器に基づき、100 回繰り返した 5 フォールドの非空間クロスバリデーションと空間 CV を計算し、箱ひげ図を用いて両方のリサンプリング戦略からの AUROC 値を比較しなさい。 + +ヒント: 非空間リサンプリング戦略を指定する必要がある。 + +追加ヒント: `mlr3::benchmark()` と `mlr3::benchmark_grid()` を使って、練習問題 4 から 6 を一度に解くことができます (詳しくは https://mlr3book.mlr-org.com/chapters/chapter10/advanced_technical_aspects_of_mlr3.html#sec-fallback を参照)。 +その際、計算には非常に時間がかかり、おそらく数日かかることを覚悟しよう。 +もちろん、これはシステムに依存する。 +自由に使える RAM とコアが多ければ多いほど、計算時間は短くなる。 + +```{r 12-ex-e4, eval=FALSE} +# データをアタッチ (E1 と E2 で行っていない場合) +data("lsl", package = "spDataLarge") # landslide points with terrain attributes + +# task を作成 +task = TaskClassifST$new( + id = "lsl_ecuador", + backend = mlr3::as_data_backend(lsl), target = "lslpts", positive = "TRUE", + coordinate_names = c("x", "y"), + coords_as_features = FALSE, + crs = 32717 +) + +# 学習機を準備 (すべての演習で使用) +# GLM +lrn_glm = lrn("classif.log_reg", predict_type = "prob") +lrn_glm$fallback = lrn("classif.featureless", predict_type = "prob") + +# SVM +# SVM 学習機を準備 (using ksvm function from the kernlab package) +lrn_ksvm = lrn("classif.ksvm", predict_type = "prob", kernel = "rbfdot", + type = "C-svc") +lrn_ksvm$fallback = lrn("classif.featureless", predict_type = "prob") + +# ネストされたリサンプリングを特定し学習器を調整 +# 5 つの空間的に非連続な分割 +tune_level = rsmp("spcv_coords", folds = 5) +# ランダムに選択した 50 のハイパーパラメータ +terminator = trm("evals", n_evals = 50) +tuner = tnr("random_search") +# ランダムに選択したハイパーパラメータの外側を定義 +ps = ps( + C = p_dbl(lower = -12, upper = 15, trafo = function(x) 2^x), + sigma = p_dbl(lower = -15, upper = 6, trafo = function(x) 2^x) +) +at_ksvm = AutoTuner$new( + learner = lrn_ksvm, + resampling = tune_level, + measure = msr("classif.auc"), + search_space = ps, + terminator = terminator, + tuner = tuner +) + +# QDA +lrn_qda = lrn("classif.qda", predict_type = "prob") +lrn_qda$fallback = lrn("classif.featureless", predict_type = "prob") + +# ハイパーパラメータのチューニングなしで SVM +vals = lrn_ksvm$param_set$values +lrn_ksvm_notune = lrn_ksvm$clone() +lrn_ksvm_notune$param_set$values = c(vals, C = 1, sigma = 1) + +# リサンプリング戦略を定義 +# リサンプリング方法を指定 例 空間 CV 100 繰り返し 5 fold +# -> 各繰り返しで、データセットを 5 倍に分割 +# 方法: repeated_spcv_coords -> 空間分割 +rsmp_sp = rsmp("repeated_spcv_coords", folds = 5, repeats = 100) +# 方法: repeated_cv -> 非空間分割 +rsmp_nsp = rsmp("repeated_cv", folds = 5, repeats = 100) + +# (空間) クロスバリデーション +#**************************** +# デザインを作成 +grid = benchmark_grid(tasks = task, + learners = list(lrn_glm, at_ksvm, lrn_qda, + lrn_ksvm_notune), + resamplings = list(rsmp_sp, rsmp_nsp)) +# クロスバリデーションを実行 +library(future) +# 外側ループを順次実行し、内側ループを並列化 +future::plan(list("sequential", "multisession"), + workers = floor(availableCores() / 2)) +set.seed(021522) +bmr = benchmark(grid, + store_backends = FALSE, + store_models = FALSE, + encapsulate = "evaluate") +# 並列化を終了 +future:::ClusterRegistry("stop") +# 結果を保存 +# saveRDS(bmr, file = "extdata/12-bmr.rds") + +# 結果をプロット +autoplot(bmr, measure = msr("classif.auc")) +``` + +E5. 二次判別分析 (quadratic discriminant analysis, QDA) を用いて地すべり感受性をモデル化。 +QDA の予測性能を評価する。 +QDA と GLM の空間クロスバリデーション平均 AUROC 値の差は? + +```{r 12-ex-e5, eval=FALSE} +# データをアタッチ (E4 で行っていない場合) +bmr = readRDS("extdata/12-bmr.rds") + +# 結果をプロット +autoplot(bmr, measure = msr("classif.auc")) +# QDA は GLM よりも平均して AUROC 値が高く、これは中程度であることを示している +# non-linear boundaries +``` + +E6. ハイパーパラメータを調整せずに SVM を実行する。 +`rbfdot` カーネルを $sigma$ = 1 と *C* = 1 で使用する。 +**kernlab** の `ksvm()` でハイパーパラメータを指定しないままにしておくと、自動的に非空間的なハイパーパラメータチューニングが初期化される。 + +```{r 12-ex-e6, eval=FALSE} +# データをアタッチ (E4 で行っていない場合) +bmr = readRDS("extdata/12-bmr.rds") +# 結果をプロット +autoplot(bmr, measure = msr("classif.auc")) +``` diff --git a/_13-ex-ja.Rmd b/_13-ex-ja.Rmd new file mode 100644 index 0000000..1785e53 --- /dev/null +++ b/_13-ex-ja.Rmd @@ -0,0 +1,72 @@ +```{r 13-ex-e0, message=FALSE} +library(sf) +library(spDataLarge) +``` + +E1. 本章で紹介した分析の多くでは、アクティブ (自転車のこと) という交通モードに焦点を当てたが、車でのトリップ についてはどうだろうか? + + - `desire_lines` オブジェクトに含まれるトリップのうち、車でのトリップの割合は? + - 直線距離が 5 km 以上の`desire_lines`の割合は? + - 長さが 5 km 以上の希望線に含まれるトリップのうち、車で移動するトリップの割合は? + - 長さが 5 km 未満で、移動の 50% 以上が車である希望線をプロットする。 + - これらの自動車に依存しながらも短い希望線の位置について、何か気づくことはあるか? + +```{r 13-e1, eval=FALSE, echo=FALSE} +sum(desire_lines$car_driver) / sum(desire_lines$all) +# 57% +desire_lines_5km_plus = desire_lines |> + filter(distance_km > 5) +# Just over are half ar 5km+, 54%: +nrow(desire_lines_5km_plus) / nrow(desire_lines) +# 71 of 5km+ trips are made by car +sum(desire_lines_5km_plus$car_driver) / sum(desire_lines_5km_plus$all) + +desire_lines_driving = desire_lines |> + mutate(`Proportion driving` = car_driver / all) |> + filter(`Proportion driving` > 0.5) +nrow(desire_lines_5km_plus_driving) / nrow(desire_lines) + +desire_lines_5km_less_50_pct_driving = desire_lines |> + filter(distance_km <= 5) |> + mutate(`Proportion driving` = car_driver / all) |> + filter(`Proportion driving` > 0.5) +desire_lines_5km_less_50_pct_driving |> + tm_shape() + + tm_lines("Proportion driving") +``` + +E2. 最後の図に示されたすべてのルート (既存の自転車道から 100 m 以上離れた区間) が建設された場合、自転車道の長さはどの程度増加するか? + +```{r 13-transport-29, eval=FALSE, echo=FALSE} +sum(st_length(route_network_no_infra)) +# 104193.6 [m] +# Just over 100 km +``` + +E3. `desire_lines` に含まれるトリップのうち、`routes_short_scenario` に含まれるトリップの割合はいくらか? + + - ボーナス: 全トリップのうち、`routes_short_scenario` を横切る希望線の割合は？ + +```{r 13-transport-30, echo=FALSE, eval=FALSE} +sum(routes_short_scenario$all) / sum(desire_lines$all) # 13% +d_intersect = desire_lines[routes_short_scenario, , op = st_crosses] +sum(d_intersect$all) / sum(desire_lines$all) # 88% +``` + +E4. 本章で紹介する分析は、ジオコンピュテーションの手法をどのように交通研究に応用できるかを教えるためのものである。 +実際に政府機関や交通コンサルタント会社でこのようなことをする場合、どの点が変わるだろうか? 大きいものから 3 点述べなさい。 + +```{r} +# Higher level of geographic resolution. +# Use cycle-specific routing services. +# Identify key walking routes. +# Include a higher proportion of trips in the analysis +``` + +E5. 最後の図で特定されたルートは、明らかに、全体像の一部を示しているに過ぎない。 +どのように分析を拡張するか? + +E6. カーフリーゾーン、駐輪ポイント、減車戦略など、場所ベースのサイクリング政策に投資するための主要な**エリア** (ルートではない) を作成することによって、シナリオを拡張したいと想像する。 +ラスタ\index{らすた@ラスタ}データセットは、この作業をどのように支援できるか? + + - ボーナス: Bristol 地域を 100 のセル（10 × 10）に分割し、それぞれの道路の平均制限速度を `bristol_ways` データセットから推定するラスタレイヤを開発しなさい (Chapter \@ref(location) 参照)。 diff --git a/_14-ex-ja.Rmd b/_14-ex-ja.Rmd new file mode 100644 index 0000000..982ce21 --- /dev/null +++ b/_14-ex-ja.Rmd @@ -0,0 +1,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") +``` diff --git a/_15-ex-ja.Rmd b/_15-ex-ja.Rmd new file mode 100644 index 0000000..85a5376 --- /dev/null +++ b/_15-ex-ja.Rmd @@ -0,0 +1,247 @@ +回答するには、以下のパッケージをアタッチすることとする (他のパッケージも必要に応じてアタッチする)。 + +```{r 15-ex-e0, message=FALSE, warning=FALSE, eval=FALSE} +library(sf) +library(terra) +library(data.table) +library(dplyr) +library(future) +library(ggplot2) +library(lgr) +library(mlr3) +library(mlr3learners) +library(mlr3spatiotempcv) +library(mlr3tuning) +library(mlr3viz) +library(progressr) +library(qgisprocess) +library(tictoc) +library(vegan) +``` + +E1. コミュニティ行列のパーセンテージデータを使用して、NMDS\index{NMDS} を実行する。 +ストレス値を報告し、存在-不在データを使用して NMDS から取得したストレス値と比較します。 +この違いを説明するものは何か? + +```{r 15-ex-e1, message=FALSE, eval=FALSE} +data("comm", package = "spDataLarge") +pa = vegan::decostand(comm, "pa") +pa = pa[rowSums(pa) != 0, ] +comm = comm[rowSums(comm) != 0, ] +set.seed(25072018) +nmds_pa = vegan::metaMDS(comm = pa, k = 4, try = 500) +nmds_per = vegan::metaMDS(comm = comm, k = 4, try = 500) +nmds_pa$stress +nmds_per$stress +``` + +```{asis 15-ex-e1-asis, message=FALSE} +存在-不在値を使った NMDS は、パーセントを使ったもの `nmds_per$stress`) よりも良い結果 (`nmds_pa$stress`) を出した。 +これは一見意外に思えるかもしれない。 +一方、パーセント行列はより多くの情報とノイズの両方を含んでいる。 +もう一つの側面は、データの収集方法である。 +フィールドにいる植物学者を想像してみてほしい。 +被度 5% の植物と被度 10% の植物を区別することは可能だと思えるかもしれない。 +しかし、3 回しか検出されず、その結果被覆率が 0.0001% など非常に小さい草本種はどうだろう。 +別の草本種が 6 回検出されたとして、そのカバー率は 0.0002% だろうか？ +ここで重要なのは、フィールドキャンペーン中に指定されたパーセントデータは、データにない精度を反映している可能性があるということである。 +これはまたノイズをもたらし、順序付けの結果を悪化させる。 +それでも、ある種が他のプロットよりも高い頻度やカバー率を持っている場合は、単なる存在-不在データと比べて貴重な情報である。 +妥協案としては、ロンド尺度のようなカテゴリ尺度を使用することである。 +``` + +E2. この章で使用したすべての予測ラスタ (集水勾配、集水面積) を計算し、`SpatRaster`オブジェクトに格納しなさい。 +そこに `dem` と `ndvi` を追加しなさい。\index{らすた@ラスタ} +次に、プロファイルと接線曲率を計算し、追加の予測ラスタとして追加しなさい (ヒント: `grass7:r.slope.aspect`)。 +最後に、応答予測行列を構築しなさい。 +最初の NMDS\index{NMDS} 軸のスコア (存在-不在コミュニティ行列を使用したときの結果) を標高に従って回転させたものが応答変数を表し、`random_points`に結合とする (内側結合を使用する)。 +応答予測行列を完成させるために、環境予測ラスタ・オブジェクトの値を `random_points` に抽出しなさい。 + + +```{r 15-ex-e2, eval = FALSE} +# まず、本章で使用した terrain 属性を計算 +library(dplyr) +library(terra) +library(qgisprocess) +library(vegan) +data("comm", "random_points", package = "spDataLarge") +dem = terra::rast(system.file("raster/dem.tif", package = "spDataLarge")) +ndvi = terra::rast(system.file("raster/ndvi.tif", package = "spDataLarge")) + +# 存在-不在行列を使う、空を削除 +pa = vegan::decostand(comm, "pa") +pa = pa[rowSums(pa) != 0, ] + +# プラグインを使用 +qgisprocess::qgis_enable_plugins(c("grassprovider", "processing_saga_nextgen")) + +# 環境予測因子 (ep) 集水傾斜、集水域を計算する。 +ep = qgisprocess::qgis_run_algorithm( + alg = "sagang:sagawetnessindex", + DEM = dem, + SLOPE_TYPE = 1, + SLOPE = tempfile(fileext = ".sdat"), + AREA = tempfile(fileext = ".sdat"), + .quiet = TRUE) +# 集水域、集水傾斜を読み込む +ep = ep[c("AREA", "SLOPE")] |> + unlist() |> + terra::rast() +# 列名を変更 +names(ep) = c("carea", "cslope") +# すべてのラスタを同じオリジンに揃える +origin(ep) = origin(dem) +# dem と ndvi を multilayer SpatRaster オブジェクトに追加 +ep = c(dem, ndvi, ep) +ep$carea = log10(ep$carea) + +# 曲率を計算 +qgis_show_help("grass7:r.slope.aspect") +curvs = qgis_run_algorithm( + "grass7:r.slope.aspect", + elevation = dem, + .quiet = TRUE) +# 曲率を ep に追加 +curv_nms = c("pcurvature", "tcurvature") +curvs = curvs[curv_nms] |> + unlist() |> + terra::rast() +curvs = terra::app(curvs, as.numeric) +names(curvs) = curv_nms +ep = c(ep, curvs) +random_points[, names(ep)] = + # terra::extract は ID 列を追加するが、不要 + terra::extract(ep, random_points) |> + select(-ID) +elev = dplyr::filter(random_points, id %in% rownames(pa)) %>% + dplyr::pull(dem) +# NMDSを高度に応じて回転させる (湿度の代理) +rotnmds = MDSrotate(nmds_pa, elev) +# 最初の 2 軸を抽出 +sc = vegan::scores(rotnmds, choices = 1:2, display = "sites") +rp = data.frame(id = as.numeric(rownames(sc)), + sc = sc[, 1]) +# 予測因子 (dem, ndvi, terrain attributes) を結合 +rp = dplyr::inner_join(random_points, rp, by = "id") +# saveRDS(rp, "extdata/15-rp_exercises.rds") +``` + +E3. 空間交差検証\index{くろすばりでーしょん@クロスバリデーション!くうかんてき@空間的 CV}を使用して、ランダムフォレスト\index{らんだむふぉれすと@ランダムフォレスト}と線形モデルのバイアス削減 RMSE を取得しなさい。 +ランダムフォレストのモデリングには、最適なハイパーパラメータ\index{はいぱーぱらめーた@ハイパーパラメータ}の組み合わせの推定 (50 回の反復によるランダム探索) を内部チューニングループに含めなさい。 +チューニングレベルを並列化しなさい。 +平均 RMSE\index{RMSE} を報告し、箱ひげ図を使用して、検索されたすべての RMSE を可視化しなさい。 +この課題は、mlr3の関数 `benchmark_grid()` と `benchmark()` (詳細は https://mlr3book.mlr-org.com/perf-eval-cmp.html#benchmarking を参照) を用いて解くのが最適。 + +```{r 15-ex-e3, message=FALSE, eval=FALSE} +library(dplyr) +library(future) +library(mlr3) +library(mlr3spatiotempcv) +library(mlr3learners) +library(mlr3viz) +library(paradox) +library(ranger) + +# 前の演習を行っていない場合、以下を実行 +# rp = readRDS("extdata/15-rp_exercises.rds") + +# task を定義 +task = mlr3spatiotempcv::as_task_regr_st( + select(rp, -id, -spri), + target = "sc", + id = "mongon") + +# 学習器を定義 +mlr3::mlr_learners +# 線形モデル +lrn_lm = mlr3::lrn("regr.lm", predict_type = "response") +# ランダムフォレスト +lrn_rf = mlr3::lrn("regr.ranger", predict_type = "response") +# ランダムフォレストの AutoTuner を定義 +search_space = paradox::ps( + mtry = paradox::p_int(lower = 1, upper = ncol(task$data()) - 1), + sample.fraction = paradox::p_dbl(lower = 0.2, upper = 0.9), + min.node.size = paradox::p_int(lower = 1, upper = 10) +) +at_rf = mlr3tuning::AutoTuner$new( + learner = lrn_rf, + # 空間分割 + resampling = mlr3::rsmp("spcv_coords", folds = 5), + # パフォーマンス計測 + measure = mlr3::msr("regr.rmse"), + search_space = search_space, + # 50 回繰り返しのランダム探索 + terminator = mlr3tuning::trm("evals", n_evals = 50), + tuner = mlr3tuning::tnr("random_search") +) +# リサンプリング戦略を定義 +rsmp_sp = mlr3::rsmp("repeated_spcv_coords", folds = 5, repeats = 100) + +# ベンチマーク設計を作成 +design_grid = mlr3::benchmark_grid( + tasks = task, + learners = list(lrn_lm, at_rf), + resamplings = rsmp_sp) +print(design_grid) +# 外部ループを順次実行し、内部ループを並列化 +future::plan(list("sequential", "multisession"), + workers = floor(future::availableCores() / 2)) +set.seed(10112022) +# verbosity を下げる +lgr::get_logger("mlr3")$set_threshold("warn") +lgr::get_logger("bbotk")$set_threshold("info") +# 注意: ベンチマークの実行には時間がかかる +# よって、結果は以下のファイルに保存しておく +# extdata/15-bmr-exercises.rds (下記参照) +tictoc::tic() +progressr::with_progress(expr = { + bmr = mlr3::benchmark( + design = design_grid, + # `resample()` と `benchmark()` の引数 `encapsulate` で、カプセル化 + # かつフォールバック学習器をそれぞれの特徴なし学習器に設定する + # これは簡便性のためだけである + # それぞれの学習器を設定したりより細かく設定することも可能 + # + encapsulate = "evaluate", + store_backends = FALSE, + store_models = FALSE) +}) +tictoc::toc() + +# 並列化を終了 +future:::ClusterRegistry("stop") +# すでの結果は保存済み +# saveRDS(bmr, file = "extdata/15-bmr_exercises.rds") +# 非常に時間がかかるため、自分で空間 CV を実行したくない場合には +# これを読み込む +# bmr = readRDS("extdata/15-bmr_exercises.rds") + +# 平均 RMSE +bmr$aggregate(measures = msr("regr.rmse")) +# あるいは、計算する +agg = bmr$aggregate(measures = msr("regr.rmse")) + +# 平均 rmse を考慮すると、lm の方がわずかに良い +purrr::map(agg$resample_result, ~ mean(.$score(msr("regr.rmse"))$regr.rmse)) +# 中央値を見ると、ランダムフォレストの方がわずかに良い +purrr::map(agg$resample_result, ~ median(.$score(msr("regr.rmse"))$regr.rmse)) + +# 箱ひげ図を作成 (autoplot を使うには mlr3viz をアタッチする!) +library(mlr3viz) +autoplot(bmr, measure = msr("regr.rmse")) + +# 自動で行わない場合 +# AUROC 値を抽出し、data.table に保存 +d = purrr::map_dfr(agg$resample_result, ~ .$score(msr("regr.rmse"))) +# 箱ひげ図を作成 +library(ggplot2) +ggplot(data = d, mapping = aes(x = learner_id, y = regr.rmse)) + + geom_boxplot(fill = c("lightblue2", "mistyrose2")) + + theme_bw() + + labs(y = "RMSE", x = "model") +``` + +```{asis 15-ex-e3-asis, message=FALSE} +実際、`lm`は、少なくともランダムフォレスト・モデルと同程度のパフォーマンスを示す。したがって、はるかに理解しやすく、計算負荷が少ない (ハイパーパラメータの適合が不要) ので、好まれるべきである。 +しかし、使用したデータ集合は，オブザベーションと予測変数が小さくm応答-予測変数の関係も比較的線形であることに留意してほしい． +```