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=が日本語を含むパスだとうまくいかないのですが、
解決方法をご存知の方ぜひコメント等教授お願いします。
sfパッケージによる空間データの操作
利用パッケージ
library(sf) library(tidyverse) library(viridis) library(jpndistrict)
データの準備
fukuoka <- jpndistrict::spdf_jpn_pref(code=40) class(fukuoka)
[1] "SpatialPolygonsDataFrame" attr(,"package") [1] "sp"
spクラスからsfクラスへの変換
fukuoka_sf <- st_as_sf(fukuoka) class(fukuoka_sf)
[1] "sf" "data.frame"
CRSをセットする。
fukuoka_sf <- fukuoka_sf %>% st_set_crs(4612)
出力（作図）する
- 通常のプロット
fukuoka_sf %>% st_geometry() %>% plot()
- ggplot2を利用したプロット
fukuoka_sf %>% ggplot() + geom_sf()
dplyr関連の操作を行う
- filter()
北九州市のみを選択する
fukuoka_sf %>% filter(city_name_=="北九州市")
Simple feature collection with 23 features and 5 fields
geometry type: POLYGON
dimension: XY
bbox: xmin: 129.9819 ymin: 33.00003 xmax: 131.1906 ymax: 34.24994
epsg (SRID): 4612
proj4string: +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
First 20 features:
pref_name city_name_ city_name city_name_full city_code
1 福岡県 北九州市 門司区 北九州市 門司区 40101
2 福岡県 北九州市 門司区 北九州市 門司区 40101
3 福岡県 北九州市 門司区 北九州市 門司区 40101
4 福岡県 北九州市 若松区 北九州市 若松区 40103
5 福岡県 北九州市 若松区 北九州市 若松区 40103
6 福岡県 北九州市 若松区 北九州市 若松区 40103
7 福岡県 北九州市 若松区 北九州市 若松区 40103
8 福岡県 北九州市 若松区 北九州市 若松区 40103
9 福岡県 北九州市 若松区 北九州市 若松区 40103
10 福岡県 北九州市 戸畑区 北九州市 戸畑区 40105
11 福岡県 北九州市 小倉北区 北九州市 小倉北区 40106
12 福岡県 北九州市 小倉北区 北九州市 小倉北区 40106
13 福岡県 北九州市 小倉北区 北九州市 小倉北区 40106
14 福岡県 北九州市 小倉北区 北九州市 小倉北区 40106
15 福岡県 北九州市 小倉北区 北九州市 小倉北区 40106
16 福岡県 北九州市 小倉北区 北九州市 小倉北区 40106
17 福岡県 北九州市 小倉北区 北九州市 小倉北区 40106
18 福岡県 北九州市 小倉北区 北九州市 小倉北区 40106
19 福岡県 北九州市 小倉南区 北九州市 小倉南区 40107
20 福岡県 北九州市 小倉南区 北九州市 小倉南区 40107
geometry
1 POLYGON((131.003197471 33.8...
2 POLYGON((131.008333333 33.8...
3 POLYGON((130.977422581 33.8...
4 POLYGON((130.744441712 33.8...
5 POLYGON((130.808888418 33.9...
6 POLYGON((130.759679559 33.9...
7 POLYGON((130.71166083 34.00...
8 POLYGON((130.72547917 34.01...
9 POLYGON((130.72615917 34.02...
10 POLYGON((130.851154955 33.9...
11 POLYGON((130.91506677 33.89...
12 POLYGON((130.855854449 33.9...
13 POLYGON((130.846127782 33.9...
14 POLYGON((130.852372218 33.9...
15 POLYGON((130.832898612 33.9...
16 POLYGON((130.813251115 34.0...
17 POLYGON((130.81524166 34.00...
18 POLYGON((130.809503606 34.0...
19 POLYGON((130.977270078 33.8...
20 POLYGON((131.038985383 33.8...
- select()
必要な列を選択する
fukuoka_sf %>% select(city_name) %>% head()
Simple feature collection with 6 features and 1 field geometry type: POLYGON dimension: XY bbox: xmin: 130.6732 ymin: 33.84096 xmax: 131.0239 ymax: 33.96988 epsg (SRID): 4612 proj4string: +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs city_name geometry 1 門司区 POLYGON((131.003197471 33.8... 2 門司区 POLYGON((131.008333333 33.8... 3 門司区 POLYGON((130.977422581 33.8... 4 若松区 POLYGON((130.744441712 33.8... 5 若松区 POLYGON((130.808888418 33.9... 6 若松区 POLYGON((130.759679559 33.9...
- rename()
列名を変更する
fukuoka_sf %>% rename(Prefecture = pref_name) %>% head()
Simple feature collection with 6 features and 5 fields
geometry type: POLYGON
dimension: XY
bbox: xmin: 130.6732 ymin: 33.84096 xmax: 131.0239 ymax: 33.96988
epsg (SRID): 4612
proj4string: +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
Prefecture city_name_ city_name city_name_full city_code
1 福岡県 北九州市 門司区 北九州市 門司区 40101
2 福岡県 北九州市 門司区 北九州市 門司区 40101
3 福岡県 北九州市 門司区 北九州市 門司区 40101
4 福岡県 北九州市 若松区 北九州市 若松区 40103
5 福岡県 北九州市 若松区 北九州市 若松区 40103
6 福岡県 北九州市 若松区 北九州市 若松区 40103
geometry
1 POLYGON((131.003197471 33.8...
2 POLYGON((131.008333333 33.8...
3 POLYGON((130.977422581 33.8...
4 POLYGON((130.744441712 33.8...
5 POLYGON((130.808888418 33.9...
6 POLYGON((130.759679559 33.9...
- mutate()
新規列の作成する。 新たに面積の列を追加する。
fukuoka_sf %>% mutate(area = st_area(fukuoka_sf)) %>% head()
Simple feature collection with 6 features and 6 fields
geometry type: POLYGON
dimension: XY
bbox: xmin: 130.6732 ymin: 33.84096 xmax: 131.0239 ymax: 33.96988
epsg (SRID): 4612
proj4string: +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
pref_name city_name_ city_name city_name_full city_code
1 福岡県 北九州市 門司区 北九州市 門司区 40101
2 福岡県 北九州市 門司区 北九州市 門司区 40101
3 福岡県 北九州市 門司区 北九州市 門司区 40101
4 福岡県 北九州市 若松区 北九州市 若松区 40103
5 福岡県 北九州市 若松区 北九州市 若松区 40103
6 福岡県 北九州市 若松区 北九州市 若松区 40103
area geometry
1 2700.535 m^2 POLYGON((131.003197471 33.8...
2 2229.020 m^2 POLYGON((131.008333333 33.8...
3 73715391.181 m^2 POLYGON((130.977422581 33.8...
4 62644478.607 m^2 POLYGON((130.744441712 33.8...
5 8136861.212 m^2 POLYGON((130.808888418 33.9...
6 5817.998 m^2 POLYGON((130.759679559 33.9...
- arrange()
並び変える
fukuoka_sf %>% mutate(area = st_area(fukuoka_sf)) %>% arrange(desc(area)) %>% head()
Simple feature collection with 6 features and 6 fields
geometry type: POLYGON
dimension: XY
bbox: xmin: 130.039 ymin: 33.10404 xmax: 130.9911 ymax: 33.86936
epsg (SRID): 4612
proj4string: +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
pref_name city_name_ city_name city_name_full city_code area
1 福岡県 <NA> 八女市 八女市 40210 482531309 m^2
2 福岡県 <NA> 朝倉市 朝倉市 40228 246641750 m^2
3 福岡県 <NA> 久留米市 久留米市 40203 229978807 m^2
4 福岡県 <NA> 糸島市 糸島市 40230 215059938 m^2
5 福岡県 <NA> 飯塚市 飯塚市 40205 213853733 m^2
6 福岡県 北九州市 小倉南区 北九州市 小倉南区 40107 169355048 m^2
geometry
1 POLYGON((130.724037224 33.3...
2 POLYGON((130.791228885 33.4...
3 POLYGON((130.666793891 33.3...
4 POLYGON((130.291798612 33.5...
5 POLYGON((130.743817224 33.7...
6 POLYGON((130.977422581 33.8...
- sample_n()
サンプリングする
fukuoka_sf %>% sample_n(size = 5, replace=T)
Simple feature collection with 5 features and 5 fields
geometry type: POLYGON
dimension: XY
bbox: xmin: 130.4423 ymin: 33.10404 xmax: 130.8875 ymax: 33.90448
epsg (SRID): 4612
proj4string: +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
pref_name city_name_ city_name city_name_full city_code
100 福岡県 朝倉郡 筑前町 朝倉郡 筑前町 40447
50 福岡県 <NA> 八女市 八女市 40210
65 福岡県 <NA> 宗像市 宗像市 40220
95 福岡県 遠賀郡 岡垣町 遠賀郡 岡垣町 40383
77 福岡県 <NA> 朝倉市 朝倉市 40228
geometry
100 POLYGON((130.662384721 33.5...
50 POLYGON((130.724037224 33.3...
65 POLYGON((130.443467471 33.9...
95 POLYGON((130.635642568 33.8...
77 POLYGON((130.791228885 33.4...
- group_by(), summarise()
以下の二つの図を比較することで group_by(), summarise()を利用可能性を理解することができる。
# 左の図 fukuoka_sf %>% filter(city_name_=="福岡市") %>% #group_by(city_name_) %>% #summarise() %>% ggplot() + geom_sf() # 右の図 fukuoka_sf %>% filter(city_name_=="福岡市") %>% group_by(city_name_) %>% summarise() %>% ggplot() + geom_sf()
- join関数
平成28年経済センサス‐活動調査結果 のデータを結合する
eco <- readxl::read_excel("28_kigyou_03_0.xlsx",skip=10,na = "-") eco
# A tibble: 8 x 34
地域 総数_企業数 総数_事業所数 総数_常用雇用者数 `300万円未満_企業数`
<chr> <dbl> <dbl> <dbl> <dbl>
1 福岡市 20737 49288 642559 1715
2 東 区 2728 5971 94885 198
3 博多区 5664 19689 292437 406
4 中央区 5317 12845 165969 457
5 南 区 2719 4630 38766 247
6 城南区 958 1341 9125 112
7 早良区 1900 2685 22685 155
8 西 区 1451 2127 18692 140
# ... with 29 more variables: `300万円未満_事業所数` <dbl>,
# `300万円未満_常用雇用者数` <dbl>, `300～500万円未満_企業数` <dbl>,
# `300～500万円未満_事業所数` <dbl>,
# `300～500万円未満_常用雇用者数` <dbl>,
# `500～1,000万円未満_企業数` <dbl>,
# `500～1,000万円未満_事業所数` <dbl>,
# `500～1,000万円未満_常用雇用者数` <dbl>,
# `1,000～3,000万円未満_企業数` <dbl>,
# `1,000～3,000万円未満_事業所数` <dbl>,
# `1,000～3,000万円未満_常用雇用者数` <dbl>,
# `3,000～5,000万円未満_企業数` <dbl>,
# `3,000～5,000万円未満_事業所数` <dbl>,
# `3,000～5,000万円未満_常用雇用者数` <dbl>,
# `5,000～1億円未満_企業数` <dbl>, `5,000～1億円未満_事業所数` <dbl>,
# `5,000～1億円未満_常用雇用者数` <dbl>, `1～3億円未満_企業数` <dbl>,
# `1～3億円未満_事業所数` <dbl>, `1～3億円未満_常用雇用者数` <dbl>,
# `3～10億円未満_企業数` <dbl>, `3～10億円未満_事業所数` <dbl>,
# `3～10億円未満_常用雇用者数` <dbl>, `10～50億円未満_企業数` <dbl>,
# `10～50億円未満_事業所数` <dbl>, `10～50億円未満_常用雇用者数` <dbl>,
# `50億円以上_企業数` <dbl>, `50億円以上_事業所数` <dbl>,
# `50億円以上_常用雇用者数` <dbl>
fukuoka_sf %>% left_join(eco, by =c("city_name"="地域")) %>% filter(city_name_=="福岡市") %>% ggplot() + geom_sf(aes(fill=総数_企業数)) + scale_fill_viridis(na.value="white")
ggplot2でグループごとの図を出力する
ggplotでグループごとに図を出力したい場合がある。 意外と簡単な方法でできるのでメモ。
library(ggplot2) library(dplyr)
res <- iris %>% group_by(Species) %>% do( plots= ggplot(data=.) + geom_point(aes(Sepal.Length, Petal.Length)) + ggtitle(.$Species) )
res
Species <fctr> plots <list> 1 setosa <S3: gg> 2 versicolor <S3: gg> 3 virginica <S3: gg>
res[[2]]
また株価で遊んでみた
ここで取得した株価で遊んでみた。 ooooooha.hatenablog.com
建設業関連の株価の終値の増加率の累積地を可視化してみた。 全体的に株価が上昇しているのがわかる。 塩漬けでも投資しておけば小遣い程度は稼げたんだろうな。 これをみると景気拡大期間バブル超えというのも納得できる。 実感としては増税等のせい？で減ってきている印象しかありません。
load(file = "output/建設2017-06-18.Rdata") code_list <- read.table("codelist.csv", sep=",", header = T, check.names = F) data_list %>% as_tibble() %>% arrange(code, date) %>% group_by(code) %>% mutate(prev_end_price = lag(end_price), diff_end_price = end_price / prev_end_price -1) %>% mutate(abs_rate = diff_end_price %>% abs() %>% max(na.rm=T)) %>% filter(abs_rate < 2) %>% replace_na(list(diff_end_price = 0)) %>% group_by(code) %>% mutate(cum_diff_end_price = cumsum(diff_end_price)) %>% ungroup %>% mutate(code = plyr::mapvalues(x = code, from = as.numeric(code_list$コード), to = code_list$銘柄名, warn_missing = F)) %>% ggplot() + geom_line(aes(x = date, y = cum_diff_end_price)) + geom_abline(slope = 0, linetype = 2) + facet_wrap(~ code, scale = "fixed") + theme(axis.text = element_text(size = 7), strip.text = element_text(size = 7)) + scale_x_date(labels=scales::date_format("%m"))
