建設総合統計の可視化
国土交通省では、国内の建設活動を出来高ベースで把握することを目的に建設総合統計を毎月発表している。 気が向いたので、ちょっと可視化してみた。 使用するデータは2017/9/19に発表された2010年4月から2017年7月分までのデータを利用する。
都道府県別の毎年4月のデータの推移から下記のことがわかる
東北三県に着目し出来高および東北3県に占める各県の出来高割合の推移をみると、 東日本大震災の復興のため3県とも建設需要が増加している。
特に宮城県で増加しており以下のようなことが考えられる。
宮城県の震災前 約500億円 –> 2000億円以上 と約4倍の建設特需はすごいと思う。
建設動向としては、 復興関連の需要は減少し、東京オリンピックに向けた需要が増加中といえる。
毎年12月の出来高の増加率の推移をみると、2013年はほとんどの都道府県で出来高が増加しており、
安倍政権になって建設需要が急増したといえる。
都道府県別の月別の出来高増加率の推移をみると、年間を通して発注・完成時期がある程度平準化している地域とそうでない地域があるのがわかる。 国交省では施工時期の平準化を進めているようであるが、これを見ればその進捗状況が一目でわかるといえる。
青森県は、雪の影響で一時期に発注が集中すること等があるのか極端に増減がある。
sf::read_sf()とrgdal::readOGR()とtmaptools::read_shp()の速度比較
Rでシェープファイルデータを読込むための関数でどれが一番速度が速いか比較する。
library(sf) library(tmaptools) library(rgdal) library(microbenchmark)
rfile <- system.file("shape/nc.shp", package="sf") microbenchmark(sf::read_sf(rfile), sf::st_read(rfile, quiet=T), rgdal::readOGR(rfile, verbose=F), tmaptools::read_shape(rfile), times=10)
Unit: milliseconds
expr min lq mean
sf::read_sf(rfile) 8.671537 8.976245 14.69445
sf::st_read(rfile, quiet = T) 6.392833 6.685836 15.94281
rgdal::readOGR(rfile, verbose = F) 567.662413 597.656352 691.19212
tmaptools::read_shape(rfile) 583.133058 596.592328 654.66730
median uq max neval cld
10.928342 14.948451 36.55517 10 a
6.852349 7.603736 54.89468 10 a
623.355303 678.476535 1090.38659 10 b
626.761277 651.405453 836.15833 10 b
ということでsfパッケージ一択でした。
microbenchmark(read_sf(rfile), st_read(rfile, quiet=T, stringsAsFactors = F), times=20)
Unit: milliseconds
expr min lq
read_sf(rfile) 8.032669 8.185589
st_read(rfile, quiet = T, stringsAsFactors = F) 5.793233 5.904619
mean median uq max neval cld
8.703159 8.319064 9.221673 11.538136 20 b
6.202890 6.021859 6.422662 8.131218 20 a
毎回引数指定めんどくさいかもしれないが、大規模なデータを読み込む場合はst_read()がよいみたい。
もっと早く読み込める関数があったら教えてください。
ついでにggplot()とtmap()の速度について比較してみる。
library(ggplot2) library(tmap)
nc <- st_read(rfile, quiet=T, stringsAsFactors = F)
microbenchmark(ggplot(nc)+geom_sf(aes(fill=CNTY_)),
tm_shape(nc)+tm_polygons(col="CNTY_"))
Unit: microseconds
expr min lq mean
ggplot(nc) + geom_sf(aes(fill = CNTY_)) 2419.920 2465.9850 2947.4016
tm_shape(nc) + tm_polygons(col = "CNTY_") 114.408 128.1895 214.9879
median uq max neval cld
2528.097 2605.502 24483.139 100 b
157.263 175.387 6131.546 100 a
ggplot2は使い慣れているから便利だけど、やはり描画速度に難があるみたい。 特に空間情報データの作図の際にはその傾向が強いと思う。
mapview::mapview()の引数オプション
mapview::mapview()で利用する引数に関するメモ。
library(mapview) library(sf) nc <- read_sf(system.file("shape/nc.shp", package="sf")) mapview(nc, zcol="AREA", # 色分けするデータ（列） #color="blue", # ラインの色 #col.regions="red", # ポリゴンの色 cex=12, # ポイントデータの大きさ(default:6) na.color="white", # NAの色 lwd = 2, # ラインの太さ map.types="OpenStreetMap.DE", # ベースマップ layer.name="test", # レイヤーの名前(右下に表示されるレイヤ名) #alpha=0.5, # ラインの透過率（0:透過、1:透過なし) #alpha.regions=0.5, # ポリゴンの透過率 #na.alpha=0.4 # NA値の透過率 popup=popupTable(nc[,"AREA"]), # 対象魏小目鳥をクリックした際に表示する項目 label=nc$AREA %>% as.character(), # マウスオーバした際に表示されるラベル legend=T # 凡例を表示するかどうか )
sfオブジェクトをggplot()する場合はgeometryカラムの名前に注意
geom_sf()関数を利用するにあたり
geometryカラムがgeometryではなくgeomの場合にはエラーが出た。
解決方法のメモ。
library(sf) library(ggplot2) library(spData)
data(world)
head(world)
Simple feature collection with 6 features and 10 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -73.41544 ymin: -55.25 xmax: 75.15803 ymax: 42.68825
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
iso_a2 name_long continent region_un subregion
1 AF Afghanistan Asia Asia Southern Asia
2 AO Angola Africa Africa Middle Africa
3 AL Albania Europe Europe Southern Europe
4 AE United Arab Emirates Asia Asia Western Asia
5 AR Argentina South America Americas South America
6 AM Armenia Asia Asia Western Asia
type area_km2 pop lifeExp gdpPercap
1 Sovereign country 652270.07 31627506 60.37446 1844.022
2 Sovereign country 1245463.75 24227524 52.26688 6955.960
3 Sovereign country 29694.80 2893654 77.83046 10698.525
4 Sovereign country 79880.74 9086139 77.36817 63830.700
5 Sovereign country 2784468.59 42980026 76.15861 18872.967
6 Sovereign country 28656.60 3006154 74.67571 7706.133
geom
1 MULTIPOLYGON(((61.210817091...
2 MULTIPOLYGON(((16.326528354...
3 MULTIPOLYGON(((20.590247430...
4 MULTIPOLYGON(((51.579518670...
5 MULTIPOLYGON(((-65.5 -55.2,...
6 MULTIPOLYGON(((43.582745802...
ggplot(world)+geom_sf() Error in FUN(X[[i]], ...) : object 'geometry' not found
オブジェクトのクラスはdata.frame, sfとなっているにもかかわらず、
先ほどのコードはエラーとなる。
原因はgeometryカラムの名前がgeomとなっていることである。
githubのggplot2にissueとして挙げられており、
geom_sf()にデータを明示的に与える必要があるようである。
もしくは下記のようにaes()にgeometry=geomとすることで解決できる。
class(world)
[1] "sf" "data.frame"
names(world)
[1] "iso_a2" "name_long" "continent" "region_un" "subregion" [6] "type" "area_km2" "pop" "lifeExp" "gdpPercap" [11] "geom"
ggplot(world)+geom_sf(aes(geometry=geom))
なお、geometryカラムの名前がgeometryの場合は問題なくggplot()により作図できる。
nc <- read_sf(system.file("shape/nc.shp", package="sf")) names(nc)
[1] "AREA" "PERIMETER" "CNTY_" "CNTY_ID" "NAME" [6] "FIPS" "FIPSNO" "CRESS_ID" "BIR74" "SID74" [11] "NWBIR74" "BIR79" "SID79" "NWBIR79" "geometry"
ggplot(nc)+geom_sf()
facet_wrap()利用時にすべてのデータを表示させておく
ggplot2においてfacet_wrap()関数は便利であり利用頻度が高いが、
当該分類データ以外のデータを背景として表示させ、
分布状況をより明快にしたい場合がある。
library(ggplot2) library(dplyr)
mtcars %>% ggplot(aes(disp, mpg)) + geom_point()+ facet_wrap(~cyl)
このようなときには引数data=にfacet_wrap()で分類する項目を NULL としたデータを与えることで
グラフを作成することが可能となる。
mtcars %>% ggplot(aes(disp, mpg)) + geom_point(data = . %>% mutate(cyl=NULL), color="darkgrey")+ geom_point()+ facet_wrap(~cyl)
特定の値をNAに置き換える（またはNAを特定の値に置き換える）
国土数値情報からダウンロードしたメッシュデータの中にはunknownという値が含まれており、
各フィールドはキャラクターとして認識され、意図している数値型として認識されていない事例に遭遇した。
皆さんどうしてるんですかね？
とりあえずunknownという値を置き換えるべくネット検索したところ、
gdataパッケージに便利なものがあったのでメモ。
詳細はpackage Vignetteに詳しい
library(gdata)
xNum <- c(0, 6, 0, 7, 8, 9, NA) xNum
[1] 0 6 0 7 8 9 NA
デフォルトではNAがunknownとして認識される。
isUnknown(x = xNum)
[1] FALSE FALSE FALSE FALSE FALSE FALSE TRUE
0をunknownとする場合は以下のようにする。
isUnknown(x = xNum, unknown = 0)
[1] TRUE FALSE TRUE FALSE FALSE FALSE FALSE
複数のunknownとして指定する場合には以下のようにする。
isUnknown(x = xNum, unknown = c(0, NA))
[1] TRUE FALSE TRUE FALSE FALSE FALSE TRUE
NAを置き換える場合にはNAToUnknown()関数を利用する。
NAToUnknown(x = xNum, unknown = 999)
[1] 0 6 0 7 8 9 999
unknownをNAに置き換える場合はunknownToNA()関数を利用する
unknownToNA(x = xNum, unknown=6)
[1] 0 NA 0 7 8 9 NA
dNum <- data.frame( a = c(1, 3, "unknown"), b = c(4, 5, NA) )
データフレームに対しても利用できる。
isUnknown(x = dNum, unknown = NA)
a b 1 FALSE FALSE 2 FALSE FALSE 3 FALSE TRUE
NAToUnknown(x = dNum, unknown = 999)
a b 1 1 4 2 3 5 3 unknown 999
unknownToNA(dNum, unknown = c("unknown"))
a b 1 1 4 2 3 5 3 <NA> NA
read_sf()とst_read()
st_read()とread_sf()
最近read_sf()という関数を知った。
st_read()との違いに関してメモ。
read_sf()の関数について確認すると以下のようになっており、
quiet=TRUE, stringsAsFactors=FALSEでデータを読み込んでいる。
library(sf) library(dplyr)
read_sf
function (..., quiet = TRUE, stringsAsFactors = FALSE)
{
if (!requireNamespace("tibble", quietly = TRUE))
stop("package tibble not available: install first?")
st_as_sf(tibble::as_tibble(as.data.frame(st_read(..., quiet = quiet,
stringsAsFactors = stringsAsFactors))))
}
<environment: namespace:sf>
nc <- read_sf(system.file("shape/nc.shp", package="sf")) glimpse(nc)
Observations: 100 Variables: 15 $ AREA <dbl> 0.114, 0.061, 0.143, 0.070, 0.153, 0.097, 0.062, 0.0... $ PERIMETER <dbl> 1.442, 1.231, 1.630, 2.968, 2.206, 1.670, 1.547, 1.2... $ CNTY_ <dbl> 1825, 1827, 1828, 1831, 1832, 1833, 1834, 1835, 1836... $ CNTY_ID <dbl> 1825, 1827, 1828, 1831, 1832, 1833, 1834, 1835, 1836... $ NAME <chr> "Ashe", "Alleghany", "Surry", "Currituck", "Northamp... $ FIPS <chr> "37009", "37005", "37171", "37053", "37131", "37091"... $ FIPSNO <dbl> 37009, 37005, 37171, 37053, 37131, 37091, 37029, 370... $ CRESS_ID <int> 5, 3, 86, 27, 66, 46, 15, 37, 93, 85, 17, 79, 39, 73... $ BIR74 <dbl> 1091, 487, 3188, 508, 1421, 1452, 286, 420, 968, 161... $ SID74 <dbl> 1, 0, 5, 1, 9, 7, 0, 0, 4, 1, 2, 16, 4, 4, 4, 18, 3,... $ NWBIR74 <dbl> 10, 10, 208, 123, 1066, 954, 115, 254, 748, 160, 550... $ BIR79 <dbl> 1364, 542, 3616, 830, 1606, 1838, 350, 594, 1190, 20... $ SID79 <dbl> 0, 3, 6, 2, 3, 5, 2, 2, 2, 5, 2, 5, 4, 4, 6, 17, 4, ... $ NWBIR79 <dbl> 19, 12, 260, 145, 1197, 1237, 139, 371, 844, 176, 59... $ geometry <simple_feature> MULTIPOLYGON(((-81.47275543..., MULTIPOLY...
ということで、今後は基本的にはこの関数を利用したほうが便利そうである。
引数のdsn=が日本語を含むパスだとうまくいかないのですが、
解決方法をご存知の方ぜひコメント等教授お願いします。
