From e1e9ed6d4b9c79ad5208c90c4a38b3ae9b671f86 Mon Sep 17 00:00:00 2001 From: babayoshihiko Date: Fri, 19 Jan 2024 16:24:52 +0900 Subject: [PATCH] 17 Jan 2024 --- 05-geometry-operations-ja.Rmd | 18 +-- 09-mapping-ja.Rmd | 223 ++++++++++++++++----------------- 10-gis-ja.Rmd | 228 ++++++++++++++++------------------ 11-algorithms-ja.Rmd | 27 ++-- 4 files changed, 239 insertions(+), 257 deletions(-) diff --git a/05-geometry-operations-ja.Rmd b/05-geometry-operations-ja.Rmd index d071e42..9041eab 100644 --- a/05-geometry-operations-ja.Rmd +++ b/05-geometry-operations-ja.Rmd @@ -529,17 +529,17 @@ polyg = st_cast(multipoint, "POLYGON") ``` ```{r single-cast, echo = FALSE, fig.cap="多点ジオメトリからキャストされた線とポリゴンの例。", warning=FALSE, fig.asp=0.3, fig.scap="Examples of casting operations."} -p_sc1 = tm_shape(st_sfc(multipoint)) + tm_symbols(shape = 1, - col = "black", +p_sc1 = tm_shape(st_sfc(multipoint, crs = "+proj=merc")) + tm_symbols(shape = 1, + col = "black", size = 0.5) + tm_title("複合点") + - tm_layout(inner.margins = c(0.15, 0.05, 0.15, 0.05), fontfamily = "HiraginoSans-W3") -p_sc2 = tm_shape(st_sfc(linestring)) + tm_lines() + - tm_title("線") -+ tm_layout(inner.margins = c(0.15, 0.05, 0.15, 0.05), fontfamily = "HiraginoSans-W3") -p_sc3 = tm_shape(st_sfc(polyg)) + tm_polygons(border.col = "black") + - tm_title("ポリゴン") -+ tm_layout(inner.margins = c(0.15, 0.05, 0.15, 0.05), fontfamily = "HiraginoSans-W3") + tm_layout(inner.margins = c(0.15, 0.05, 0.15, 0.05)) +p_sc2 = tm_shape(st_sfc(linestring, crs = "+proj=merc")) + tm_lines() + + tm_title("複合線") + + tm_layout(inner.margins = c(0.15, 0.05, 0.15, 0.05)) +p_sc3 = tm_shape(st_sfc(polyg, crs = "+proj=merc")) + tm_polygons(border.col = "black") + + tm_title("ポリゴン") + + tm_layout(inner.margins = c(0.15, 0.05, 0.15, 0.05)) tmap_arrange(p_sc1, p_sc2, p_sc3, ncol = 3) ``` diff --git a/09-mapping-ja.Rmd b/09-mapping-ja.Rmd index 86ae70c..7a4b2d1 100644 --- a/09-mapping-ja.Rmd +++ b/09-mapping-ja.Rmd @@ -2,6 +2,10 @@ # R で地図を作る {#adv-map} +```{r, include=FALSE} +source("code/before_script.R") +``` + ## 必須パッケージ {- #prerequisites-09} - この章では、すでに使用している以下のパッケージが必要である。 @@ -14,10 +18,18 @@ library(spData) library(spDataLarge) ``` -- また、以下の可視化パッケージを使用する (インタラクティブな地図アプリケーションを開発する場合は shiny もインストールしておこう)。 +- 本章での主要なパッケージは **tmap** である。 +CRAN よりも頻繁に更新されている [r-universe](https://r-universe.dev/) 版をお勧めする。訳註: macOS では、**tmap** は文字化けすることがある。[macOS で ragg を使用することで文字化けが解消する](https://uribo.hatenablog.com/entry/2021/03/29/202756)。 + +```{r} +#| eval: false +install.packages("tmap", repos = c("https://r-tmap.r-universe.dev", + "https://cloud.r-project.org")) +``` + +- 以下の可視化に関するパッケージを使用する (動的な地図アプリを開発したい場合は、**shiny** もインストールしよう)。 ```{r 08-mapping-2, message=FALSE} -remotes::install_github("r-tmap/tmap") library(tmap) # 静的地図と動的地図 library(leaflet) # 動的地図 library(ggplot2) # tidyverse データ可視化パッケージ @@ -79,9 +91,6 @@ Base R のアプローチは拡張可能で、`plot()` は何十もの引数を また、`tmap_mode()` を介して、同じコードで静的な地図とインタラクティブな地図を生成するユニークな機能を備えている。 最後に、(**sf** オブジェクトと **terra** オブジェクトを含む) 空間クラスを受け入れることができる点では、**ggplot2** などよりも優れている。 - - - ### tmap の基礎知識 {#tmap-basics} \index{tmap (package)!basics} @@ -122,16 +131,8 @@ source("https://github.com/Robinlovelace/geocompr/raw/main/code/09-tmshape.R", p - `tm_raster()`: ラスタデータの色付きのセル (３レイヤのあるラスタには `tm_rgb()` もある) - `tm_text()`: (複合) 点、(複合) 線、(複合) ポリゴンのテキスト - - - -レイヤ種別の一覧は、 `help("tmap-element")` とすることで表示される。 Figure \@ref(fig:tmshape) の右側のパネルでは、塗りつぶし (fill) レイヤの上に境界 (borders) を重ねた結果を示す。 - - - - ### 地図オブジェクト {#map-obj} @@ -163,8 +163,6 @@ map_nz1 = map_nz + tm_shape(nz_elev) + tm_raster(col_alpha = 0.7) ``` - - 先に作成した `map_nz` オブジェクトをベースに、新しいマップオブジェクト `map_nz1` を作成する。このオブジェクトには、ニュージーランド全土の平均標高を表す別の図形 (`nz_elev`) が含まれている (Figure \@ref(fig:tmlayers) 左図)。 さらに図形やレイヤを追加することもできる。以下のコードでは、ニュージーランドの[領海](https://en.wikipedia.org/wiki/Territorial_waters)を表す `nz_water` を作成し、作成した線を既存の地図オブジェクトに追加している。 @@ -237,18 +235,18 @@ Base R のプロットと同様に、美観を定義する引数もまた、様 ```{r 08-mapping-9, eval=FALSE} plot(st_geometry(nz), col = nz$Land_area) # 成功 -tm_shape(nz) + tm_fill(col = nz$Land_area) # 失敗 +tm_shape(nz) + tm_fill(fill = nz$Land_area) # 失敗 #> Error: palette should be a character value ``` -代わりに、`col` (および、ラインレイヤのための `lwd`、ポイントレイヤのための `size` など、異なることがある他の美観) は、プロットされるジオメトリに関連する属性を指定する文字列を必要とする。 -したがって、次のように望ましい結果を得ることができる (Figure \@ref(fig:tmcol) の右側のパネルにプロットされている)。 +代わりに、`fill` (および、ラインレイヤのための `lwd`、ポイントレイヤのための `size` など、異なることがある他の美観) は、プロットされるジオメトリに関連する属性を指定する文字列を必要とする。 +したがって、次のように望ましい結果を得ることができる (Figure \@ref(fig:tmcol) 右図)。 ```{r 08-mapping-10, fig.show='hide', message=FALSE} tm_shape(nz) + tm_fill(fill = "Land_area") ``` -```{r tmcol, message=FALSE, fig.cap="数値色フィールドの base (左) と tmap (右) の処理方法の比較。", fig.scap="Comparison of base graphics and tmap", echo=FALSE, out.width="45%", fig.show='hold', warning=FALSE, message=FALSE} +```{r tmcol, message=FALSE, fig.cap="数値色フィールドの base (左) と tmap (右) の処理方法の比較。", fig.scap="Comparison of base graphics and tmap", echo=FALSE, fig.show='hold', warning=FALSE, message=FALSE, fig.height=6, out.width="45%"} plot(nz["Land_area"]) tm_shape(nz) + tm_fill(fill = "Land_area") ``` @@ -256,14 +254,14 @@ tm_shape(nz) + tm_fill(fill = "Land_area") 視覚化の変数には、`.scale`、`.legend`、`.free` という文字列を後ろにつけた３つの追加引数がある。 例えば、`tm_fill()` には `fill`、`fill.scale`、`fill.legend`、`fill.free` といった引数がある。 `.scale` 引数は、地図と凡例での表示方法を指定し (Section \@ref(scales))、`.legend` はタイトル、方向、位置を指定する (Section \@ref(legends))。 -`.free` 引数は、多くのファセットをもつ地図で、ファセットによって縮尺や凡例が変わる場合などに使用する (Section \@ref(faceted-maps))。 +`.free` 引数は、多くのファセットをもつ地図で、ファセットによって縮尺や凡例が変わる場合などに使用する。 ### 縮尺 {#scales} \index{tmap (package)!scales} 縮尺は、地図と凡例で値がどのように表示されるかを制御し、視覚的変数に大きく依存する。 例えば、`col` という変数の場合、`col.scale` が空間オブジェクトの色が値にどう対応するかを制御する。`size` という変数の場合、`size.scale` は大きさが値にどう対応するかを制御する。 -デフォルトでは、`tm_scale()` が使われ、データ種別 (因子型、数値型、整数型) によって自動的に設定を選択する。 +デフォルトでは、`tm_scale()` が使われ、入力データ種別 (因子型、数値型、整数型) によって自動的に設定を選択する。 \index{tmap (package)!color breaks} ポリゴンの塗りつぶしの設定によって縮尺がどのように機能するか見てみよう。 @@ -328,39 +326,39 @@ The `tm_scale_intervals()` 関数は、入力データ値を間隔セットに `style` は **tmap** 関数の引数ではあるが、これはもともと `classInt::classIntervals()` の引数である。よって、詳細はこの関数のヘルプを参照。 ``` -```{r break-styles, message=FALSE, fig.cap="tmap の style 引数で設定するビン方法の違いによる説明。", , fig.scap="Illustration of different binning methods using tmap.", echo=FALSE} -source("https://github.com/geocompx/geocompr/raw/main/code/09-break-styles.R", print.eval = TRUE) +```{r break-styles, message=FALSE, fig.cap="tmap の style 引数で設定するビン方法の違いによる説明。", , fig.scap="Illustration of different binning methods using tmap.", echo=FALSE, warning=FALSE, fig.width=8} +source("code/09-break-styles.R", print.eval = TRUE) ``` \index{tmap (package)!continuous scale} -`tm_scale_continuous()` 関数は、連続色フィールドの多くの色を提示する。連続ラスタによく用いられる。 +`tm_scale_continuous()` 関数は、連続色フィールドの色を提示する。連続ラスタによく用いられる (Figure \@ref(fig:concat) 左図)。 分布が偏っている場合に対して、`tm_scale_continuous_log()` や `tm_scale_continuous_log1p()` といった派生もある。 \index{tmap (package)!categorical scale} -最後に、`tm_scale_categorical()` は、カテゴリ値を代表し、各カテゴリに固有の色が割り当てられる (Figure \@ref(fig:concat))。 +最後に、`tm_scale_categorical()` は、カテゴリ値を代表し、各カテゴリに固有の色が割り当てられる (Figure \@ref(fig:concat) 右図)。 -```{r concat, message=FALSE, fig.cap="tmap における連続スケールとカテゴリスケールの例", echo=FALSE} +```{r concat, message=FALSE, fig.cap="tmap における連続スケールとカテゴリスケールの例", echo=FALSE, fig.width=8} library(tmap) library(spData) library(spDataLarge) m_cont1 = tm_shape(nz) + - tm_polygons(fill = "Median_income", fill.scale = tm_scale_continuous()) + + tm_polygons(fill = "Median_income", fill.scale = tm_scale_continuous(n = 5)) + tm_title('tm_scale_continuous()', fontfamily = "monospace") + - tm_layout(legend.position = tm_pos_auto_in()) + tm_layout(legend.position = tm_pos_auto_in(), scale = 0.9) m_cat1 = tm_shape(nz) + tm_polygons(fill = "Island", fill.scale = tm_scale_categorical()) + tm_title('tm_scale_categorical()', fontfamily = "monospace") + - tm_layout(legend.position = tm_pos_auto_in()) + tm_layout(legend.position = tm_pos_auto_in(), scale = 0.9) tmap_arrange(m_cont1, m_cat1) ``` \index{tmap (package)!palettes} パレットは、ビンに関連付けられ、前述の `breaks`、`n`、`style` 引数で決定される色域を定義する。 -デフォルトの色パレットは `tm_layout()` で指定されている (詳しくは Section \@ref(layouts) を参照)。しかし、引数 `palette` を使えば、すぐに変更することができる。 この引数には色ベクトルまたは新しい色パレット名を与えるが、`tmaptools::palette_explorer()` で対話的に選択することができる。 プレフィックスとして `-` を付けると、パレットの順序を逆にすることができる。 ```{block2 visual-vars-values, type='rmdnote'} -入力データに応じた色パレットなどの美観に関わるデフォルト`値`は、`tmap_options()$values.var` で確認することができる。 +入力データに応じた色パレットなどの美観に関わるデフォルト`値`は、`tmap_options() で確認することができる。 +例として、`tmap_options()$values.var` を実行してみよう。 ``` 色パレット\index{map making!color palettes}は大きく分けて三種類ある。カテゴリ、連続、発散である (Figure \@ref(fig:colpal))。目的に応じてこの三種類を使い分ける。^[ @@ -369,12 +367,7 @@ tmap_arrange(m_cont1, m_cat1) ] カテゴリパレットは、区別しやすい色で構成されており、州名や土地被覆クラスなど、特定の順序を持たないカテゴリデータに最適である。 色は直感的にわかるように、例えば川は青、牧草地は緑にする。 -カテゴリを増やしすぎると、大きな凡例や多くの色を使った地図は理解できないことがある。 - - - - - +カテゴリを増やしすぎると、大きな凡例や多くの色を使った地図は理解できないことがある。^[国別に色分けするように、個別のポリゴンが多い場合、`fill = "MAP_COLORS"` としてみよう。隣接するポリゴンにユニークな色を設定する。] 次に、連続色パレットをみてみよう。 連続色パレットは、例えば明るい色から暗い色への勾配に沿っており (明るい色は低い値を表す傾向がある)、連続的な (数値) 変数に適している。 @@ -425,10 +418,6 @@ many_palette_plotter(c(all_default_pals$div, all_default_pals$seq, all_default_p 例えば、緑は植物や低地を表し、青は水や涼しさを連想させる色である。 また、情報を効果的に伝えるために、カラーパレットは分かりやすいものが望ましい。 どの数値が低く、どの数値が高いかが明確で、色も徐々に変化することが望ましい。 - - - 第二に、色の変化は、多くの人がアクセスできるものでなければならない。 そのため、色弱者用のパレットをできるだけ多く使うことが大切である。^[`cols4all::c4a_gui()` の "Color Blind Friendliness" パネルの "Color vision" オプションを参照。] @@ -442,8 +431,9 @@ many_palette_plotter(c(all_default_pals$div, all_default_pals$seq, all_default_p 以下のコードは、変数名 `Land_area` よりも魅力的な名前を提供することで、この機能を示している (`expression()` は上付き文字を設定するために使用)。 ```{r 08-mapping-11} +#| eval: false legend_title = expression("Area (km"^2*")") -map_nza = tm_shape(nz) + +tm_shape(nz) + tm_polygons(fill = "Land_area", fill.legend = tm_legend(title = legend_title)) ``` @@ -451,10 +441,26 @@ map_nza = tm_shape(nz) + 凡例の位置は、`position` 引数で設定する。 ```{r} +#| eval: false +tm_shape(nz) + + tm_polygons(fill = "Land_area", + fill.legend = tm_legend(title = legend_title, + orientation = "landscape", + position = tm_pos_out("center", "bottom"))) +``` + +```{r} +#| eval: false +#| echo: false +legend_title = expression("Area (km"^2*")") +map_nza = tm_shape(nz) + + tm_polygons(fill = "Land_area", fill.legend = tm_legend(title = legend_title), position = tm_pos_out("right", "top")) map_nza2 = tm_shape(nz) + - tm_polygons(fill = "Land_area", fill.legend = tm_legend(title = legend_title, - orientation = "landscape", - position = tm_pos_out("center", "bottom"))) + tm_polygons(fill = "Land_area", + fill.legend = tm_legend(title = legend_title, + orientation = "landscape", + position = tm_pos_out("center", "bottom"))) +tmap_arrange(map_nza, map_nza2) ``` 凡例の位置 (その他の **tmap** の地図要素の位置も) は、いくつかの関数でカスタマイズできる。 @@ -465,13 +471,13 @@ map_nza2 = tm_shape(nz) + - `tm_pos_in()`: は、地図フレーム内に配置する。 最初の引数は `"left"`、`"center"`、`"right"` のいずれかで、２番目の引数は `"bottom"`、`"center"`、`"top"` のいずれかである。 -または、２つの値のベクトル (または 0 から 1 の２つの値) を与えても良い。この場合、凡例は地図フレーム領域の外に配置される。 +または、2つの値のベクトル (または 0 から 1 の 2 つの値) を与えても良い。この場合、凡例は地図フレーム領域の外に配置される。 ### レイアウト {#layouts} \index{tmap (package)!layouts} 地図レイアウトとは、すべての地図要素を組み合わせて、まとまりのある地図にすることである。 -マップ要素には、マップされるオブジェクト、タイトル、縮尺バー、地図グリッド、マージンなどがあり、前のセクションで説明したカラー設定は、マップの見え方に影響を与えるパレットとブレークポイントに関連している。 +地図要素には、マップされるオブジェクト、地図グリッド (メッシュ)、縮尺バー、タイトル、マージンなどがあり、前のセクションで説明したカラー設定は、地図の見え方に影響を与えるパレットとブレークポイントに関連している。 どちらも微妙な変化をもたらすだろうが、地図が残す印象には同じように大きな影響を与える。 経緯線網\index{tmap (package)!graticules}、北矢印\index{tmap (package)!north arrows}、スケールバー\index{tmap (package)!scale bars}、タイトルなどの追加要素には、それぞれ関数がある。`tm_graticules()`、`tm_compass()`、`tm_scalebar()`、`tm_title()` である (Figure \@ref(fig:na-sb))。^[この他、`tm_grid()`、`tm_logo()`、`tm_credits()` がある。]訳注：mac で title に日本語を使用する場合は、fontfamily = "HiraginoSans-W3" を追加する。 @@ -499,14 +505,13 @@ source("code/09-layout1.R", print.eval = TRUE) `tm_layout()` の他の引数は、マップが配置されるキャンバスとの関係で、マップの多くの側面を制御する。 ここでは、便利なレイアウト設定をご紹介する (一部、Figure \@ref(fig:layout2))。 -- `outer.margin` を含むマージン設定と `inner.margin` +- `inner.margin` と `outer.margin` のマージン設定 - `fontface` で制御されるフォント設定と `fontfamily` - 凡例設定には、`legend.show` (凡例を表示すかどうか)、`legend.only` (地図を省略するか)、`legend.outside` (凡例を地図の外に出すか) などの二値オプションや、次のような多肢選択設定がある。 `legend.position` - フレーム幅 (`frame.lwd`) と二重線 (`frame.double.line`) を許可するオプション - `sepia.intensity` (マップのセピア度合) と `saturation` (色・グレースケール) を制御する色設定 -```{r layout2, message=FALSE, fig.cap="選択したレイアウトオプションの説明図。", echo=FALSE, fig.asp=0.56, warning=FALSE} -# todo: add more useful settings to this plot +```{r layout2, message=FALSE, fig.cap="選択したレイアウトオプションの説明図。", echo=FALSE, fig.asp=0.56, warning=FALSE, fig.width=8} source("code/09-layout2.R", print.eval = TRUE) ``` @@ -623,7 +628,7 @@ ins_vp = viewport(width = ins_dim[1] * 0.5, height = ins_dim[2] * 0.5, 最後に、新規にキャンバスを作り、主地図を表示し、主地図の領域内に差し込み地図を配置する。 -```{r insetmap1, message=FALSE, fig.cap="差し込み地図で背景を説明 - ニュージーランドの南アルプスの中央部の位置。", fig.scap="Inset map providing a context."} +```{r insetmap1, message=FALSE, fig.cap="差し込み地図で背景を説明 - ニュージーランドの南アルプスの中央部の位置。", fig.scap="Inset map providing a context.", fig.width=9} grid.newpage() print(nz_height_map, vp = main_vp) pushViewport(main_vp) @@ -635,10 +640,11 @@ print(nz_map, vp = ins_vp) また、差し込み地図は、非連続なエリアを1つの地図にするために使用する。 おそらく最もよく使われる例は、アメリカ合衆国の地図で、アメリカ本土とハワイ、アラスカで構成されている。 このようなケースでは、個々の差し込みに最適なプロジェクションを見つけることが非常に重要である (詳しくは Chapter \@ref(reproj-geo-data) を参照)。 -`tm_shape()` の引数 `projection` に US National Atlas Equal Area の EPSG コードを指定すれば、米国本土の地図に US National Atlas Equal Area を使用することができる。 +`tm_shape()` の引数 `crs` に US National Atlas Equal Area の EPSG コードを指定すれば、米国本土の地図に US National Atlas Equal Area を使用することができる。 ```{r 08-mapping-19} -us_states_map = tm_shape(us_states, projection = "EPSG:2163") + tm_polygons() + +us_states_map = tm_shape(us_states, crs = "EPSG:9311") + + tm_polygons() + tm_layout(frame = FALSE) ``` @@ -656,7 +662,7 @@ alaska_map = tm_shape(alaska) + tm_layout(frame = FALSE, bg.color = NA) ``` -これら3つの地図を組み合わせて配置することで、最終的な地図ができあがる。 +これら 3 つの地図を組み合わせ、サイズを調整し配置することで、最終的な地図ができあがる。 ```{r insetmap2, message=FALSE, fig.cap="アメリカ合衆国の地図。", warning=FALSE} us_states_map @@ -664,7 +670,7 @@ print(hawaii_map, vp = grid::viewport(0.35, 0.1, width = 0.2, height = 0.1)) print(alaska_map, vp = grid::viewport(0.15, 0.15, width = 0.3, height = 0.3)) ``` -上記で紹介したコードはコンパクトで、他の差し込み地図のベースとして使用することができるが、その結果、Figure \@ref(fig:insetmap2)、ハワイとアラスカの位置がうまく表現されていないことがわかる。 +上記で紹介したコードはコンパクトで、他の差し込み地図のベースとして使用することができるが、その結果、Figure \@ref(fig:insetmap2)、ハワイとアラスカの位置とサイズがうまく表現されていないことがわかる。 より詳細なアプローチについては、**geocompkg** の [ `us-map` ](https://geocompx.github.io/geocompkg/articles/us-map.html) vignette を参照。 ## アニメーションマップ {#animated-maps} @@ -684,7 +690,7 @@ R でアニメーションを生成する方法はいくつかあり、**ggplot2 Figure \@ref(fig:urban-animated) は、アニメーションマップの簡単な例である。 ファセットプロットとは異なり、複数の地図を一画面に押し込むことはなく、世界で最も人口の多い集積地の空間分布が時間とともにどのように進化していくかを見ることができる (アニメーション版は同書のウェブサイトを参照)。 -```{r urban-animated, message=FALSE, fig.cap="1950年から2030年までの、国連による人口予測に基づく都市集積の上位30位を示したアニメーション地図。アニメーション版は、geocompr.robinlovelace.net で見ることができる。", fig.scap="Animated map showing the top 30 largest 'urban agglomerations'.", echo=FALSE} +```{r urban-animated, message=FALSE, fig.cap="1950年から2030年までの、国連による人口予測に基づく都市集積の上位30位を示したアニメーション地図。アニメーション版は、geocompr.robinlovelace.net で見ることができる。", fig.scap="Animated map showing the top 30 largest 'urban agglomerations'.", echo=FALSE, fig.height=3.3} if (knitr::is_latex_output()){ knitr::include_graphics("figures/urban-animated.png") } else if (knitr::is_html_output()){ @@ -697,7 +703,7 @@ source("https://github.com/Robinlovelace/geocompr/raw/main/code/09-urban-animati ``` Figure \@ref(fig:urban-animated) に示したアニメーションマップは、Section \@ref(faceted-maps) で示したファセット・マップを生成するのと同じ **tmap** 技術を使用して作成することができる。 -ただし、`tm_facets()` の引数に関連して、2つの違いがある。 +ただし、`tm_facets_wrap()` の引数に関連して、2 つの違いがある。 - `by = "year"` の代わりに `along = "year"` が使用される。 - `free.coords = FALSE`、これは各マップ反復のためのマップ範囲を維持する。 @@ -707,7 +713,7 @@ Figure \@ref(fig:urban-animated) に示したアニメーションマップは ```{r 08-mapping-22} urb_anim = tm_shape(world) + tm_polygons() + tm_shape(urban_agglomerations) + tm_symbols(size = "population_millions") + - tm_facets(by = "year", nrow = 1, ncol = 1, free.coords = FALSE) + tm_facets_wrap(by = "year", nrow = 1, ncol = 1, free.coords = FALSE) ``` その結果、`urb_anim`、各年度の個別の地図のセットを表している。 @@ -720,7 +726,7 @@ tmap_animation(urb_anim, filename = "urb_anim.gif", delay = 25) アニメーション・マップの威力を示すもう一つの例が、Figure \@ref(fig:animus) にある。 これは、アメリカにおける州の発達を示すもので、最初は東部で形成され、その後徐々に西部へ、最後は内陸部へと発展していった。 -このマップを再現するためのコードは、スクリプト `09-usboundaries.R` に記載されている。 +このマップを再現するためのコードは、本書の GitHub レポジトリのスクリプト `code/09-usboundaries.R` に記載されている。 ```{r 08-mapping-24, echo=FALSE, eval=FALSE} source("https://github.com/Robinlovelace/geocompr/raw/main/code/09-usboundaries.R") @@ -746,7 +752,7 @@ if (knitr::is_latex_output()){ より高度なインタラクティブ機能としては、下記の **mapdeck** の例で示したように、地図を傾けたり回転させたりする機能や、ユーザーがパンやズームをすると自動的に更新される「動的にリンクした」サブプロット [@pezanowski_senseplace3_2018] を提供する機能などが挙げられる。 しかし、インタラクティブ性の最も重要なタイプは、インタラクティブまたは「スリッピー」なウェブマップ上での地理データの表示である。 -2015年にリリースされた **leaflet** (leaflet JavaScript ライブラリを使用) パッケージは、R 内からインタラクティブな Web マップの作成に革命をもたらし、多くのパッケージがこれらの基盤の上に新機能を追加し (例：**leaflet.extras**)、Web マップの作成を静的マップ作成と同じくらいシンプルにしている (例：**mapview** や **tmap** な)）。 +2015年にリリースされた **leaflet** (leaflet JavaScript ライブラリを使用) パッケージは、R 内からインタラクティブな Web マップの作成に革命をもたらし、多くのパッケージがこれらの基盤の上に新機能を追加し (例：**leaflet.extras2**)、Web マップの作成を静的マップ作成と同じくらいシンプルにしている (例：**mapview** や **tmap** など)）。 ここでは、各アプローチを逆の順番で説明する。 ここでは、**tmap** (すでに学習済みの構文)、**mapview**\index{mapview (package)}、**mapdeck**\index{mapdeck (package)}そして最後に **leaflet** \index{leaflet (package)} (対話型マップの低レベル制御を提供) を使って、動く地図を作成する方法を探究する。 @@ -769,11 +775,8 @@ if (knitr::is_latex_output()){ # file.copy("tmview-1.html", "~/geocompr/geocompr.github.io/static/img/tmview-1.html") knitr::include_url("https://geocompx.org/static/img/tmview-1.html") } - ``` - - インタラクティブモードが「オン」になったので、**tmap** で作成したすべての地図が起動する (インタラクティブな地図を作成する別の方法として、`tmap_leaflet` 機能がある)。 このインタラクティブモードの特筆すべき点は、以下のデモのように `tm_basemap()` (または `tmap_options()`) でベースマップを指定できることである (結果は表示していない)。 @@ -782,7 +785,7 @@ map_nz + tm_basemap(server = "OpenTopoMap") ``` あまり知られていないが、**tmap** の表示モードは、ファセット・プロットにも対応している。 -この場合、`tm_facets()` の引数 `sync` を使用すると、以下のコードで作成した Figure \@ref(fig:sync) のように、ズームとパンの設定が同期した複数のマップを作成することができる。 +この場合、`tm_facets_wrap()` の引数 `sync` を使用すると、以下のコードで作成した Figure \@ref(fig:sync) のように、ズームとパンの設定が同期した複数のマップを作成することができる。 ```{r 08-mapping-27, eval=FALSE} world_coffee = left_join(world, coffee_data, by = "name_long") @@ -914,7 +917,7 @@ leaflet(data = cycle_hire) |> ``` ```{r leaflet, message=FALSE, fig.cap="ロンドン市内の自転車レンタルポイントを紹介したリーフレットパッケージの実例。インタラクティブ版は[オンライン](https://geocompr.github.io/img/leaflet.html)を参照。", fig.scap="The leaflet package in action.", echo=FALSE} -if (knitr::is_latex_output() | knitr::is_html_output()){ +if (knitr::is_latex_output() || knitr::is_html_output()){ knitr::include_graphics("figures/leaflet-1.png") } else { # pre-generated for https://github.com/ropensci/stplanr/issues/385 @@ -952,15 +955,15 @@ Section \@ref(interactive-maps) で示されたインタラクティブなウェ R の場合は幸運なことに、**shiny**\index{shiny (package)} を使って、ウェブ地図アプリケーションを素早く作成できるようになった。 オープンソース本 [Mastering Shiny](https://mastering-shiny.org/) で説明されているように、 **shiny** は、R コードをインタラクティブなウェブアプリに変換する R パッケージでありフレームワークである [@wickham_mastering_2021]。 -`tmap::renderTmap()` と [`leaflet::renderLeaflet()`](https://rstudio.github.io/leaflet/shiny.html) を使うことで、shiny アプリにインタラクティブ地図を追加することができる。 -このセクションでは、ウェブ地図の観点から **shiny** の基本を学び、100行未満のコードで全画面マッピングアプリケーションを完成させることができる。 + [`leaflet::renderLeaflet()`](https://rstudio.github.io/leaflet/shiny.html) を使うことで、shiny アプリにインタラクティブ地図を追加することができる。 +このセクションでは、ウェブ地図の観点から **shiny** の基本を学び、100 行未満のコードで全画面マッピングアプリケーションを完成させることができる。 -**shiny** の仕組みは、[shiny.posit.co](https://shiny.posit.co/) に詳しく書かれている。これは「フロントエンド」 (ユーザーが見る部分) と「バックエンド」コードという２つの構成要素を強調している。 +**shiny** の仕組みは、[shiny.posit.co](https://shiny.posit.co/) に詳しく書かれている。これは「フロントエンド」 (ユーザーが見る部分) と「バックエンド」コードという 2 つの構成要素を強調している。 **shiny** アプリでは、これらの要素は通常、`app フォルダ`内にある `app.R` という R スクリプト内の `ui` と `server` というオブジェクトで作成される。 -これにより、ウェブマッピングのアプリケーションを1つのファイルで表現することができ、例えば、本書の GitHub リポジトリにある [`CycleHireApp/app.R`](https://github.com/Robinlovelace/geocompr/blob/main/apps/CycleHireApp/app.R) ファイルのような単一のファイルで表現できる。 +これにより、ウェブマッピングのアプリケーションを 1 つのファイルで表現することができ、例えば、本書の GitHub リポジトリにある [`CycleHireApp/app.R`](https://github.com/Robinlovelace/geocompr/blob/main/apps/CycleHireApp/app.R) ファイルのような単一のファイルで表現できる。 ```{block2 shiny, type = 'rmdnote'} -**shiny** アプリでは、これらはしばしば `ui.R` (ユーザーインターフェースの略) と `server.R` ファイルに分けられる。この命名規則は、一般向けの Web サイトで shiny アプリを提供するサーバーサイド Linux アプリケーション、 `shiny-server` で使用されている。 +**shiny** アプリでは、これらは `ui.R` (ユーザーインターフェースの略) と `server.R` ファイルに分けられることが多い。この命名規則は、一般向けの Web サイトで shiny アプリを提供するサーバーサイド Linux アプリケーション、 `shiny-server` で使用されている。 `shiny-server` は、'app フォルダ' 内にある `app.R` という単一ファイルで定義されるアプリを提供することもある。 詳細は https://github.com/rstudio/shiny-server 。 ``` @@ -1050,9 +1053,9 @@ knitr::include_graphics("figures/09_cycle_hire_app.png") 最も成熟した選択肢は、コアな空間パッケージである **sf** と **terra** が提供する `plot()` メソッドを使用することである。 これらのセクションで触れていないのは、ベクタとラスタのオブジェクトのプロットメソッドは、結果が同じプロットエリアに描画される場合、組み合わせることができるということである (**sf** プロットのキーやマルチバンドのラスタなどの要素はこれを邪魔する)。 この動作は、Figure \@ref(fig:nz-plot) を生成する次のコードチャンクで説明される。 -`plot()` には他にも多くのオプションがあり、`?plot` のヘルプページと **sf** vignette [`sf5`](https://cran.r-project.org/package=sf/vignettes/sf5.html) のリンクをたどって調べることができる (訳注：[vignette 日本語版](https://www.uclmail.net/users/babayoshihiko/R/))。 +`plot()` には他にも多くのオプションがあり、`?plot` のヘルプページと **sf** 5番目の vignette [`sf5`](https://cran.r-project.org/package=sf/vignettes/sf5.html) のリンクをたどって調べることができる (訳注：[vignette 日本語版](https://www.uclmail.net/users/babayoshihiko/R/))。 -```{r nz-plot, message=FALSE, fig.cap="plot() で作成したニュージーランドの地図。右の凡例は標高 (海抜 1000 m) を示している。", fig.scap="Map of New Zealand created with plot()."} +```{r nz-plot, message=FALSE, fig.cap="plot() で作成したニュージーランドの地図。右の凡例は標高 (海抜 1000 m) を示している。", fig.scap="Map of New Zealand created with plot().", fig.height=3.5} g = st_graticule(nz, lon = c(170, 175), lat = c(-45, -40, -35)) plot(nz_water, graticule = g, axes = TRUE, col = "blue") terra::plot(nz_elev / 1000, add = TRUE, axes = FALSE) @@ -1097,7 +1100,7 @@ ggplot() + scale_fill_continuous(na.value = NA) ``` -```{r nz-gg2, message=FALSE, fig.cap="ggplot2 のみ (左) と ggplot2 と ggspatial (右) で生成したニュージーランド地図の比較。", out.width="45%", fig.show='hold', echo=FALSE} +```{r nz-gg2, message=FALSE, fig.cap="ggplot2 のみ (左) と ggplot2 と ggspatial (右) で生成したニュージーランド地図の比較。", out.width="45%", fig.show='hold', echo=FALSE, warning=FALSE, fig.height=7} library(ggplot2) g1 = ggplot() + geom_sf(data = nz, aes(fill = Median_income)) + geom_sf(data = nz_height) + @@ -1135,8 +1138,6 @@ Table \@ref(tab:map-gpkg) は、さまざまな地図製作パッケージが利 特に注目すべきは **mapsf** で、コロプレス図、プロポーショナルシンボル地図、フロー地図など、さまざまな地理的視覚化を生成することができる。 これらは、[`mapsf`](https://cran.r-project.org/package=mapsf/vignettes/mapsf.html)\index{mapsf (package)} vignette に記載されている。 - - Table \@ref(tab:map-spkg) に示すように、いくつかのパッケージは、特定のマップタイプに焦点を当てている。 地理空間を歪めたカルトグラムの作成、ラインマップの作成、ポリゴンの正六角形グリッドへの変換、複雑なデータを地理的トポロジーを表すグリッド上に可視化し、３次元表現をするパッケージである。 @@ -1151,9 +1152,10 @@ knitr::kable(map_spkg_df, ``` + しかし、前述のパッケージはいずれも、データ準備や地図作成のアプローチが異なっている。 -次の段落では、**cartogram** パッケージ \index{cartogram (package)} にのみ焦点を当てる。 -そのため、[linemap](https://github.com/rCarto/linemap)\index{linemap (package)}、[geogrid](https://github.com/jbaileyh/geogrid)\index{geogrid (package)}、 [geofacet](https://github.com/hafen/geofacet) \index{geofacet (package)}、[rayshader](https://github.com/tylermorganwall/rayshader)\index{rayshader (package)} のドキュメントを読んで、より詳しく知ることを勧める。 +次の段落では、**cartogram** パッケージ [@R-cartogram]\index{cartogram (package)} にのみ焦点を当てる。 +そのため、[geogrid](https://github.com/jbaileyh/geogrid)\index{geogrid (package)}、 [geofacet](https://github.com/hafen/geofacet)\index{geofacet (package)}、[linemap](https://github.com/rCarto/linemap)\index{linemap (package)}、[tanaka](https://github.com/riatelab/tanaka)\index{tanaka (package)}、[rayshader](https://github.com/tylermorganwall/rayshader)\index{rayshader (package)} のドキュメントを読んで、より詳しく知ることを勧める。 カルトグラムとは、地図変数を表現するために一定の幾何学的な歪みを持たせた地図のことである (訳注：訳語は[東北大学空間計画科学研究室](http://www.plan.civil.tohoku.ac.jp/inoue/research/cartogram/)に準じた)。 このような地図の作成は、R では **cartogram** を用いることで、連続・非連続の面積カルトグラム (area cartogram) を作成することが可能である。 @@ -1162,7 +1164,7 @@ knitr::kable(map_spkg_df, `cartogram_cont()` 機能は、連続した面積のカルトグラムを作成する。 入力として、`sf` オブジェクトと変数名 (列) を受け取る。 さらに、`intermax` 引数 (カルトグラム変換の最大反復回数) を変更することが可能である。 -例えば、ニュージーランドの地域の所得の中央値を連続カルトグラム (Figure \@ref(fig:cartomap1) の右側のパネル) で表すと、次のようになる。 +例えば、ニュージーランドの地域の所得の中央値を連続カルトグラム (Figure \@ref(fig:cartomap1) 右図) で表すと、次のようになる。 ```{r 08-mapping-39, fig.show='hide', message=FALSE} library(cartogram) @@ -1170,15 +1172,16 @@ nz_carto = cartogram_cont(nz, "Median_income", itermax = 5) tm_shape(nz_carto) + tm_polygons("Median_income") ``` -```{r cartomap1, echo=FALSE, message=FALSE, fig.cap="標準地図 (左) と連続範囲 (右) の比較。", fig.scap="Comparison of standard map and continuous area cartogram.", warning=FALSE} -carto_map1 = tm_shape(nz) + +```{r cartomap1, echo=FALSE, message=FALSE, fig.cap="標準地図 (左) と連続範囲 (右) の比較。", fig.scap="Comparison of standard map and continuous area cartogram.", warning=FALSE, fig.width=9} +carto_map1 = tm_shape(nz) + tm_polygons("Median_income", fill.scale = tm_scale(values = "Greens"), - fill.legend = tm_legend(title = "Median income (NZD)")) -carto_map2 = tm_shape(nz_carto) + + fill.legend = tm_legend_hide()) +carto_map2 = tm_shape(nz_carto) + tm_polygons("Median_income", fill.scale = tm_scale(values = "Greens"), - fill.legend = tm_legend(title = "Median income (NZD)")) + fill.legend = tm_legend(title = "Median income (NZD)", + position = c("right", "bottom"))) tmap_arrange(carto_map1, carto_map2) ``` @@ -1188,41 +1191,30 @@ tmap_arrange(carto_map1, carto_map2) 以下のコードは、米国各州の人口の非連続面積とサークルエリア (Dorling) カルトグラムの作成例である (Figure \@ref(fig:cartomap2))。 ```{r 08-mapping-40, fig.show='hide', message=FALSE} -us_states2163 = st_transform(us_states, "EPSG:2163") -us_states2163_ncont = cartogram_ncont(us_states2163, "total_pop_15") -us_states2163_dorling = cartogram_dorling(us_states2163, "total_pop_15") +us_states9311 = st_transform(us_states, "EPSG:9311") +us_states9311_ncont = cartogram_ncont(us_states9311, "total_pop_15") +us_states9311_dorling = cartogram_dorling(us_states9311, "total_pop_15") ``` - - - - - - - - - - - - - - - - - ## 演習 @@ -1231,6 +1223,3 @@ res = knitr::knit_child('_09-ex.Rmd', quiet = TRUE, options = list(include = FALSE, eval = FALSE)) cat(res, sep = '\n') ``` - - - diff --git a/10-gis-ja.Rmd b/10-gis-ja.Rmd index 0d5a4d1..22b3b5d 100644 --- a/10-gis-ja.Rmd +++ b/10-gis-ja.Rmd @@ -1,8 +1,12 @@ # GISソフトウェアへのブリッジ {#gis} +```{r, include=FALSE} +source("code/before_script.R") +``` + ## 必須パッケージ {- #prerequisites-10} -- 本章では、QGIS\index{QGIS}、SAGA\index{SAGA}、GRASS\index{GRASS} がインストールされていること、および以下のパッケージが添付されていることを条件とする。 +- 本章では、QGIS\index{QGIS}、SAGA\index{SAGA}、GRASS GIS\index{GRASS GIS} がインストールされていること、および以下のパッケージが添付されていることを条件とする。 ```{r 09-gis-1, message=FALSE, echo=FALSE} library(sf) @@ -38,16 +42,16 @@ has_qgis_plugins = RStudio や VS Code などの対話型開発環境を使った作業では、通常ソースエディタでソースファイルに書き込み、`Ctrl+Enter` などのショートカットキーでコードの対話的実行を制御する。 CLI は R だけのものではない。初期のコンピュータ環境は、おおむねコマンドライン「シェル」に依存しており、GUI\index{graphical user interface} が一般的になったのは、コンピューターのマウスが普及した1990年代以降のことである。 -例えば、最も長い歴史を持つ GIS\index{GIS} プログラムの一つである GRASS は、洗練された GUI [@landa_new_2008] を獲得するまでは、主にコマンドラインでの対話に依存していた。 +例えば、最も長い歴史を持つ GIS\index{GIS} プログラムの一つである GRASS GIS は、洗練された GUI [@landa_new_2008] を獲得するまでは、主にコマンドラインでの対話に依存していた。 ほとんどの GIS\index{GIS} パッケージは、グラフィカルユーザーインターフェイス (Grafical User Interface, GUI\index{graphical user interface} を採用している。 -QGIS\index{QGIS}、SAGA\index{SAGA}、GRASS\index{GRASS}、gvSIG を、システムターミナルや組み込み CLI\index{command-line interface} から操作することも可能であるが、「マウス操作」が一般的である。 +QGIS\index{QGIS}、SAGA\index{SAGA}、GRASS GIS\index{GRASS GIS}、gvSIG を、システムターミナルや組み込み CLI\index{command-line interface} から操作することも可能であるが、「マウス操作」が一般的である。 これは多くの GIS\index{GIS} ユーザーがコマンドラインの利点を見逃していることを意味する。 QGIS\index{QGIS} [@sherman_desktop_2008] の作成者によれば、 > 「近代的」GIS ソフトウェアの発展に伴い、マウス操作を好む人がほとんどである。それはいいことであるが、コマンドラインには、とてつもない量の柔軟性とパワーが待っている。繰り返しコマンドラインで作業をする場合、同じことを GUI で行う場合よりも短い時間で行うことができる。 「CLI vs GUI」\index{graphical user interface} の議論は、敵対的になる必要はない。作業の内容 (フィーチャを描く時は GUI の方が優れている)、再現可能性、ユーザーのスキルセットに応じて、どちらの選択肢も利点はある。 -GRASS は、CLI をベースとしながらも GUI がある GIS の良い例である。 +GRASS GIS は、CLI をベースとしながらも GUI がある GIS の良い例である。 同様に、R は CLI\index{command-line interface} であり、そして RStudio\index{RStudio} のような IDE\index{IDE} が GUI を提供することで、アクセシビリティを向上している。 このように、ソフトウエアは CLI または GUI というように完全に分けられるものではない。 しかしながら、CLI は以下のような重要な点があることは押さえておきたい。 @@ -61,8 +65,7 @@ GRASS は、CLI をベースとしながらも GUI がある GIS の良い例で 一方、GUI ベースの GIS\index{GIS} システムにも有利な点がある。 - 学習曲線が「浅い」ので、新しい言語を何時間も学ぶことなく、地理データを探索し、可視化することができる -- トレース、スナップ、トポロジーツールなど、「デジタイジング」 (新しいベクタデータセットの作成) のための優れたサポートを提供する -(**mapedit** パッケージは、R から開いたブラウザウィンドウでフィーチャを素早く編集することができるが、専門的で大規模な地図製作のデジタイズはできない。)。 +- トレース、スナップ、トポロジーツールなど、「デジタイジング」 (新しいベクタデータセットの作成) のための優れたサポートを提供する^[**mapedit** パッケージは、R から開いたブラウザウィンドウでフィーチャを素早く編集することができるが、専門的で大規模な地図製作のデジタイズはできない。] - 地上基準点によるジオリファレンス (georeference、ラスタ画像と既存地図とのマッチング)、オルソレクティフィケーション (orthorectification) が可能である - 立体視マッピング (LiDAR、Structure from Motion など) に対応する @@ -83,19 +86,19 @@ R は GIS として設計されたものではない。しかし、専用の GIS GIS ブリッジを通すことで、R はさまざまな作業を実行でき、さらには CLI の持っている再現可能性、拡張性、生産性といった付加価値をもたらす。 さらに、R は、インタラクティブ/アニメーションマップ作成 (Chapter \@ref(adv-map) 参照) や空間統計モデリング ( Chapter \@ref(spatial-cv) 参照) など、ジオコンピューティング\index{geocomputation}のいくつかの分野では GIS を凌駕する性能を有している。 -この章では、3つの成熟したオープンソース GIS 製品への「ブリッジ」(Table \@ref(tab:gis-comp) 参照) に焦点を当てる。 +この章では、Table \@ref(tab:gis-comp) にまとめた 3 つの成熟したオープンソース GIS 製品への「ブリッジ」に焦点を当てる。 - QGIS\index{QGIS} (**qgisprocess**\index{qgisprocess (package)}; Section \@ref(rqgis)) - SAGA\index{SAGA} (**Rsagacmd**\index{Rsagacmd (package)}; Section \@ref(saga)) -- GRASS\index{GRASS} (**rgrass**\index{rgrass (package)}; Section \@ref(grass)) +- GRASS GIS\index{GRASS GIS} (**rgrass**\index{rgrass (package)}; Section \@ref(grass)) -ここでは取り上げないパッケージもある。例えば、R パッケージの **RPyGeo** は、R から商用 GIS へアクセスすることができるが、もうメンテナンスされていない。 -また、QGIS\index{QGIS} ([docs.qgis.org](https://docs.qgis.org/3.28/en/docs/training_manual/processing/r_intro.html) 参照) や GRASS GIS\index{GRASS} ([grasswiki.osgeo.org](https://grasswiki.osgeo.org/wiki/R_statistics/rgrass#R_within_GRASS) 参照) では、GIS から R を実行する環境も開発されている。 + +また、QGIS\index{QGIS} ([docs.qgis.org](https://docs.qgis.org/3.28/en/docs/training_manual/processing/r_intro.html) 参照) や GRASS GIS\index{GRASS GIS} ([grasswiki.osgeo.org](https://grasswiki.osgeo.org/wiki/R_statistics/rgrass#R_within_GRASS) 参照) など GIS から R を実行する環境も開発されている。 ```{r gis-comp, echo=FALSE, message=FALSE} library(dplyr) d = tibble( - "GIS" = c("QGIS", "SAGA", "GRASS"), + "GIS" = c("QGIS", "SAGA", "GRASS GIS"), "First release" = c("2002", "2004", "1982"), "No. functions" = c(">1000", ">600", ">500"), "Support" = c("hybrid", "hybrid", "hybrid") @@ -110,26 +113,28 @@ knitr::kable( caption.short = "Comparison between three open-source GIS.", booktabs = TRUE ) #|> -# kableExtra::add_footnote(label = "Comparing downloads of different providers is rather difficult (see http://spatialgalaxy.net/2011/12/19/qgis-users-around-the-world), and here also useless since every Windows QGIS download automatically also downloads SAGA and GRASS.", notation = "alphabet") +# kableExtra::add_footnote(label = "Comparing downloads of different providers is rather difficult (see http://spatialgalaxy.net/2011/12/19/qgis-users-around-the-world), and here also useless since every Windows QGIS download automatically also downloads SAGA and GRASS GIS.", notation = "alphabet") ``` また、R-GIS ブリッジを補完するために、３つのブリッジの後で空間ライブラリへのインタフェース (Section \@ref(gdal))、空間データベース\index{spatial database} (Section \@ref(postgis))、地球観測データのクラウド処理 (Section \@ref(cloud)) について簡単に紹介する。 ## **qgisprocess**: QGIS へのブリッジなど {#rqgis} -QGIS\index{QGIS} は、最も人気のあるオープンソース GIS である [Table \@ref(tab:gis-comp); @graser_processing_2015]。 -QGIS は、統一されたインターフェイスで、QGIS自身のジオアルゴリズム、GDAL\index{GDAL}、さらにインストールされている場合にはGRASS\index{GRASS}、SAGA\index{SAGA} を利用することができる [@graser_processing_2015]。 +QGIS\index{QGIS} は、最も人気のあるオープンソース GIS である (Table \@ref(tab:gis-comp); @graser_processing_2015)。 +QGIS は、統一されたインターフェイスで、QGIS 自身のジオアルゴリズム、GDAL\index{GDAL}、さらにインストールされている場合には GRASS GIS\index{GRASS GIS}、SAGA\index{SAGA} などのプロバイダを利用することができる [@graser_processing_2015]。 バージョン 3.14 (2020年夏にリリース) 以降、QGIS は、さまざまジオコンピュテーション機能にアクセスできる `qgis_process` コマンドラインを提供している。 -`qgis_process` は、標準的 QGIS に備わる 300 以上のジオアルゴリズムと、GRASS や SAGA など 1000 以上の外部プロバイダへのアクセスを提供している。 +`qgis_process` は、標準的 QGIS に備わる 300 以上のジオアルゴリズムと、GRASS GIS や SAGA など 1000 以上の外部プロバイダへのアクセスを提供している。 **qgisprocess**\index{qgisprocess (package)} パッケージは、R からアクセスすることも可能である。 -このパッケージは、 システム上に最低でも QGIS、また本章で使用する関連するプラグインである GRASS や SAGA を必要とする。 +このパッケージは、 システム上に最低でも QGIS、また本章で使用する関連するプラグインである GRASS GIS や SAGA を必要とする。 インストールに関しては、**qgisprocess** [ドキュメント](https://r-spatial.github.io/qgisprocess/)を参照。 -あるいは、Docker をインストール済みの方は、本プロジェクトの `qgis` イメージから **qgisprocess** を使うこともできる。 + +```{block2 qgisdocker, type='rmdnote'} +あるいは、Docker がインストール済みであれば、本プロジェクトの `qgis` イメージから **qgisprocess** を使うこともできる。 Docker をインストールしており、十分な実行環境のある方は、以下のコマンドで **qgisprocess** および関連プラグインを実行できる ([geocompx/docker](https://github.com/geocompx/docker) レポジトリを参照)。 -```{bash, eval=FALSE} -docker run -e DISABLE_AUTH=true -p 8786:8787 ghcr.io/geocompx/docker:qgis +`docker run -e DISABLE_AUTH=true -p 8786:8787 ghcr.io/geocompx/docker:qgis` + ``` ```{r, eval=has_qgis_plugins} @@ -156,7 +161,7 @@ qgis_plugins() #> 4 processing_saga_nextgen FALSE ``` -プラグイン GRASS (`grassprovider`) と SAGA (`processing_saga_nextgen`) が存在しているが、有効になっていないことがわかる。 +プラグイン GRASS GIS (`grassprovider`) と SAGA (`processing_saga_nextgen`) が存在しているが、有効になっていないことがわかる。 この二つは本章の後半で使用するので、有効化しよう。 ```{r enable-providers, eval=has_qgis_plugins} @@ -183,11 +188,11 @@ qgis_providers() #> 7 sagang SAGA Next Gen 509 ``` -出力表から、QGIS のジオアルゴリズム (`native`, `qgis`, `3d`) と、サードパーティプロバイダの GDAL、SAGA、GRASS の外部アルゴリズムを QGIS インターフェースを通して使用できることが確認できた。 +出力表から、QGIS のジオアルゴリズム (`native`, `qgis`, `3d`) と、サードパーティプロバイダの GDAL、SAGA、GRASS GIS の外部アルゴリズムを QGIS インターフェースを通して使用できることが確認できた。 これで、R から QGIS などの地理計算をする準備ができた。 -それでは、2つの事例を試してみよう。 -最初のものは、異なる境界線を持つ2つのポリゴンデータセットを和集合を作成 (union)\index{union} する方法を示している (Section \@ref(qgis-vector))。 +それでは、2 つの事例を試してみよう。 +最初のものは、異なる境界線を持つ 2 つのポリゴンデータセットを和集合を作成 (union)\index{union} する方法を示している (Section \@ref(qgis-vector))。 もう一つは、ラスタ (Section \@ref(qgis-raster)) で表現された数値標高モデルから新しい情報を導き出すことに重点を置いている。 ### ベクタデータ {#qgis-vector} @@ -203,7 +208,7 @@ incongr_wgs = st_transform(incongruent, "EPSG:4326") aggzone_wgs = st_transform(aggregating_zones, "EPSG:4326") ``` -```{r uniondata, echo=FALSE, fig.cap="二つの単位の図: 不一致 (黒い線) と集合ゾーン (赤い境界)。"} +```{r uniondata, echo=FALSE, fig.cap="二つの単位の図: 不一致 (黒い線) と集合ゾーン (赤い境界)。", fig.width=8} #| message: FALSE #| results: hide library(tmap) @@ -214,7 +219,8 @@ tm_shape(incongr_wgs) + tm_add_legend(type = "lines", labels = c("incongr_wgs", "aggzone_wgs"), col = c("grey5", "red"), lwd = 3, position = tm_pos_in("right", "top")) + tm_scalebar(position = c("left", "bottom"), breaks = c(0, 0.5, 1)) + - tm_layout(frame = FALSE, legend.text.size = 1) + tm_layout(frame = FALSE, legend.text.size = 1, + inner.margins = c(0.02, 0.02, 0.03, 0.02)) ``` 最初に、二つのベクタをマージ (merge) するアルゴリズムを探そう。 @@ -271,9 +277,9 @@ union_arguments これはアルゴリズムの実行時間を長くしてしまう。 **qgisprocess** の主な機能は `qgis_run_algorithm()` であり、QGIS に入力を送り、出力を受け取る。 -これは、使用するアルゴリズム名とヘルプに示される名前付き引数のセットを受け取り、期待される計算を実行する。 -今回のケースでは、`INPUT`、`OVERLAY`、`OUTPUT` の3つの引数が重要だと思われる。 -最初の `INPUT` は、主なベクタオブジェクト `incongr_wgs` であり、2番目の `OVERLAY` は、`aggzone_wgs` である。 +これは、アルゴリズム名とヘルプに示される名前付き引数のセットを受け取り、期待される計算を実行する。 +今回のケースでは、`INPUT`、`OVERLAY`、`OUTPUT` の 3 つの引数が重要だと思われる。 +最初の `INPUT` は、主なベクタオブジェクト `incongr_wgs` であり、2 番目の `OVERLAY` は、`aggzone_wgs` である。 最後の引数、`OUTPUT` は出力ファイル名だが、指定されていない場合、**qgisprocess** は自動的に `tempdir()` に一時ファイルを作成する。 ```{r 09-gis-7, eval=has_qgis_plugins} @@ -292,7 +298,7 @@ union union_sf = st_as_sf(union) ``` -QGIS\index{QGIS} の和集合 (union \index{vector!union}) の操作は、2つの入力レイヤーの交差 (intersect)\index{vector!intersection} と対称差 (symmetrical difference) を用いて、2つの入力レイヤを1つのレイヤにマージすることに注意 (ちなみに、これは GRASS\index{GRASS} と SAGA\index{SAGA} でユニオン操作をするときのデフォルトでもある)。 +QGIS\index{QGIS} の和集合 (union \index{vector!union}) の操作は、2つの入力レイヤーの交差 (intersect)\index{vector!intersection} と対称差 (symmetrical difference) を用いて、2つの入力レイヤを1つのレイヤにマージすることに注意 (ちなみに、これは GRASS GIS\index{GRASS GIS} と SAGA\index{SAGA} でユニオン操作をするときのデフォルトでもある)。 これは `st_union(incongr_wgs, aggzone_wgs)` とは**違う** (演習参照)！ その結果である `union_sf` は、2つの入力オブジェクトよりも多くのフィーチャを持つポリゴンとなる。 @@ -303,9 +309,13 @@ QGIS\index{QGIS} の和集合 (union \index{vector!union}) の操作は、2つ ```{r, eval=has_qgis_plugins} qgis_search_algorithms("clean") +#> # A tibble: 1 × 5 +#> provider provider_title group algorithm algorithm_title +#> +#> 1 grass7 GRASS Vector (v.*) grass7:v.clean v.clean ``` -今回見つかったアルゴリズム (`v.clean`) は、QGIS ではなく、GRASS GIS\index{GRASS} に含まれている。 +今回見つかったアルゴリズム (`v.clean`) は、QGIS ではなく、GRASS GIS\index{GRASS GIS} に含まれている。 GRASS GIS の `v.clean` は、空間ベクタデータのトポロジーをクリーニング\index{topology cleaning}する強力なツールである。 重要なのは、**qgisprocess** を通して使用できることである。 @@ -388,7 +398,7 @@ qgis_search_algorithms("wetness") |> #> 2 SAGA Next Gen sagang:topographicwetnessindexonestep ``` -上記のコードの出力は、目的のアルゴリズムが SAGA GIS\index{SAGA} ソフトウェアに存在することを示唆するものである。^[TWI は、 GRASS GIS 関数 `r.topidx` でも計算可能。] +上記のコードの出力は、目的のアルゴリズムが SAGA\index{SAGA} ソフトウェアに存在することを示唆するものである。^[TWI は、GRASS GIS 関数 `r.topidx` でも計算可能。] SAGA はハイブリッド GIS であるが、主にラスタ処理、ここでは特にデジタル標高モデル \index{digital elevation model} (土壌特性、地形属性、気候パラメータ) に重点を置いている。 したがって、SAGA は大規模な (高解像度の) ラスタ\index{raster}データセットの高速処理に特に優れている [@conrad_system_2015]。 @@ -396,14 +406,11 @@ SAGA はハイブリッド GIS であるが、主にラスタ処理、ここで ```{r, eval=FALSE} qgis_show_help("sagang:sagawetnessindex") -#> Saga wetness index (sagang:sagawetnessindex) -#> Description ``` ここでは、デフォルトの引数を使用する。 与える引数は、入力となる `DEM` だけである。 -このアルゴリズムを使う際は、デフォルト値が研究の目的にあっているか確認する必要がある。^[`"sagang:sagawetnessindex"` の追加引数は、 https://gis.stackexchange.com/a/323454/20955 で詳しく解説されている。] -`"sagang:sagawetnessindex"` は、集水域、集水勾配、修正集水域、地形湿潤指数という4つのラスタを返す。 +このアルゴリズムを使う際は、パラメータ値が研究の目的にあっているか確認する必要がある。^[`"sagang:sagawetnessindex"` の追加引数は、https://gis.stackexchange.com/a/323454/20955 で詳しく解説されている。] QGIS から SAGA アルゴリズムを使う前に、デフォルトのラスタ形式を `.tif` から、SAGA のデフォルトの `.sdat` に変更しておこう。 これで、指定しなければ保存形式は `.sdat` となる。 @@ -416,8 +423,8 @@ dem_wetness = qgis_run_algorithm("sagang:sagawetnessindex", ) ``` - -`qgis_as_terra()` 関数で出力名を指定することで、選択した出力を読み出すことができる。 +`"sagang:sagawetnessindex"` は、集水域、集水勾配、修正集水域、地形湿潤指数という 4 つのラスタを返す。`qgis_as_terra()` 関数で出力名を指定することで、選択した出力を読み出すことができる。 +選択された出力は、`qgis_as_terra()` 関数に出力名を与えることで読むことができる。 QGIS から SAGA を使う作業は終わったので、デフォルト形式を `.tif` に戻す。 ```{r, eval=has_qgis_plugins} @@ -453,7 +460,7 @@ dem_geomorph = qgis_run_algorithm("grass7:r.geomorphon", ``` 出力される `dem_geomorph$forms` は、10個のカテゴリーからなるラスタファイルで、それぞれが地形形状を表している。 -これを `qgis_as_terra()` で R に読み込んで可視化したり (Figure \@ref(fig:qgis-raster-map) の右パネル)、その後の計算で使うことができる。 +これを `qgis_as_terra()` で R に読み込んで可視化したり (Figure \@ref(fig:qgis-raster-map) 右図)、その後の計算で使うことができる。 ```{r, eval=has_qgis_plugins} dem_geomorph_terra = qgis_as_terra(dem_geomorph$forms) @@ -466,13 +473,13 @@ TWI 値が最も大きいのは谷や窪地であり、最も小さいのは予 knitr::include_graphics("figures/10-qgis-raster-map.png") ``` -## SAGA GIS {#saga} +## SAGA {#saga} System for Automated Geoscientific Analyses (SAGA \index{SAGA} ; Table \@ref(tab:gis-comp) ) は、コマンドラインインタフェース \index{command-line interface} ( Windows では `saga_cmd.exe`、Linux では単に `saga_cmd` ) を介して SAGA モジュールを実行する可能性を提供する ([SAGA wiki on modules](https://sourceforge.net/p/saga-gis/wiki/Executing%20Modules%20with%20SAGA%20CMD/) を参照)。 また、Pythonインターフェース (SAGA Python API \index{API} ) も用意されている。 **Rsagacmd**\index{Rsagacmd (package)} は、前者を使って R 内で SAGA\index{SAGA} を実行している。 -この Section では、**Rsagacmd** を使用して、SAGA GIS の seeded region growing アルゴリズムを使って、2000年9月22日のペルーの Mongón 調査地域の正規化差分植生指数 (normalized difference vegetation index, NDVI) の値が類似した地域を抽出する (Figure \@ref(fig:sagasegments) 左)。^[リモートセンシング画像から NDVI を算出する方法については、Section \@ref(local-operations) を参照。] +この Section では、**Rsagacmd** を使用して、SAGA の seeded region growing アルゴリズムを使って、2000年9月22日のペルーの Mongón 調査地域の正規化差分植生指数 (normalized difference vegetation index, NDVI) の値が類似した地域を抽出する (Figure \@ref(fig:sagasegments) 左図)。^[リモートセンシング画像から NDVI を算出する方法については、Section \@ref(local-operations) を参照。] ```{r, eval=FALSE} ndvi = rast(system.file("raster/ndvi.tif", package = "spDataLarge")) @@ -481,7 +488,7 @@ ndvi = rast(system.file("raster/ndvi.tif", package = "spDataLarge")) **Rsagacmd** を始めるには、`saga_gis()` 関数を実行する必要がある。 この関数は主に2つの目的がある。 -- 有効な SAGA-GIS ライブラリやツールへのリンクを含む新しいオブジェクトを動的に作成すること^[つまり、SAGA GIS バージョンによってう使用できるライブラリが変わる。] +- 有効な SAGA ライブラリやツールへのリンクを含む新しいオブジェクトを動的に作成すること^[つまり、SAGA バージョンによってう使用できるライブラリが変わる。] - `raster_backend` (ラスタデータを扱う際に用いる R パッケージ)、`vector_backend` (ベクタデータを扱う際に用いる R パッケージ)、`cores` (処理に用いる CPU コアの最大数、デフォルトは all) など、一般的なパッケージオプションを設定すること ```{r, eval=FALSE} @@ -514,7 +521,7 @@ ndvi_seeds = sg(ndvi, band_width = 2) この出力は、3つのオブジェクトからなるリストである。`variance` -- 局所分散のラスタマップ、 `seed_grid` -- 生成されたシードを含むラスタマップ、 `seed_points` -- 生成されたシードを含む空間ベクタオブジェクト。 -2つ目の SAGA GISツールは `seeded_region_growing` である。^[詳細は https://saga-gis.sourceforge.io/saga_tool_doc/8.3.0/imagery_segmentation_3.html を参照] +2 つ目の SAGA ツールは `seeded_region_growing` である。^[詳細は https://saga-gis.sourceforge.io/saga_tool_doc/8.3.0/imagery_segmentation_3.html を参照] `seed_region_growing` ツールは、前のステップで計算した `seed_grid` と `ndvi` ラスタオブジェクトの 2 つの入力を必要とする。 さらに、入力フィーチャを標準化するための `normalize` や `neighbour` (4または8-neighborhood)、 `method` などのパラメータを指定することができる。 最後のパラメータには、`0` または `1` を指定することができる (ラスタセルの値とその位置に基づいて領域を成長させるか、値のみを成長させるか)。 @@ -530,7 +537,7 @@ plot(ndvi_srg$segments) このツールは、3つのオブジェクトのリストを返す。このツールは、`segments`、`similarity`、`table` という 3 つのオブジェクトのリストを返す。 `similarity` オブジェクトは、シードと他のセルとの類似性を示すラスタであり、`table` は入力シードに関する情報を格納したデータフレームである。 -最後に、`ndvi_srg$segments` は、結果として得られた領域 (Figure \@ref(fig:sagasegments) の右側のパネル) を表すラスタである。 +最後に、`ndvi_srg$segments` は、結果として得られた領域 (Figure \@ref(fig:sagasegments) 右図) を表すラスタである。 これをポリゴンに変換するには、`as.polygons()` と `st_as_sf()` を使用する (Section \@ref(spatial-vectorization))。 ```{r, eval=FALSE} @@ -552,35 +559,35 @@ R には、似たような値を持つポリゴン (いわゆるセグメント) ## GRASS GIS {#grass} -米国陸軍建設工学研究所 (U.S. Army - Construction Engineering Research Laboratory, USA-CERL) は、1982年から1995年にかけて、地理資源解析支援システム (Geographical Resources Analysis Support System, GRASS)\index{GRASS} の中核となるシステムを作成した [Table \@ref(tab:gis-comp); @neteler_open_2008]。 +米国陸軍建設工学研究所 (U.S. Army - Construction Engineering Research Laboratory, USA-CERL) は、1982年から1995年にかけて、地理資源解析支援システム (Geographical Resources Analysis Support System, GRASS GIS)\index{GRASS GIS} の中核となるシステムを作成した [Table \@ref(tab:gis-comp); @neteler_open_2008]。 アカデミアは1997年からこの作業を継続した。 -SAGA\index{SAGA} と同様、グラスも当初はラスタ処理に注力し、その後、GRASS 6.0 以降、高度なベクタ機能を追加している [@bivand_applied_2013]。 +SAGA\index{SAGA} と同様、グラスも当初はラスタ処理に注力し、その後、GRASS GIS 6.0 以降、高度なベクタ機能を追加している [@bivand_applied_2013]。 -GRASS は、入力データを GRASS GIS データベースに格納する。 -ベクタデータに関して、GRASS はデフォルトでトポロジカル GIS、すなわち隣接するフィーチャのジオメトリを一度だけ保存する。 +GRASS は、入力データを内部データベースに格納する。 +ベクタデータに関して、GRASS GIS はデフォルトでトポロジカル GIS、すなわち隣接するフィーチャのジオメトリを一度だけ保存する。 ベクタ属性の管理にはデフォルトで SQLite を用い、属性はキーによってジオメトリ、すなわち GRASS GIS データベースにリンクされる ([GRASS GIS vector management](https://grasswiki.osgeo.org/wiki/Vector_Database_Management#GRASS_GIS_vector_management_model))。 -GRASS を使う前に、GRASS GIS データベース\index{spatial database}を (R からも) セットアップする必要があるが、このプロセスに少し戸惑うかもしれない。 -まず、GRASS のデータベースは専用のディレクトリを必要とし、そのディレクトリには location を置く必要がある (詳しくは [grass.osgeo.org](https://grass.osgeo.org/grass-stable/manuals/index.html) の [GRASS GIS Database](https://grass.osgeo.org/grass-stable/manuals/grass_database.html) ヘルプページを参照)。 +GRASS GIS を使う前に、GRASS GIS データベース\index{spatial database}を (R からも) セットアップする必要があるが、このプロセスに少し戸惑うかもしれない。 +まず、GRASS GIS のデータベースは専用のディレクトリを必要とし、そのディレクトリには location を置く必要がある (詳しくは [grass.osgeo.org](https://grass.osgeo.org/grass-stable/manuals/index.html) の [GRASS GIS Database](https://grass.osgeo.org/grass-stable/manuals/grass_database.html) ヘルプページを参照)。 location には、1 つのプロジェクトまたは 1 つの領域のジオデータが格納される。 -通常、1つの場所の中に、異なるユーザや異なるタスクを参照するいくつかのマップセットを存在させることができる。 +通常、1 つの場所の中に、異なるユーザや異なるタスクを参照するいくつかのマップセットを存在させることができる。 各 location には、PERMANENT マップセット (自動的に作成される必須のマップセット) もある。 プロジェクトのすべてのユーザと地理データを共有するために、データベース所有者は PERMANENT マップセットに空間データを追加することができる。 さらに、PERMANENT マップセットには、ラスタデータの投影法、空間範囲、およびデフォルトの解像度が格納される。 まとめると、GRASS GIS データベースは多くの location を含み (1つのロケーションのデータはすべて同じ CRS を持つ)、それぞれの location は多くのマップセット (データセットのグループ) を格納することができる。 -GRASS 空間データベース\index{spatial database}システムの詳細は、@neteler_open_2008 と [GRASS GIS quick start](https://grass.osgeo.org/grass-stable/manuals/helptext.html) を参照。 -こkでは R から手軽に GRASS を使うため **link2GI** パッケージを使う。しかし、GRASS GIS データベースを順番に作ることもできる。 +GRASS GIS 空間データベース\index{spatial database}システムの詳細は、@neteler_open_2008 と [GRASS GIS quick start](https://grass.osgeo.org/grass-stable/manuals/helptext.html) を参照。 +こkでは R から手軽に GRASS GIS を使うため **link2GI** パッケージを使う。しかし、GRASS GIS データベースを順番に作ることもできる。 作り方は [GRASS within R](https://grasswiki.osgeo.org/wiki/R_statistics/rgrass#GRASS_within_R) を参照。 -以下の段落で解説するコードは、初めて GRASS を使う方には難しいかもしれないが、コードを１行１行実行して途中の結果を確認することで、コードの元となる理由が明らかになる。 +以下の段落で解説するコードは、初めて GRASS GIS を使う方には難しいかもしれないが、コードを 1 行 1 行実行して途中の結果を確認することで、コードの元となる理由が明らかになる。 ここでは、GIScience における最も興味深い問題の一つである巡回セールスマン問題 \index{traveling salesman} を用いた **rgrass** \index{rgrass (package)} を紹介する。 ある巡回セールスマンが24件の顧客を訪問したいとする。 さらに、自宅を起点とし、最短距離で25カ所を回りたいということであった。 この問題に対する最適解は一つであるが、考えられる解をすべてチェックすることは、現代のコンピュータでは (ほとんど) 不可能である [@longley_geographic_2015]。 -この場合、可能な解の数は `(25 - 1)! / 2`、すなわち24の階乗を2で割った数に相当する (2で割るのは、順方向と逆方向を区別しないため)。 -1回の繰り返しがナノ秒でも、`r format(factorial(25 - 1) / (2 * 10^9 * 3600 * 24 * 365))` 年間に相当する。 +この場合、可能な解の数は `(25 - 1)! / 2`、すなわち24の階乗を2で割った数に相当する (2 で割るのは、順方向と逆方向を区別しないため)。 +1 回の繰り返しがナノ秒でも、`r format(factorial(25 - 1) / (2 * 10^9 * 3600 * 24 * 365))` 年間に相当する。 幸いなことに、この想像を絶する時間のごく一部で実行できる、巧妙でほぼ最適なソリューションがある。 -GRASS GIS\index{GRASS} は、これらの解決策の一つを提供する (詳細は、 [v.net.salesman](https://grass.osgeo.org/grass-stable/manuals/v.net.salesman.html)を参照)。 +GRASS GIS\index{GRASS GIS} は、これらの解決策の一つを提供する (詳細は、 [v.net.salesman](https://grass.osgeo.org/grass-stable/manuals/v.net.salesman.html)を参照)。 今回の使用例では、ロンドンの街角にある最初の25の自転車ステーション (顧客の代わり) 間の最短経路\index{shortest route}を見つけたい (最初の自転車ステーションは、巡回セールスマン\index{traveling salesman}の自宅に相当すると単純に仮定する)。 ```{r 09-gis-24} @@ -603,11 +610,11 @@ london_streets = london_streets[["osm_lines"]] london_streets = select(london_streets, osm_id) ``` -これでデータが揃ったので、次に GRASS\index{GRASS} のセッションを開始する。 -幸い、**link2GI** パッケージの `linkGRASS()` を使えば、たった一行のコードで GRASS 環境をセットアップできる。 +これでデータが揃ったので、次に GRASS GIS\index{GRASS GIS} のセッションを開始する。 +幸い、**link2GI** パッケージの `linkGRASS()` を使えば、たった一行のコードで GRASS GIS 環境をセットアップできる。 空間オブジェクトは、空間データベースの投影と範囲を決定するものである。 -まず、`linkGRASS()` は、あなたのコンピュータにインストールされている全ての GRASS\index{GRASS} を検索する。 -ここでは `ver_select` を `TRUE` に設定しているので、見つかった GRASS-installation の中から対話的に一つを選択することができる。 +まず、`linkGRASS()` は、あなたのコンピュータにインストールされている全ての GRASS GIS\index{GRASS GIS} を検索する。 +ここでは `ver_select` を `TRUE` に設定しているので、見つかった GRASS GIS-installation の中から対話的に一つを選択することができる。 もし、インストールが一つしかない場合は、`linkGRASS()` が自動的にそれを選択する。 次に、`linkGRASS()` は、GRASS GIS への接続を確立する。 @@ -616,10 +623,10 @@ library(rgrass) link2GI::linkGRASS(london_streets, ver_select = TRUE) ``` -GRASS の geoalgorithm を使用する前に、GRASS の空間データベースにデータを追加する必要がある。 +GRASS GIS の geoalgorithm を使用する前に、GRASS GIS の空間データベースにデータを追加する必要がある。 幸いなことに、便利な関数 `write_VECT()` がこれを代行してくれる。 (ラスタデータには `write_RAST()` を使用する。) -この例では、最初の属性カラムのみを使用して、道路と自転車レンタルポイントデータを追加し、GRASS で `london_streets` と `points` という名前を付けている。 +この例では、最初の属性カラムのみを使用して、道路と自転車レンタルポイントデータを追加し、GRASS GIS で `london_streets` と `points` という名前を付けている。 ```{r 09-gis-31, eval=FALSE} write_VECT(terra::vect(london_streets), vname = "london_streets") @@ -632,8 +639,8 @@ write_VECT(terra::vect(points[, 1]), vname = "points") 現在、両方のデータセットが GRASS GIS のデータベースに存在している。 ネットワーク\index{network}の解析を行うには、トポロジカルクリーン\index{topology cleaning}な道路ネットワークが必要である。 -GRASS の `"v.clean"` は、重複、小角、ダングルの除去などを行う。 -ここでは、後続のルーティングアルゴリズムが実際に交差点で右折または左折できるように、各交差点で改行し、その出力を `streets_clean` という名前の GRASS オブジェクトに保存している。 +GRASS GIS の `"v.clean"` は、重複、小角、ダングルの除去などを行う。 +ここでは、後続のルーティングアルゴリズムが実際に交差点で右折または左折できるように、各交差点で改行し、その出力を `streets_clean` という名前の GRASS GIS オブジェクトに保存している。 ```{r 09-gis-32, eval=FALSE} execGRASS( @@ -663,8 +670,8 @@ execGRASS( 得られたクリーンなデータセットは `"v.net.salesman"` アルゴリズムの入力となり、最終的にすべての自転車レンタルステーション間の最短経路を見つけることができる。 その引数の一つが `center_cats` で、これは入力として数値の範囲を必要とする。 この範囲は、最短ルートを計算するためのポイントを表している。 -ここでは、すべての自転車ステーション間の経路を計算したいので、`1-25`に設定しておく。 -巡回セールスマンアルゴリズムの GRASS ヘルプページを参照するには、 `execGRASS("g.manual", entry = "v.net.salesman")` を実行する。 +ここでは、すべての自転車ステーション間の経路を計算したいので、`1-25` に設定しておく。 +巡回セールスマンアルゴリズムの GRASS GIS ヘルプページを参照するには、`execGRASS("g.manual", entry = "v.net.salesman")` を実行する。 ```{r 09-gis-33, eval=FALSE} execGRASS( @@ -674,11 +681,7 @@ execGRASS( ) ``` -作成したクリーンなデータセットは `"v.net.salesman"` アルゴリズムの入力となる。 -その引数の一つは `center_cats` であり、入力として数値範囲を必要とする。 -この範囲は最短ルートを計算する点を表す。 -すべてのサイクルステーションのルートを計算したいので、`1-25` に設定する。 -巡回セールスマンの GRASS ヘルプページにアクセスするには、`execGRASS("g.manual", entry = "v.net.salesman")`を実行する。 +結果を見るには、結果を R に入れ、ジオメトリのみ保持した `sf` オブジェクトに変換する。これを **mapview** で可視化する (Figure \@ref(fig:grass-mapview) と Section \@ref(interactive-maps))。 ```{r 09-gis-34, eval=FALSE} route = read_VECT("shortest_route") |> @@ -705,51 +708,40 @@ mapview::mapshot(m_1, その際、いくつか注意すべき点がある。 -- GRASSの空間データベース\index{spatial database}を使えば、より高速に処理できる。 -つまり、最初に地理データを書き出しただけである。 +- GRASS GIS の空間データベース\index{spatial database}を使えば、より高速に処理できた。 +しかし、ここでは最初に地理データを書き出した。 そして、新しいオブジェクトを作成し、最終結果だけを R にインポートした。 現在利用可能なデータセットを調べるには、`execGRASS("g.list", type = "vector,raster", flags = "p")` を実行する。 -- また、R から既にある GRASS 空間データベースにアクセスすることも可能である。 +- また、R から既にある GRASS GIS 空間データベースにアクセスすることも可能であった。 R にデータをインポートする前に、いくつかの (空間) 部分集合を作成したい場合がある。 ベクタデータには `"v.select"` と `"v.extract"` を使用する。 `"db.select"` を使用すると、対応するジオメトリを返さずにベクタレイヤの属性テーブルの部分集合を選択することができる。 -- また、実行中の GRASS のセッションから R を起動することもできる [詳細は @bivand_applied_2013 を参照]。 -- GRASS で提供されて入るジオアルゴリズム\index{geoalgorithm}の素晴らしいドキュメントは、 [GRASS online help](https://grass.osgeo.org/grass-stable/manuals/) または `execGRASS("g.manual", flags = "i")` を参照。 +- また、実行中の GRASS GIS\index{GRASS GIS} のセッションから R を起動することもできる [詳細は @bivand_applied_2013 を参照]。 +- GRASS GIS で提供されて入るジオアルゴリズム\index{geoalgorithm}の素晴らしいドキュメントは、 [GRASS GIS online help](https://grass.osgeo.org/grass-stable/manuals/) または `execGRASS("g.manual", flags = "i")` を参照。 ## いつ、何を使うべきか？ R-GIS のインターフェースは、個人の好みや作業内容、GIS\index{GIS} の使い方に依存するため、一概にお勧めすることはできないし、研究分野にもよるだろう。 -前述の通り、SAGA \index{SAGA} は大規模 (高解像度) ラスタ\index{raster}データセットの高速処理に特に優れており、水文学者、気候学者、土壌学者に頻繁に利用されている [@conrad_system_2015]。 -一方、GRASS GIS\index{GRASS} は、トポロジーに基づく空間データベースをサポートする唯一の GIS であり、ネットワーク分析だけでなくシミュレーション研究にも特に有用である。 -QGIS\index{QGIS} は、GRASS-GIS や SAGA-GIS と比較して、特に初めて GIS を使う方にとって使いやすく、おそらく最も人気のあるオープンソースのGISだと思われる。 +前述の通り、SAGA\index{SAGA} は大規模 (高解像度) ラスタ\index{raster}データセットの高速処理に特に優れており、水文学者、気候学者、土壌学者に頻繁に利用されている [@conrad_system_2015]。 +一方、GRASS GIS\index{GRASS GIS} は、トポロジーに基づく空間データベースをサポートする唯一の GIS であり、ネットワーク分析だけでなくシミュレーション研究にも特に有用である。 +QGIS\index{QGIS} は、GRASS GIS や SAGA と比較して、特に初めて GIS を使う方にとって使いやすく、おそらく最も人気のあるオープンソースの GIS だと思われる。 したがって、**qgisprocess**\index{qgisprocess (package)} は、ほとんどのユースケースに適切な選択である。 その主なメリットは -- 複数の GIS に統一的にアクセスできるため、重複した機能を含む 1000 以上のジオアルゴリズム ( Table \@ref(tab:gis-comp) ) を提供。例えば、QGIS\index{QGIS}、SAGA\index{SAGA}、GRASS\index{GRASS} などのジオアルゴリズムを使ってオーバーレイ操作を実行することが可能である。 -- データフォーマットの自動変換 (SAGAは `.sdat` グリッドファイル、GRASS は独自のデータベースフォーマットを使用するが、対応する変換は QGIS が行う。) +- 複数の GIS に統一的にアクセスできるため、重複した機能を含む 1000 以上のジオアルゴリズム ( Table \@ref(tab:gis-comp) ) を提供。例えば、QGIS\index{QGIS}、SAGA\index{SAGA}、GRASS GIS\index{GRASS GIS} などのジオアルゴリズムを使ってオーバーレイ操作を実行することが可能である。 +- データフォーマットの自動変換 (SAGAは `.sdat` グリッドファイル、GRASS GIS は独自のデータベースフォーマットを使用するが、対応する変換は QGIS が行う。) - 地理的な R オブジェクトを QGIS ジオアルゴリズム\index{geoalgorithm}に自動的に渡し、R に戻すことができる。 - 名前付き引数、デフォルト値の自動取得をサポートする便利な機能 (**rgrass**\index{rgrass (package)} からインスパイアされた) もちろん、他の R-GIS ブリッジを使用した方が良いケースもある。 QGIS は、複数の GIS\index{GIS} ソフトウェアパッケージへの統一インターフェースを提供する唯一の GIS であるが、対応するサードパーティのジオアルゴリズムのサブセットへのアクセスしか提供しない (詳細については、@muenchow_rqgis:_2017 を参照)。 -したがって、SAGA と GRASS の関数一式を使用するには、**RSAGA**\index{RSAGA (package)} と **rgrass** 以外は使わない方が良い。 -また、ジオデータベース\index{spatial database}を用いてシミュレーションを行いたい場合 [@krug_clearing_2010]、**qgisprocess** が、呼び出しごとに常に新しい GRASS セッションを開始するので、**rgrass** を直接使用してみよう。 -最後に、地形データ、空間データベース管理機能 (マルチユーザーアクセスなど) が必要な場合は、GRASSの利用を勧める。 +したがって、SAGA と GRASS GIS の関数一式を使用するには、**RSAGA**\index{RSAGA (package)} と **rgrass** 以外は使わない方が良い。 +また、ジオデータベース\index{spatial database}を用いてシミュレーションを行いたい場合 [@krug_clearing_2010]、**qgisprocess** が、呼び出しごとに常に新しい GRASS GIS セッションを開始するので、**rgrass** を直接使用してみよう。 +最後に、地形データ、空間データベース管理機能 (マルチユーザーアクセスなど) が必要な場合は、GRASS GIS の利用を勧める。 -なお、スクリプティング・インターフェースを持つ GIS ソフトウェアパッケージは以下のように数多くあるが、これらを利用できる専用の R パッケージはない：gvSig、OpenJump、Orfeo Toolbox。^[注記：**link2GI** は Orfeo Toolbox\index{Orfeo Toolbox} を部分的に統合しており、**qgisprocess** から Orfeo Toolbox アルゴリズムへアクセスすることもできる。TauDEM\index{TauDEM} は R パッケージの **traudem** からアクセスできる。] +なお、スクリプティング・インターフェースを持つ GIS ソフトウェアパッケージは以下のように数多くあるが、これらを利用できる専用の R パッケージはない: gvSig、OpenJump、Orfeo Toolbox。^[注記: **link2GI** は Orfeo Toolbox\index{Orfeo Toolbox} を部分的に統合しており、**qgisprocess** から Orfeo Toolbox アルゴリズムへアクセスすることもできる。TauDEM\index{TauDEM} は R パッケージの **traudem** からアクセスできる。] -## その他のブリッジ - -この章では、当面、デスクトップ GIS \index{GIS} ソフトウェアへの R インタフェースに焦点を当てる。 -これらの橋を強調したのは、専用の GIS ソフトウェアがよく知られており、地理データを理解するための一般的な「入り口」であることがある。 -また、多くのジオアルゴリズム\index{geoalgorithm}にアクセスすることができる。 - -その他の「ブリッジ」には、空間ライブラリへのインタフェース (Section \@ref(gdal) は R から GDAL\index{GDAL} CLI\index{command-line interface} にアクセスする方法を示している)、空間データベース\index{spatial database} ( Section \@ref(postgis) 参照)、ウェブマッピングサービス ( Chapter \@ref(adv-map) 参照) などがある。 -このセクションでは、可能なことのほんの一端を紹介する。 -システムから他のプログラムを呼び出したり、他の言語と統合したり (特に **Rcpp** と **reticulate**\index{reticulate (package)} を介して)、R\index{R} の柔軟性のおかげで、他の多くのブリッジが可能になっている。 -目的は包括的なものではなく、章の冒頭で @sherman_desktop_2008 が引用した「柔軟性とパワー」にアクセスする他の方法を示すことである。 - -### GDAL へのブリッジ {#gdal} +## GDAL へのブリッジ {#gdal} Chapter \@ref(read-write) で述べたように、GDAL\index{GDAL} は多くの地理データ形式をサポートする低レベルのライブラリである。 GDAL は非常に効果的なので、ほとんどのGISプログラムは、車輪の再発明や特注の読み書きコードを使用するのではなく、地理データのインポートとエクスポートのためにバックグラウンドで GDAL\index{GDAL} を使用している。 @@ -794,24 +786,24 @@ system(cmd) GDAL ツールの全リストとそのヘルプファイルは https://gdal.org/programs/ 。 **link2GI** が提供する GDAL への「リンク」は、R やシステムの CLI からより高度な GDAL の作業を行うための基盤として利用することができるだろう。 -TauDEM (http://hydrology.usu.edu/taudem) や Orfeo Toolbox (https://www.orfeo-toolbox.org/) は、コマンドラインインタフェースを提供する空間データ処理ライブラリ/プログラムである -- 上記の例は、Rを介してシステムのコマンドラインからこれらのライブラリにアクセスする方法である。 +TauDEM (http://hydrology.usu.edu/taudem) や Orfeo Toolbox (https://www.orfeo-toolbox.org/) は、コマンドラインインタフェースを提供する空間データ処理ライブラリ/プログラムである -- 上記の例は、R を介してシステムのコマンドラインからこれらのライブラリにアクセスする方法である。 これは、新しい R パッケージという形で、これらのライブラリへの適切なインタフェースを作成するための出発点となる可能性がある。 しかし、新しいブリッジを作成するプロジェクトに飛び込む前に、既存の R パッケージのパワーと、`system()` の呼び出しがプラットフォームに依存しない (一部のコンピュータで失敗する) 可能性があることを認識しておくことが重要である。 -さらに、**sf** は GDAL\index{GDAL}、GEOS\index{GEOS}、PROJ\index{PROJ} が提供するパワーのほとんどを **Rcpp** が提供する R/C++ \index{C++} インターフェイスを介して R にもたらし、`system()` の呼び出しを回避しているのである。 +一方、**sf** は GDAL\index{GDAL}、GEOS\index{GEOS}、PROJ\index{PROJ} が提供するパワーのほとんどを **Rcpp** が提供する R/C++\index{C++} インターフェイスを介して R にもたらし、`system()` の呼び出しを回避している。 ### 空間データベースへのブリッジ {#postgis} \index{spatial database} -空間データベース管理システム (空間DBMS) は、空間および非空間データを構造化して保存する。 +空間データベース管理システム (空間 DBMS) は、空間および非空間データを構造化して保存する。 大規模なデータの集合を、一意の識別子 (主キーと外部キー) および暗黙のうちに空間を介して関連するテーブル (エンティティ) に整理することができる (たとえば、空間結合を考えてみてみよう)。 地理的なデータセットはすぐに大きくなったり、乱雑になったりする傾向があるため、この機能は便利である。 データベースは、空間および非空間フィールドに基づく大規模なデータセットの保存とクエリを効率的に行うことができ、マルチユーザーアクセスとトポロジー\index{topological relations}のサポートを提供する。 最も重要なオープンソースの空間データベース\index{spatial database}は PostGIS \index{PostGIS} である [@obe_postgis_2015]。^[ -SQLite/SpatiaLiteも確かに重要であるが、GRASS \index{GRASS} はバックグラウンドでSQLiteを使っているので、暗黙のうちにこの方法をすでに導入している (Section \@ref(grass) を参照)。 +SQLite/SpatiaLite も確かに重要であるが、GRASS GIS\index{GRASS GIS} はバックグラウンドで SQLite を使っているので、暗黙のうちにこの方法をすでに導入している (Section \@ref(grass) を参照)。 ] -PostGIS \index{PostGIS} のような空間 DBMS への R ブリッジは重要で、数ギガバイトの地理データを RAM にロードすることなく、、R セッションをクラッシュさせる可能性があるような巨大なデータストアにアクセスできる。 +PostGIS\index{PostGIS} のような空間 DBMS への R ブリッジは重要で、数ギガバイトの地理データを RAM にロードすることなく、R セッションをクラッシュさせる可能性があるような巨大なデータストアにアクセスできる。 このセクションの残りの部分では、*PostGIS in Action, Second Edition* の "Hello real world" に基づいて、R から PostGIS を呼び出す方法を紹介する [@obe_postgis_2015]。^[ この例の使用を許可してくださった Manning Publications、Regina Obe、Leo Hsu に感謝する。 ] @@ -936,19 +928,18 @@ tm_shape(buf) + border.alpha = 0.3, labels = "35km buffer") + tm_add_legend(type = "symbols", fill = "#F10C26", labels = "Restaurants", shape = 21) + + tm_scalebar() + tm_layout(frame = FALSE) ``` PostGIS とは異なり、**sf** は空間ベクタデータのみをサポートしている。 - - PostGIS データベースに格納されたラスタデータを照会・操作するには、**rpostgis** パッケージ [@bucklin_rpostgis_2018]、または PostGIS\index{PostGIS} インストールの一部に含まれる `rastertopgsql` などのコマンドラインツールを使用する必要がある。 このサブセクションでは、PostgreSQL/PostGIS の簡単な紹介にとどめる。 それでも、地理的および非地理的データを空間 DBMS で保存しながら、さらなる (地理) 統計解析に必要なそれらのサブセットだけを R のグローバル環境にアタッチするという実践を奨励したい。 提示された SQL クエリのより詳細な説明と PostgreSQL/PostGIS 一般のより包括的な紹介は @obe_postgis_2015 を参照。 PostgreSQL/PostGIS は、非常に難解なオープンソースの空間データベースである。 -しかし、軽量なデータベースエンジンである SQLite/SpatiaLite や、バックグラウンドで SQLite を使用する GRASS\INDEX{GRASS} も同様と言える (Section \@ref(grass) 参照)。 +しかし、軽量なデータベースエンジンである SQLite/SpatiaLite や、バックグラウンドで SQLite を使用する GRASS GIS\INDEX{GRASS GIS} も同様と言える (Section \@ref(grass) 参照)。 データセットが PostgreSQL/PostGIS では大きすぎる場合、大規模な空間データ管理とクエリ性能を必要とする場合、分散コンピューティングシステム上での大規模な地理クエリを検討する価値があるかもしれない。 このようなシステムは本書の範囲外ではあるが、この機能を提供するオープンソースソフトウェアが存在することは触れておく価値がある。 @@ -991,29 +982,26 @@ items = s |> ``` クラウドストレージはローカルのハードディスクとは異なり、従来の画像ファイル形式はクラウドベースのジオプロセシングではうまく機能しない。 -クラウドに最適化された GeoTIFF は、GeoTIFF の一種であり、クラウドストレージから画像の一部分のみを効率的に読み込むことができる。 -そのため、画像の矩形部分や低解像度の画像の読み込みが非常に効率的になる。 +クラウドに最適化された GeoTIFF は、画像の矩形部分や低解像度の画像の読み込みが非常に効率的になる。 [GDAL](https://gdal.org)\index{GDAL} (およびそれを使ったパッケージ) はすでに COG を扱うことができるので、R ユーザーであれば COG を扱うために何かをインストールする必要はない。 ただし、データ提供者のカタログを閲覧する際には、COG が利用可能であることが大きなプラスになることを覚えておこう。 領域が大きい時、要求された画像を扱うのはまだ比較的困難である。それらは異なる地図投影を使用することがあり、空間的に重なることがあり、空間解像度はしばしばスペクトルバンドに依存する。 -**gdalcubes** パッケージ\index{gdalcubes (package)} [@appel_gdalcubes_2019] は、個々の画像から抽象化し、画像コレクションを4次元データキューブ\index{data cube}として作成し処理するために使用することができる。 +**gdalcubes** パッケージ\index{gdalcubes (package)} [@appel_gdalcubes_2019] は、個々の画像から抽象化し、画像コレクションを 4 次元データキューブ\index{data cube}として作成し処理するために使用することができる。 以下のコードは、前回の STAC-API 検索で返された Sentinel-2 画像から、低解像度 (250m) の最大 NDVI コンポジットを作成する最小限の例を示している。 ```{r 09-gdalcubes-example, eval=FALSE} library(gdalcubes) # クラウドカバーで画像をフィルタし、画像コレクションを生成 -collection = stac_image_collection(items$features, - property_filter = function(x) { +cloud_filter = function(x) { x[["eo:cloud_cover"]] < 10 - } -) +} +collection = stac_image_collection(items$features, + property_filter = cloud_filter) # データキューブの範囲、解像度 (250m、毎日)、CRS を定義 -v = cube_view( - srs = "EPSG:3857", extent = collection, dx = 250, dy = 250, - dt = "P1D" -) # "P1D" は ISO 8601 期間文字列 +v = cube_view(srs = "EPSG:3857", extent = collection, dx = 250, dy = 250, + dt = "P1D") # "P1D" は ISO 8601 期間文字列 # データキューブを生成し処理 cube = raster_cube(collection, v) |> select_bands(c("B04", "B08")) |> @@ -1026,11 +1014,11 @@ cube = raster_cube(collection, v) |> クラウドカバーによる画像のフィルタリングを行うために、画像コレクションを作成する際に各 STAC\index{STAC} 結果アイテムに適用されるプロパティフィルタ関数を提供している。 この関数は、画像の利用可能なメタデータを入力リストとして受け取り、関数がTRUEを返す画像のみを考慮するような単一の論理値を返す。 この場合、10% 以上のクラウドカバーがある画像は無視する。 -詳しくは、こちらの[OpenGeoHubサマースクール2021で発表したチュートリアル](https://appelmar.github.io/ogh2021/tutorial.html)を参照。 +詳しくは、こちらの [OpenGeoHub サマースクール 2021 で発表したチュートリアル](https://appelmar.github.io/ogh2021/tutorial.html)を参照。^[STAC\index{STAC} を扱う別のパッケージとして、**rsi** がある。これは、指定した時間と位置の STAC 空間データを取得する。] STAC\index{STAC}、COGs\index{COG}、data cube\index{data cube} を組み合わせて、衛星画像の (大規模) コレクションをクラウド上\index{cloud computing}で解析するクラウドネイティブワークフローを形成する. -これらのツールは、例えば、大規模な地球観測データの土地利用や土地被覆の分類を可能にする **sits** R パッケージ\index{sits (package)}のバックボーンを既に形成している。 -このパッケージは、クラウドサービスで利用可能な画像コレクションからEOデータキューブを構築し、様々な機械学習アルゴリズムを用いてデータキューブの土地分類を実行するものである。 +これらのツールは、例えば、大規模な地球観測データの土地利用や土地被覆の分類を可能にする **sits** パッケージ\index{sits (package)}のバックボーンを既に形成している。 +このパッケージは、クラウドサービスで利用可能な画像コレクションから EO データキューブを構築し、様々な機械学習と真相学習アルゴリズムを用いてデータキューブの土地分類を実行するものである。 **sits** の詳細については、https://e-sensing.github.io/sitsbook/ または関連記事 [@rs13132428] を参照。 ### openEO @@ -1043,7 +1031,7 @@ OpenEO [@schramm_openeo_2021]\index{openEO} は、データ処理のための共 その後、ユーザーは画像コレクションのロード、プロセスの適用と連鎖、ジョブの送信、結果の探索とプロットを行うことができる。 以下のコードは、[openEO platform backend](https://openeo.cloud/) に接続し、利用可能なデータセット、プロセス、出力フォーマットを要求し、Sentinel-2 データから最大 NDVI 画像を計算するプロセスグラフを定義し、最後にバックエンドにログインした後にグラフを実行する。 -openEO\index{openEO} プラットフォームのバックエンドには無料版があり、既存の機関やソーシャルプラットフォームのアカウントから登録することが可能である。 +openEO\index{openEO} プラットフォームのバックエンドには無料版があり、既存の機関やインターネットフォームのアカウントから登録することが可能である。 ```{r 09-openeo-example, eval=FALSE} library(openeo) diff --git a/11-algorithms-ja.Rmd b/11-algorithms-ja.Rmd index e1cc1e0..427ef8e 100644 --- a/11-algorithms-ja.Rmd +++ b/11-algorithms-ja.Rmd @@ -1,9 +1,13 @@ # スクリプト、アルゴリズム、関数 {#algorithms} +```{r, include=FALSE} +source("code/before_script.R") +``` + ## 必須パッケージ {- #prerequisites-11} この章では、主に base R を使用するため、必要なソフトウェアは最小限である。 -これから開発するアルゴリズムの結果を確認するために **sf**\index{sf} パッケージを使用する。 +これから開発するアルゴリズムの結果を確認するために **sf**\index{sf} パッケージだけを使用する。 Chapter \@ref(spatial-class) で紹介した地理クラスと、それを使ってさまざまな入力ファイル形式を表現する方法 (Chapter \@ref(read-write) 参照) について理解していることを前提にしている。 ## イントロダクション {#intro-algorithms} @@ -69,9 +73,9 @@ Rscript code/11-hello.R しかし、守るべき慣習もある。 - 順番に書く：映画の脚本と同じように、スクリプトも「設定」「データ処理」「結果保存」といった明確な順番が必要である (映画でいうところの「始まり」「中間」「終わり」にほぼ相当する)。 -- 他の人(と未来の自分)が理解できるように、スクリプトにコメントを追加してみよう。 +- 他の人(と未来の自分)が理解できるように、スクリプトにコメントを追加してみよう 最低限、コメントにはスクリプトの目的 (Figure \@ref(fig:codecheck) 参照) と (長いスクリプトの場合) セクションに分けることが必要である。 -これは、例えば、RStudio\index{RStudio} で、「折りたたみ可能な」コードセクションの見出しを作成するショートカット `Ctrl+Shift+R` を使って行うことができる。 +これは、例えば、RStudio\index{RStudio} で、「折りたたみ可能な」コードセクションの見出しを作成するショートカット `Ctrl+Shift+R` を使って行うことができる - 特に、スクリプトは再現可能であるべきである。どんなコンピュータでも動作する自己完結型のスクリプトは、調子の良い日に自分のコンピュータでしか動作しないスクリプトよりも有用である。 これには、必要なパッケージを最初に添付し、データを永続的なソース (信頼できるウェブサイトなど) から読み込み、前のステップが実行されたことを確認することが含まれる^[ 前のステップは、コメントまたは `if(!exists("x")) source("x.R")` のような if 文で参照できる (オブジェクト `x` が見つからない場合、スクリプトファイル `x.R` を実行する)。 @@ -95,7 +99,8 @@ knitr::include_graphics("figures/codecheck.png") ジオコンピュテーションのためのスクリプトで特に考慮すべき点は、GDAL など外部ライブラリへの依存が多い点である。実際、Chapter \@ref(read-write) のデータ入出力では GDAL をたくさん使用した。 GIS ソフトウェアの依存は、Chaptger \@ref(gis) で解説したようにより多くの特別なジオアルゴリズムを実行するために必要になる。 地理データを扱うスクリプトは、入力データセットが特定のファイル形式であることを必要とする。 -このような依存関係は、スクリプトのコメントとして、またはスクリプトの一部であるプロジェクトの他の場所で言及されるべきである。例えば、[`11-centroid-alg.R` ](https://github.com/geocompx/geocompr/blob/main/code/11-centroid-alg.R) で示している。 +このような依存関係は、スクリプトのコメントとして、またはスクリプトの一部であるプロジェクトの他の場所でコメントするか、**renv** や Docker などで依存性として記述するべきである。 + 「保守的」プログラミング技術と適切なエラーメッセージは、要件が満たされていない時に依存性を確認し、ユーザと対話する時間を節約する。 R では `if ()` で表される If 文を使用し、特定の条件が揃った時のみにメッセージを送ったりコードを実行するようにコードを書く。 以下のコード例は、特定のファイルがない場合にユーザにメッセージを送る (訳註: メッセージに日本語を使わないことを推奨する。)。 @@ -155,10 +160,10 @@ Chapter \@ref(gis) で見たような ジオアルゴリズム\index{geoalgorith 本節では、視覚化しやすいアプローチを選択する。ポリゴンを多くの三角形に分割し、それぞれの重心\index{centroid}を求める。このアプローチは、他の重心アルゴリズムと並んで @kaiser_algorithms_1993 によって議論された [簡単な説明は @orourke_computational_1998]。 コードを書く前に、このアプローチをさらに個別のタスクに分解するのに役立つ (以降、ステップ 1 からステップ 4 と呼ぶが、これらは模式図や疑似コード\index{pseudocode} として提示すこともできる)。 -1. ポリゴンを連続した三角形に分割する。 -2. 各三角形の重心\index{centroid}を求める。 -3. それぞれの三角形の面積を求める。 -4. 三角形の中心点の面積加重平均\index{centroid}を求める。 +1. ポリゴンを連続した三角形に分割する +2. 各三角形の重心\index{centroid}を求める +3. それぞれの三角形の面積を求める +4. 三角形の中心点の面積加重平均\index{centroid}を求める 一見、簡単そうに見えるが、言葉をコードに変換するには、入力に制約がある場合でも、試行錯誤を繰り返しながら作業を進める必要がある。 このアルゴリズムは、180°以上の内角を持たない凸ポリゴンに対してのみ動作し、星形は使用できない (パッケージの **decido** と **sfdct** は外部ライブラリを使用して非凸ポリゴンを三角測量できる。[geocompx.org](https://geocompx.org/) の [algorithm](https://geocompx.github.io/geocompkg/articles/algorithm.html) vignetteに示されている。)。 @@ -173,8 +178,8 @@ source("code/11-centroid-setup.R") ```{r 10-algorithms-4} # ポリゴンを表現する行列を作成 -x_coords = c(10, 0, 0, 12, 20, 10) -y_coords = c(0, 0, 10, 20, 15, 0) +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) ``` @@ -276,7 +281,7 @@ source("code/11-centroid-alg.R") アルゴリズム\index{algorithm}開発は大変な作業である。 このことは、base R\index{R} を使った重心\index{centroid}アルゴリズム\index{algorithm}の開発に費やした作業量から明らかである。このアルゴリズムは、実世界での応用が限られている問題に対する一つの、むしろ非効率的なアプローチに過ぎない。というのも、実際には凸ポリゴンは珍しいからである。 -この経験は、GEOS\index{GEOS} (`sf::st_centroid()` の基盤となっている) や CGAL\index{CGAL} (計算幾何学アルゴリズムライブラリ) など、高速に動作し、かつ幅広い入力ジオメトリタイプに対応する低レベル地理ライブラリの理解につながるはずである。 +この経験は、GEOS\index{GEOS} や CGAL\index{CGAL} (計算幾何学アルゴリズムライブラリ, Computational Geometry Algorithms Library) など、高速に動作し、かつ幅広い入力ジオメトリタイプに対応する低レベル地理ライブラリの理解につながるはずである。 このようなライブラリのオープンソース化の大きな利点は、そのソースコード\index{source code}が、研究、理解、(技術と自信があれば) 改変のために容易に利用できることである。^[ 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 で見ることができる。