メッシュデータについて
GISデータでメッシュ単位でデータを公開している場合があるが、 メッシュデータにはいくつかの決まり事があるらしい。
詳しくは総務省統計局、パスコのHPやArcGIS Tipsがわかりやすい。
代表的なものとして第1次~第3次メッシュがあり、以下のようになっている。
名前 | メッシュ単位 | メッシュコード |
---|---|---|
第一次メッシュ | 約80km四方 | 4桁 |
第二次メッシュ | 約10km四方 | 6桁 |
第三次メッシュ | 約1km四方 | 8桁 |
第一次メッシュは1度ごとの経線と2/3度ごとの緯線によって全国を分割して作られたものだそうである。それらを8等分したのが第2次メッシュ、さらに10等分して作られたのが第3次メッシュだそうである。また、メッシュ単位に 約とついているとおり、地域(緯度経度)によって多少メッシュの大きさが異なるそうである。
メッシュコードは緯度経度をもとに算出されているそうであるが、 ここでは割愛する。
市町村別メッシュコードは総務省統計局で公開している。
Rから国土数値情報をWebAPIを利用するkokudosuuchiパッケージが公開されている。
これらを使えば容易にデータが取得できそうである。 ただ、市町村別メッシュコードの更新年度が2010年なので少し注意が必要かもしれない。
mapviewパッケージによるインタラクティブな作図
leaflet::leaflet()
で基本図の設定を行ったうえで、
mapview::garnishMap()
関数を用いることにより様々な図を
重ねていくことができる。
最後にaddLayersControl()
により、レイヤーの構成を調整することができる。
library(mapview) library(sf) nc <- st_read(system.file("shape/nc.shp", package = "sf")) nc1 <- nc[1:50, ] nc2 <- nc[51:100, ]
m_base <- leaflet((options =leafletOptions(crs = leafletCRS(code="EPSG:4267"))) %>% addTiles(group = "base") %>% addProviderTiles(provider = "Stamen.TonerHybrid", group = "base_steman") m <- garnishMap(m_base, addPolygons, data = as(nc2, "Spatial"), weight = 2, group = "nc2") m <- garnishMap(m, addPolygons, data = as(nc1, "Spatial"), weight = 2, group = "nc2") m <- addLayersControl(m, baseGroups = c("base", "base_steman"), overlayGroups = c("nc1", "nc2")) m
ポリゴン内のラインデータの抽出(切取り)
何回も調べたりしているため、メモ。
ポリゴン内のライン長さを出したいときがある。 その際はラインをポリゴンでインターセクトする必要がある。
関数はrgeos::gIntersection()
を利用する。
sp::over()
関数で行うと意図と異なる結果となるため注意。
(適切な方法で行えば実行できるだろうが。。。)
library(stplanr) library(magrittr) data(routes_fast) data(zones)
- よい例
plot(zones[3,], col = "yellow") plot(routes_fast, col ="grey90", add = T, lwd = 5) rgeos::gIntersection(spgeom1 = routes_fast, spgeom2 = zones[3,]) %>% plot(add = T, col = "blue")
- ダメな例
plot(zones[3,], col = "yellow") plot(routes_fast, col ="grey90", add = T, lwd = 5) tmp <- sp::over(x = routes_fast, y = zones[3, ]) plot(routes_fast[!is.na(tmp$geo_code), ], add = T, col = "red", lwd = 1)
sf classデータの操作及びプロット
sf
パッケージ(クラス)は空間データを簡単に操作することができ非常に便利である。
sp
パッケージ(クラス)はそれらを扱う多様なパッケージが存在するが、
自分自身で多少の編集を行うのは少し面倒であった。
Simple Features for Rだけあってデータ操作は簡単に行うことができる。
dplyr
関係のものもサポートされており、
ggplot2パッケージではこれらのクラスを扱うgeom_sf()
関数も開発されており、
今後も開発が活発に行われることが想像される。
簡単な例題を残しておく
- mutate
nc <- sf::st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
nc %>% mutate(A2 = AREA/ max(AREA)) %>% ggplot(data = .)+ geom_sf(aes(fill = AREA))
- filter
nc %>% filter(AREA > 0.1) %>% ggplot(data = .) + geom_sf(aes(fill = AREA))
- color gradient
nc %>% ggplot()+ geom_sf(aes(fill = AREA))+ theme_bw() + scale_fill_distiller(palette = "YlOrRd", direction = 1)
- annotate
ggplot(nc) + geom_sf() + annotate("point", x = -80, y = 35, colour = "red", size = 4) ggplot(nc) + geom_sf() + annotate("point", x = -80, y = 35, colour = "red", size = 4) + geom_text(aes(x= -80, y = 35, label = "annotate"), hjust = "inward", vjust = "inward")
ggpairs()できれいなペアプロット
pairs()をggplot2風に書けるらしい。
library(tidyverse) library(GGally) iris %>% GGally::ggpairs(data = ., aes(color = Species))
rmarkdwonでの参考文献
ヘッダーは下記のようにする。
--- output: html_document: bibliography: bibdb.bib csl: japan-journal-of-industrial-and-applied-mathematics.csl ---
文章中の参考としたところには[@****]
と表示してやればよい。
参考文献のスタイルを指定するcslファイルはZetoro style repositoryから自分に合うものを探す。
そして最終行を下記のようにする。(書かなくても表示はされる)
# References
参考したサイト
sfクラスからSpatial*****DataFrameクラスへの変換
Rにおけるジオメトリデータの操作は今後sf
パッケージが主流になるらしい。
http://notchained.hatenablog.com/entry/2017/01/06/213333
開発が活発なためせっかく整備したコードが今後も使えるという保証はないため, 利用に少しためらいもあるが,最近少し使っている。
sfクラスからSpatial***DataFrameクラスに変換する際のメモである。
通常は以下のようにas(., "Spatial")
で従来のSpatial****DataFrameに変換することができる。
> nc <- st_read(system.file("shape/nc.shp", package="sf"), quiet = TRUE) > > class(nc) [1] "sf" "data.frame" > methods::as(nc, "Spatial") %>% class [1] "SpatialPolygonsDataFrame" attr(,"package") [1] "sp"
しかし,tbl_df
クラスも持つ場合にはエラーが生じてしまう。
その場合は以下のようにして,クラスをsf
, data.frame
にすると解決できる。
nc %>% `class<-`(c("sf", "data.frame")) %>% methods::as(., "Spatial")