Tokyo.R #107がジオデータ特集をすると聞いて、GDALを思い出そう思ったのと、折角だし、sfパッケージも勉強して、rgdalとの対比をLTしようともくろんだ。
でも、ちょっと気が変わったので、その準備だけ残しておくことにした。
次回、第107回 #TokyoR は2023年7月15日（土）開催ジオデータ特集です！
— Tokyo.R (@TokyoRCommunity) June 16, 2023
LT発表者募集を開始しました！ジオデータに限らずテーマもレベルも任意です。以下のフォームから奮ってご応募ください！ https://t.co/HAFnJOywhl
なお、発表者以外の方の参加受付は後日開始します。しばらくお待ち下さい。
GISデータ
GISデータは地理的な位置や形状、属性などの情報を持つデータで、ベクター形式やラスター形式などの種類がある。
GISデータは複数のファイルから構成されており、シェープファイル（shpファイル）が図形の情報、shxファイルが図形のインデックス情報、dbfファイルが図形の属性情報をそれぞれ格納している。
shpファイルはESRI社が策定した業界標準フォーマットで、ポイントやラインやポリゴンなどの図形の情報と属性情報を格納したベクター形式データ。
RでGISデータの分析するパッケージ
調べると、RでのGISデータ処理が変わり目にあると知った。
| パッケージ | 特徴 | 利点 | 欠点 |
|---|---|---|---|
| rgdal | GISデータの読み書きや変換を行うGeospatial Data Abstraction Library (GDAL) を使う | ラスターとベクターの地理空間データの取り込み、spオブジェクトの出力 | 2023年10月に廃止される、sfパッケージに移行することになる |
| sf | Simple FeaturesというGISデータの標準仕様に対応 | 地理空間データをデータフレームに変換、tidyverseとの互換性 | sfオブジェクトを扱うので、spオブジェクトを使えない |
| spdep | 空間自己相関や空間回帰などの空間統計分析 | 空間統計の定番パッケージ、モラン散布図や地域指標などの分析 | sfオブジェクトに対応していない、as_Spatial関数で変換する必要がある |
| gstat | 地理空間データの統計モデリング | 点データや面データに対する重回帰分析、kriging法のモデリング | モデリングに必要な前処理や後処理が多い |
| raster | 地表面を格子状に区切ったラスター形式のGISデータを扱う | ラスターデータの処理、他の地理空間データとの統合 | データフレームを扱わない、tidyverseパッケージとの互換性が低い |
| tmap | 地理空間データの可視化 | ggplot2に似た文法を使える、静的、動的に処理できる | カスタマイズには多くのオプションが必要 |
RでGISデータを使ってみる
データをダウンロードする
少し前に流行った日本の都道府県のGISデータにCOVID-19のデータを重ねてみる
日本の都道府県のGISデータ
wd <- "d:/Desktop/analysis_gis_data/" setwd(wd) # 日本の都道府県のGISデータをダウンロード x <- "https://nlftp.mlt.go.jp/ksj/gml/data/N03/N03-20230101_GML.zip" file_name <- sub(".*/", "./", x) download.file(x, file_name) # ダウンロードしたファイルを解凍 unzip(file_name) shp_file <- dir(pattern = "\\.shp$", full.names = TRUE)
NHKがまとめたCOVID-19のデータ
# データをダウンロード x <- "https://www3.nhk.or.jp/n-data/opendata/coronavirus/nhk_news_covid19_prefectures_daily_data.csv" file_name <- sub(".*/", "./", x) download.file(x, file_name) # データを読み込む covid <- read.csv(file_name) covid <- subset(covid, select = c(日付, 都道府県名, 各地の死者数_累計)) # データを表示する head(covid)
$\ast$ JSON形式でも取得できる
# # データをダウンロード # x <- "https://www3.nhk.or.jp/news/special/coronavirus/data/latest-pref-data.json" # file_name <- sub(".*/", "./", x) # download.file(x, file_name) # # データを読み込む # covid <- jsonlite::fromJSON(file_name)
Rのパッケージを使う
rgdalパッケージ
# rgdalパッケージを読み込む library(rgdal) # 解凍したファイルからGISデータを読み込む gis_data <- readOGR(shp_file) # GISデータをプロットする plot(gis_data)
sfパッケージ
# sfパッケージとggplot2パッケージをインストール install.packages(c("sf", "ggplot2")) # sfパッケージとggplot2パッケージを読み込む library(sf) library(ggplot2) # 解凍したファイルからGISデータを読み込む pref <- st_read("N03-20200101_21_GML/N03-21_200101.shp") # GISデータをggplot2でプロットする ggplot(pref) + geom_sf() + theme_minimal()
spdepパッケージ
# spdepパッケージとtmaptoolsパッケージとtidyverseパッケージをインストール install.packages(c("spdep", "tmaptools", "tidyverse")) # spdepパッケージとtmaptoolsパッケージとtidyverseパッケージを読み込む library(spdep) library(tmaptools) library(tidyverse) # 解凍したファイルからGISデータを読み込む pref <- read_shape("N03-20200101_21_GML/N03-21_200101.shp") # GISデータとcovid-19死亡者数と人口密度のデータフレームを結合する pref_covid <- pref %>% left_join(covid, by = c("N03_001" = "pref")) # GISデータをspオブジェクトに変換する pref_covid <- as_Spatial(pref_covid) # 空間的自己相関分析を行うための隣接行列を作成する nb <- poly2nb(pref_covid, queen = TRUE) # 空間的自己相関分析を行うための重み行列を作成する w <- nb2listw(nb, style = "W") # 死亡者数の空間的自己相関分析を行う moran.test(pref_covid$death, w) # モラン散布図をプロットする moran.plot(pref_covid$death, w)
gstatパッケージ
# gstatパッケージとspパッケージとtidyverseパッケージをインストール install.packages(c("gstat", "sp", "tidyverse")) # gstatパッケージとspパッケージとtidyverseパッケージを読み込む library(gstat) library(sp) library(tidyverse) # GISデータとcovid-19死亡者数と人口密度と都道府県庁所在地の緯度と経度のデータフレームを結合する pref_covid <- pref %>% left_join(covid, by = c("N03_001" = "pref")) # GISデータをspオブジェクトに変換する coordinates(pref_covid) <- ~lon + lat # 死亡者数に対して重回帰分析を行う lm1 <- lm(death ~ density + N03_001, data = pref_covid) # 重回帰分析の結果を表示する summary(lm1) # 死亡者数に対してkriging法を行う kr1 <- krige(death ~ density + N03_001, pref_covid, pref_covid) # kriging法の結果を表示する summary(kr1)
rasterパッケージ
# rasterパッケージとrgdalパッケージをインストール install.packages(c("raster", "rgdal")) # rasterパッケージとrgdalパッケージを読み込む library(raster) library(rgdal) # 解凍したファイルからラスターデータを読み込む dem <- raster("FG-GML-6441-36-DEM5A.tif") # ラスターデータをプロットする plot(dem)
tmapパッケージ
# tmapパッケージとsfパッケージとtidyverseパッケージをインストール install.packages(c("tmap", "sf", "tidyverse")) # tmapパッケージとsfパッケージとtidyverseパッケージを読み込む library(tmap) library(sf) library(tidyverse) # 解凍したファイルからGISデータを読み込む pref <- st_read("N03-20200101_21_GML/N03-21_200101.shp") # GISデータとcovid-19死亡者数と人口密度のデータフレームを結合する pref_covid <- pref %>% left_join(covid, by = c("N03_001" = "pref")) # GISデータをtmapでプロットする tm_shape(pref_covid) + tm_fill(col = "death") + tm_borders() + tm_layout(title = "日本の都道府県別のcovid-19死亡者数")
