Dimension Planet Adventure 最終章 最終話『栄光なる未来』

魔神ベベカル!!もうこの世界を貴様の好きにはさせん!!ここで俺が倒す!!!

【追記あり】sfパッケージでシェープファイルを読み込んでmapviewパッケージで可視化するまで

お久しぶりです。 位置情報データ分析をやっていくことになりました。 素人なので備忘録を書いていきます。 GISを使えませんが、Rなら少々できますので、限界が来るまでRでチャレンジしてみます。 周りにRで位置情報分析をする人がいない環境なので、ツッコミを欲しています。 何卒よろしくお願いいたします。


今日の環境はMacです。

まず、位置情報データを入手し、可視化することにチャレンジしてみます。 今回のお題は地域メッシュです。とにかくシェープファイルを入手することからスタートです。 位置情報系のデータは、国土数値情報ダウンロードサービスに多くのデータがあるので、今回もこちらを使用します。下記URLです。

国土数値情報 土地利用3次メッシュデータの詳細

全国分をダウンロードすると容量が大変なことになりそうなので、東京都のデータだけにしておきます。 zipを展開して任意の位置に配置します。

シェープファイルの取り扱いには sf パッケージを使用するのが 2017年現在ナウいようです。詳しくは いつもお世話になっている @yutannihilation さんの下記ブログをご参照ください。

notchained.hatenablog.com

では、まずはシェープファイルを読み込みます。

> library(sf)
> d1 <- sf::read_sf("./Downloads/L03-a-14_5339-tky_GML/L03-a-14_5339-tky.shp")
 make.names(vnames, unique = TRUE) でエラー: 
   '<83><81><83>b<83>V<83><85>' に不正なマルチバイト文字があります 
 wrapup 中にエラーが起こりました:  '<83><81><83>b<83>V<83><85>' に不正なマルチバイト文字があります 

怒られました。

sfパッケージにおけるマルチバイト文字を含むシェープファイルを読み込む方法が結局わからなかった*1ので、今回は従来の方法のひとつである rgdal パッケージで読み込んでから sf パッケージで使えるクラスに変換する方針でやっていきましょう。

> library(rgdal)
> d1 <- rgdal::readOGR("./Downloads/L03-a-14_5339-tky_GML/L03-a-14_5339-tky.shp", encoding = "Shift_JIS") %>%
+   sf::st_as_sf() # Convert to an sf object
OGR data source with driver: ESRI Shapefile 
Source: "./Downloads/L03-a-14_5339-tky_GML/L03-a-14_5339-tky.shp", layer: "L03-a-14_5339-tky"
with 6300 features
It has 13 fields
Integer64 fields read as strings:  田 他農用地 森林 荒地 建物用地 道路 鉄道 他用地 河川湖沼 海浜 海水域 ゴルフ場 
> d1 %>% head
Simple feature collection with 6 features and 13 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 139 ymin: 35.33333 xmax: 139.075 ymax: 35.34167
epsg (SRID):    NA
proj4string:    +proj=longlat +ellps=bessel +no_defs
  メッシュ    田 他農用地    森林   荒地 建物用地 道路 鉄道 他用地 河川湖沼 海浜 海水域 ゴルフ場
0 53390000     0        0 1050436      0        0    0    0      0        0    0      0        0
1 53390001     0        0  976895      0        0    0    0      0        0    0      0    73530
2 53390002     0        0  808820 241595        0    0    0      0        0    0      0        0
3 53390003     0        0 1050404      0        0    0    0      0        0    0      0        0
4 53390004     0        0 1039890      0    10504    0    0      0        0    0      0        0
5 53390005 42015   136550  829804  21008        0    0    0  21007        0    0      0        0
                        geometry
0 POLYGON((139 35.33333333333...
1 POLYGON((139.0125 35.333333...
2 POLYGON((139.025 35.3333333...
3 POLYGON((139.0375 35.333333...
4 POLYGON((139.05 35.33333333...
5 POLYGON((139.0625 35.333333...

いい感じです。


【追記01: 2017/06/27 20時】

コメントにてsfパッケージで読み込む方法を教えていただきました。 Macでもマルチバイト文字がきちんと表示されていました。 id:aaaazzzz036 さん、ありがとうございます。

d1 <- sf::st_read("./Downloads/L03-a-14_5339-tky_GML/L03-a-14_5339-tky.shp", options=c("ENCODING=CP932"))

【追記02: 2017/06/29 11時】

sf::read_sf 関数でも同じオプションでOKです。

d1 <- sf::read_sf("./Downloads/L03-a-14_5339-tky_GML/L03-a-14_5339-tky.shp", options=c("ENCODING=CP932"))

では、可視化してみましょう。

位置情報データの可視化には mapview パッケージを使ってみます。 今後の目標としては sfmapview をある程度使えるようになりたい というものがあります。

今回はお試しなので、新宿近辺の4メッシュほどを可視化してみましょう。

> library(mapview)
> library(dplyr)
> d1 %>% 
+   dplyr::filter(メッシュ %in% c(53394525, 53394526, 53394535, 53394536)) %>% 
+   mapview::mapview()

f:id:ksmzn:20170627145320p:plain

こちらもいい感じです。左のボタンからOpenStreetMapに切り替えると見やすいですね。

今の例のように、sf パッケージも mapview パッケージも、tidyverse系ライブラリと組み合わせて使えるのが便利すぎます。

今後も少しずつ sf パッケージと mapview パッケージをやっていきます。

SFが読みたい! 2017年版

SFが読みたい! 2017年版

*1:ご存知のお方は、教えていただけると幸いです…