Rで調査地点を示した地図をつくる
まえがき
自分用のメモとして残します。RもGISも初心者でほぼ独学なので、参考にされる方は自己責任でお願いします。
調査地点を示した地図を作製する方法は多岐にわたります。筑波大生であれば大学がライセンスを持っているためArcGISを利用することができますが、私は完成した図のデザインや操作性があまり好きではありません。そこで、Rを用いて作製してみることにしました。QGISを用いて作製する場合は下の記事が大変参考になると思います。
なお、私はWindows10でRStudioを使用しています。このあたりの情報は何を提供するべきなのかわかっていません。
目的
調査地点を点で示した地図を作製することです。私は水生生物の調査を行うことが多いので、河川と湖沼も地図で示します。しかし、河川を面として表示する方法を知らないので、川幅が表示できなくとも違和感がない、都道府県以上のスケールで地図を作製することをオススメします。河川データが必要ない場合はより小さいスケールで違和感のない地図が作成できます。
今回は茨城県で適当に選んだいくつかの調査地点を表示した地図を作成します。時間とお金が無限にあれば茨城県を調査しつくしたいです。
準備
作製に用いるデータを準備していきます。まずは調査地点の座標を取得します。Googleマップで調査地点にピンを指して座標を取得しました。以下の画像のようにExcelに打ち込み、csv形式で保存しておきます。Excelは勝手に四捨五入するので情報を減らされないように注意。ファイル名は「site」です。
続いて地図の作製に用いるシェープファイルをダウンロードします。シェープファイルとは、地図と色々なデータが紐づけされたファイルのことです。
『政府統計の総合窓口(e-Stat)』から茨城県の境界データを選びます。
地図>境界データダウンロード>小地域>国勢調査>2020年>小地域(町丁・字等)(JGD2011)>世界測地系緯度経度・Shapefile>茨城県>茨城県全域
これでダウンロードできました。ダウンロードしたフォルダを解凍すると4つのファイルが含まれていることがわかりますが、この4つが一緒のフォルダに入っていることが重要らしいので、SHPファイルのみを移動させようとしてはいけません。フォルダ名が長いので「茨城県境界」に変更してわかりやすい場所に移動させておきます。
続いて、『国土数値情報ダウンロードサイト』から河川と湖沼のシェープファイルをダウンロードします。
国土数値情報>水域>河川(ライン)(ポイント)>茨城
国土数値情報>水域>湖沼(ポリゴン)>全国
同様に解凍して名前を変更して移動させておきます。「茨城県河川」、「湖沼」というフォルダ名にしました。湖沼の都道府県ごとのデータをください。
地図の作製
ここからはRStudio上での作業です。地図の作製に用いるパッケージは「dplyr」、「ggplot2」、「sf」です。これらがインストールされていない場合はインストールし、呼び出しておきます。
install.packages("dylyr")
install.packages("ggplot2")
install.packages("sf")
library(dplyr)
library(ggplot2)
library(sf)
シェープファイルとCSVファイルを読み込みます。
ibaraki <- read_sf("R/茨城調査/茨城県境界/r2ka08.shp")
river <- read_sf("R/茨城調査/茨城県河川/W05-08_08-g_Stream.shp")
lake <- read_sf("R/茨城調査/湖沼/W09-05-g_Lake.shp")
site <- read.csv("R/茨城調査/site.csv")
ちらっと確認。この作業が本当に必要なのかは知りません。
glimpse(ibaraki)
glimpse(river)
glimpse(lake)
glimpse(site)
茨城県の境界データをそのまま表示すると市区町村よりも細かい境界が出現してしまうので、県単位で結合(Dissolve)します。groub_by(CITY_NAME)とすれば市区町村単位で結合することもできますが、変な線が入ってしまいます。sf_use_s2(FALSE)としないと結合できませんでしたが、仕組みはわかりません。
sf_use_s2(FALSE)
ibaraki2 <- ibaraki %>% group_by() %>% summarise()
glimpse(ibaraki2)
座標系の変換および設定を行います。これをやらないと調査地点が上手く地図上に表示できません。ここがわからず私は2時間くらい詰まりました。
ibaraki3 <- st_transform(ibaraki2, 4326)
river2 <- st_set_crs(river, 4326)
lake2 <- st_set_crs(lake, 4326)
EPSGコードについては以下を参照。おそらく4326しか使わないと思います。
ggplotで気合で作図します。繰り返し表示して微調整します。
ggplot() +
#基本的なレイヤーの重ね合わせ
#colorで境界線の色を指定しfillで塗りつぶしの色を指定
geom_sf(data = ibaraki3, color = NA, fill = "gray") +
geom_sf(data = river2, color = "dimgray", linewidth = 0.3) +
geom_sf(data = lake2, color = NA, fill = "dimgray") +
#aesでxとyを指定して調査地点の図示
geom_point(data = site, aes(x = lng, y = lat), size = 3) +
#表示範囲の指定
scale_x_continuous(limits = c(139.6,141)) +
scale_y_continuous(limits = c(35.7,37)) +
#テーマの指定
theme_test() +
#軸ラベルを無くして数値の大きさを調節
theme(axis.title = element_blank(), text = element_text(size = 12))
RStudioのPlots画面では文字化けしていましたが、PDFにエクスポートすると解消されました。
最後に残っていた茨城県外の池を画像編集ソフトで消します。地点のラベルも入れます。
できました。とりあえず納得の出来です。茨城県を白く表示することもできますが、この配色が私は好きです。
おわりに
もっと簡単で綺麗に作図する方法がある気がします。ネット上で探しましたが私には理解できないものが多かったので、誰かの参考になれば幸いです。何かご指摘がありましたらコメントでお願いします。
参考資料の後にコード全体を載せています。上で断片的に見せていたコードを全て繋げただけです。作法などは知らないのであしからず。
参考資料
Rをはじめよう生命科学のためのRStudio入門
この本でRの基礎を学びました。Rを使ってシェープファイルから白地図を描く
最初はここを参考にしました。Rを使った地理空間データの可視化と分析
ここで色々学びましたが、8割以上理解できていません。Rで生態学データのGIS
一番詰まった座標系の変換や指定のことに気づけました。Vector operations in R
ディゾルブのやり方を学びました。Loop 0 is not valid: Edge x has duplicate vertex with edge y.
この記事のおかげで茨城県のディゾルブに成功しました。
rm(list=ls())
library(dplyr)
library(ggplot2)
library(sf)
#読み込み
ibaraki <- read_sf("R/茨城調査/茨城県境界/r2ka08.shp")
river <- read_sf("R/茨城調査/茨城県河川/W05-08_08-g_Stream.shp")
lake <- read_sf("R/茨城調査/湖沼/W09-05-g_Lake.shp")
site <- read.csv("R/茨城調査/site.csv")
#確認
glimpse(ibaraki)
glimpse(river)
glimpse(lake)
glimpse(site)
#茨城県のディゾルブ
#group_byで属性ごとにディゾルブできる
sf_use_s2(FALSE)
ibaraki2 <- ibaraki %>% group_by() %>% summarise()
glimpse(ibaraki2)
#座標系の変換および設定
ibaraki3 <- st_transform(ibaraki2, 4326)
river2 <- st_set_crs(river, 4326)
lake2 <- st_set_crs(lake, 4326)
#地図の表示
ggplot() +
#基本的なレイヤーの重ね合わせ
#colorで境界線の色を指定しfillで塗りつぶしの色を指定
geom_sf(data = ibaraki3, color = NA, fill = "gray") +
geom_sf(data = river2, color = "dimgray", linewidth = 0.3) +
geom_sf(data = lake2, color = NA, fill = "dimgray") +
#aesでxとyを指定して調査地点の図示
geom_point(data = site, aes(x = lng, y = lat), size = 3) +
#表示範囲の指定
scale_x_continuous(limits = c(139.6,141)) +
scale_y_continuous(limits = c(35.7,37)) +
#テーマの指定
theme_test() +
#軸ラベルを無くし数値の大きさを調節
theme(axis.title = element_blank(), text = element_text(size = 12))