Rでポリゴンを融合する（準備編） | R, GIS
この記事は、Rでポリゴンを融合する為の準備に関する記事です。
記事の概要です。
今回と次回で、茨城県を５つの地域に分けた地図を作成します。
５つの地域は「県北・県央・鹿行（ろっこう）・県南・県西」です。
この記事では、茨城県のシェープファイルに地域区分の属性情報を付与します。
シェープファイルに地域区分の属性情報を付与する部分です。
shp <- mutate(shp, N03_002 = case_when(N03_004 %in% kenhoku ~ "県北",
N03_004 %in% kenou ~ "県央",
N03_004 %in% rokkou ~ "鹿行",
N03_004 %in% kennan ~ "県南",
N03_004 %in% kensei ~ "県西",
TRUE ~ NA_character_ ))kenhokuやkenouなどは、県北や県央に属する市町村名のベクトルです。
こちらはソースコード全体です。
## データ読み込み: 茨城県のシェープファイル
library(sf)
shp <- st_read("./data/N03-20150101_08_GML", options = "ENCODING=SHIFT-JIS",
stringsAsFactors = FALSE)
shp # 確認：シェープファイルの中身
## N03_002カラムに各市町村(N03_004)が属する地域を格納
kenhoku <- c("日立市", "高萩市", "北茨城市", "常陸太田市", "常陸大宮市", "大子町") # 県北
kenou <- c("水戸市", "笠間市", "小美玉市", "茨城町", "大洗町", "城里町",
"ひたちなか市", "那珂市", "東海村") # 県央
rokkou <- c("鹿嶋市", " 神栖市", " 潮来市", " 行方市", " 鉾田市") # 鹿行
kennan <- c("土浦市", "石岡市", "龍ケ崎市", "取手市", "牛久市", "かすみがうら市",
"阿見町", "つくば市", "守谷市", "つくばみらい市", "稲敷市", "美浦村",
"河内町", "利根町") # 県南
kensei <- c("古河市", "常総市", "坂東市", "五霞町", "境町", "結城市", "下妻市",
"筑西市", "桜川市", "八千代町") # 県西
library(dplyr)
shp <- mutate(shp, N03_002 = case_when(N03_004 %in% kenhoku ~ "県北",
N03_004 %in% kenou ~ "県央",
N03_004 %in% rokkou ~ "鹿行",
N03_004 %in% kennan ~ "県南",
N03_004 %in% kensei ~ "県西",
TRUE ~ NA_character_ ))
shp # 確認：N03_002カラムの内容この記事で使用したRとパッケージのバージョンを記載しておきます。
R 3.6.1
sf 0.7-6
dplyr 0.8.3
茨城県の地域区分
茨城県は、県北・県央・鹿行・県南・県西の５つの地域に区分されます。
具体的には、こちらの記事をご確認ください。
シェープファイルに地域区分情報を付与
まず、茨城県のシェープファイルを読み込みます。
シェープファイルが無い場合は、こちらを参考にダウンロードして下さい。
library(sf)
shp <- st_read("./data/N03-20150101_08_GML", options = "ENCODING=SHIFT-JIS")読み込んだ内容を確認しましょう。
shp # 確認：シェープファイルの中身
Simple feature collection with 162 features and 5 fields
geometry type: POLYGON
dimension: XY
bbox: xmin: 139.6877 ymin: 35.73914 xmax: 140.852 ymax: 36.9453
epsg (SRID): NA
proj4string: +proj=longlat +ellps=GRS80 +no_defs
First 10 features:
N03_001 N03_002 N03_003 N03_004 N03_007 geometry
1 茨城県 2 水戸市 08201 POLYGON ((140.3894 36.46392...
2 茨城県 1 日立市 08202 POLYGON ((140.6364 36.47452...
3 茨城県 1 日立市 08202 POLYGON ((140.6301 36.48646...
4 茨城県 1 日立市 08202 POLYGON ((140.6305 36.49491...
5 茨城県 1 日立市 08202 POLYGON ((140.6388 36.49908...
6 茨城県 1 日立市 08202 POLYGON ((140.6353 36.50807...
7 茨城県 1 日立市 08202 POLYGON ((140.6463 36.5421,...
8 茨城県 1 日立市 08202 POLYGON ((140.6467 36.54221...
9 茨城県 1 日立市 08202 POLYGON ((140.6476 36.54248...
10 茨城県 1 日立市 08202 POLYGON ((140.6506 36.54408...特に問題は無いようです。
では、QGISの記事に合わせ、N03_002カラムに各市町村が属する地域を格納しましょう。
プログラムを見やすくするため、市町村を地域毎にまとめたベクトルを作成します。
kenhoku <- c("日立市", "高萩市", "北茨城市", "常陸太田市", "常陸大宮市", "大子町") # 県北
kenou <- c("水戸市", "笠間市", "小美玉市", "茨城町", "大洗町", "城里町",
"ひたちなか市", "那珂市", "東海村") # 県央
rokkou <- c("鹿嶋市", " 神栖市", " 潮来市", " 行方市", " 鉾田市") # 鹿行
kennan <- c("土浦市", "石岡市", "龍ケ崎市", "取手市", "牛久市", "かすみがうら市",
"阿見町", "つくば市", "守谷市", "つくばみらい市", "稲敷市", "美浦村",
"河内町", "利根町") # 県南
kensei <- c("古河市", "常総市", "坂東市", "五霞町", "境町", "結城市", "下妻市",
"筑西市", "桜川市", "八千代町") # 県西作成できたら、N03_002カラムに各市町村（N03_004）が属する地域名を格納します。
library(dplyr)
shp <- mutate(shp, N03_002 = case_when(N03_004 %in% kenhoku ~ "県北",
N03_004 %in% kenou ~ "県央",
N03_004 %in% rokkou ~ "鹿行",
N03_004 %in% kennan ~ "県南",
N03_004 %in% kensei ~ "県西",
TRUE ~ NA_character_ ))case_whenのところでは、市町村（N03_004）一つ一つに対して、その市町村が属する地域名を返しています。
市町村がどの地域にも属していない場合、地域名は欠損値（NA_character_）となります。
簡単に確認しましょう。
shp # 確認：N03_002カラムの内容
Simple feature collection with 162 features and 5 fields
geometry type: POLYGON
dimension: XY
bbox: xmin: 139.6877 ymin: 35.73914 xmax: 140.852 ymax: 36.9453
epsg (SRID): NA
proj4string: +proj=longlat +ellps=GRS80 +no_defs
First 10 features:
N03_001 N03_002 N03_003 N03_004 N03_007 geometry
1 茨城県 県央 水戸市 08201 POLYGON ((140.3894 36.46392...
2 茨城県 県北 日立市 08202 POLYGON ((140.6364 36.47452...
3 茨城県 県北 日立市 08202 POLYGON ((140.6301 36.48646...
4 茨城県 県北 日立市 08202 POLYGON ((140.6305 36.49491...
5 茨城県 県北 日立市 08202 POLYGON ((140.6388 36.49908...
6 茨城県 県北 日立市 08202 POLYGON ((140.6353 36.50807...
7 茨城県 県北 日立市 08202 POLYGON ((140.6463 36.5421,...
8 茨城県 県北 日立市 08202 POLYGON ((140.6467 36.54221...
9 茨城県 県北 日立市 08202 POLYGON ((140.6476 36.54248...
10 茨城県 県北 日立市 08202 POLYGON ((140.6506 36.54408...特に問題は無いようです。
本当はもっと細かく確認すべきだと思います。
ただ、もし問題があれば、次回行う描画時に発見できるはず。
なので、今回はこれくらいで。
おわりに
お疲れさまでした。
次回は実践編です！