sf vignette 6
Rのsfパッケージの vignette6について、 これまでのvignetteで書けなかったことなどを雑多に書いてあるので、 あまりまとまりがない感じである。
EPSGとは
EPSGは空間座標参照システムとしてよく知られており、 IOGP (International Association of Oil & Gass Producers)によってメンテナンスされている。 下記の関数を利用することでEPSG一覧を見ることができる。
rgdal::make_EPSG()
sfで第二のgeometry columnを扱うには
sfオブジェクトは一つ以上のgeometry list-columnを持つことができるが、有効となるのは一つだけであり、
st_geometry()関数によって取得できる。
printメソッドを利用することでgeometry columnが有効か確認することができる。
library(sf) demo(nc, ask =F, echo = F)
nc$geom2 = st_centroid(st_geometry(nc)) print(nc, n = 2)
Simple feature collection with 100 features and 14 fields
Attribute-geometry relationship: 0 constant, 8 aggregate, 6 identity, 1 NA's
Active geometry column: geom
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
epsg (SRID): 4267
proj4string: +proj=longlat +datum=NAD27 +no_defs
First 2 features:
AREA PERIMETER CNTY_ CNTY_ID NAME FIPS FIPSNO CRESS_ID BIR74
1 0.114 1.442 1825 1825 Ashe 37009 37009 5 1091
2 0.061 1.231 1827 1827 Alleghany 37005 37005 3 487
SID74 NWBIR74 BIR79 SID79 NWBIR79 geom
1 1 10 1364 0 19 MULTIPOLYGON (((-81.47276 3...
2 0 10 542 3 12 MULTIPOLYGON (((-81.23989 3...
geom2
1 POINT (-81.49826 36.4314)
2 POINT (-81.12515 36.49101)
geometry columnを変更するにはst_geometry<-またはst_set_geometry関数を利用する。
plot(st_geometry(nc))
st_geometry(nc) <- "geom2" plot(st_geometry(nc))
st_simplifyでのトポロジー関係の維持について
st_simpligyはトポロジー保持する関数であるが、これは個々のfeature geometrieで行われる。
つまり、ポリゴンに対して実行したものはポリゴンとなる。
しかし、2つのfeaturesが境界を共有しているときにst_simplifyを実行した場合には2つのポリゴンが同じ境界を共有しているとは限らない。これは個々のポリゴンごとに単純化が行われるためである。
sfオブジェクトに対してdplyr verbsが使えない
多くの開発者はパッケージをロードせずに、sf::のように関数の属するパッケージを明示して下記のように記述する。
i = sf::st_intersects(sf1, sf2)
しかし、sfオブジェクトに対して、dplyr::selectは利用できません。
このケースではsfパッケージをロードする必要がある。
表示されるerror/warning/messageについて
座標系は緯度経度であるが xxxは平面を仮定しています（although coordinates are longitude/latitude, xxx assumes that they are planar）
sfで利用されるgeometryに関する計算においてはGEOSライブラリを利用している。このライブラリは2次元のフラットなユクリッド距離の空間座標を考慮する。しかし、例えば極を含む緯度経度データを扱い場合にはこれだけでは不十分です。
下記に例を示すが、これは間違った結果を返すこととなる。
polygon = st_sfc( st_polygon( list(rbind(c(0, 80), c(120, 80), c(240, 80), c(0, 80)))), crs =4326) pole = st_sfc(st_point(c(0, 90)), crs = 4326) st_intersects(polygon, pole)
Sparse geometry binary predicate list of length 1, where the predicate was `intersects' 1: (empty)
st_centroidは緯度経度データに対しては正しい値を返しません(
st_centroid does not give correct centroids for longitude/latitude data)
先ほどの例と同じですが、中心は2次元空間を前提としています。
st_centroid(polygon)[[1]]
この場合中心は極である必要があります。
距離は10進数を前提としています(
dist is assumed to be in decimal degrees (arc_degrees).)
このメッセージはsfでは距離の値が度で与えられていることを示しています。このメッセージを避けるには正しいユニット（単位）を与えてやる必要があります。
pt = st_sfc(st_point(c(0,0)), crs = 4326) buf = st_buffer(polygon, 1) buf = st_buffer(polygon, units::set_units(1, degree))
両方の式のオペランドがunitオブジェクトもしくはアトミックベクトルである必要があります(
both operands of the expression should be "units" objects or _Error: $ operator is invalid for atomic vectors")
このメッセージはsfパッケージをロードした後にspatstatパッケージをロードした時に現れる。パッケージ単位で作成された単位オブジェクトの場合、unitクラスも実装するspatstatのメソッドが選択される。この問題はspatstatの後にunitパッケージをロードすることで解決できる。なお、バージョン1.53-2.012のspatstatパッケージではこれが解決されています。
メルカリを利用してみて
いらないものがあったので、この機会にメルカリを利用してみた。
販売価格：8000円
これで売れたら、臨時収入ウキウキと思っていたが、 いろいろ手数料がとられて手間ばかりの割に。。。という状態に。
| メルカリに販売手数料 | 10％（800円） |
| 送料 | 900円 |
| 売上金の振込手数料 | 210円 |
| 最終振込額 | 6090円 |
8000円から諸々で1910円も取られてる。
ちなみにクロネコに持ち込んだ時に一緒にいた人もメルカリ関係での発送だったぽい。
今回は比較て販売額が大きかったからまだよかったと思うが、販売金額が小さいものだとこのもろもろひかれる額が大きすぎるような気がする。 近所の中古屋に売るより高いだろうけど、 なんでそんなにはやるのかいまいちわからなかった。
sfクラスのデータに関する隣接関係の計算に関するメモ
rにおける空間情報を扱うにあたり、最近はsfパッケージが速くて便利である。
さて、sfクラスののデータに対して隣接関係を用いて分析等を行いたいときはどうすればよいのだろうか。 旧来のspパッケージでいうとspdepパッケージがそれらに関する関数が充実していた。
sfクラスのデータに対してはどのように行うかについてspdepのvignetteに書かれている。 github.com
結論から言うと、spdepパッケージは効率よく計算できるように設計してあり、 現状ではsf --> sp クラスに変換して行うことが処理速度上好ましいようである。
prophetでガソリン価格の予測
prophetパッケージを利用してみようということで、 とりあえずvignetteにあるものをデータを変えて実行してみた。
ガソリン価格の予測を行う
資源エネルギー庁のHPから給油所小売価格調査（ガソリン、軽油、灯油）の週次データを取得し、 レギュラーガソリン価格の推移を利用して今後のガソリン価格の予測を行う。 当たり前だが、データ依存による影響がみられており、適切なデータ選定が必要なことがわかる。
library(tidyverse) library(prophet) library(magrittr) file_url <- "https://github.com/ogwf/gasoline_value/blob/master/gasoline_prices.xls?raw=true" tmp <- paste0(getwd(), "/gasoline_prices.xls") download.file(url = file_url, destfile = tmp, mode = "wb") gasoline <- readxl::read_xls(path = tmp, sheet = "レギュラー", col_types = "numeric") names(gasoline) %<>% stringr::str_replace_all(pattern = "[:space:]", "") gasoline <- gasoline %>% filter(!is.na(調査日)) %>% mutate(調査日 = 調査日 %>% as.Date(origin = "1900-01-01"))
これまでのガソリン価格の推移はこのような状況である。
ggplot(gasoline, aes(調査日, 全国)) + geom_path()
これをprophetを利用して予測してみると。。。
res <- gasoline %>% transmute(ds = 調査日, y = 全国) %>% prophet()
Initial log joint probability = -9.47537 Optimization terminated normally: Convergence detected: relative gradient magnitude is below tolerance
future <- make_future_dataframe(res, periods = 365) plot(res, predict(res, future))
prophet_plot_components(res, predict(res))
2010年以降は、価格の変動が激しいこともありうまくトレースできておらず、近年は下落傾向となっている。
直近一年間のガソリン価格を使って予測してみる。
近年のガソリン価格の変動が激しくうまく予測できていないようなので、 期間を直近一年間に絞り同様に予測してみる。
直近のガソリン価格の推移は下記の通り。
gasoline %>% arrange(desc(調査日)) %>% slice(1:52) %>% ggplot(aes(調査日, 全国)) + geom_path()
予測結果は
res <- gasoline %>% arrange(desc(調査日)) %>% slice(1:52) %>% transmute(ds = 調査日, y = 全国) %>% prophet()
Initial log joint probability = -2.03264 Optimization terminated normally: Convergence detected: relative gradient magnitude is below tolerance
future <- make_future_dataframe(res, periods = 60) plot(res, predict(res, future))
prophet_plot_components(res, predict(res))
tmapパッケージはしばらく更新がないらしい
先日tmap(version 1.11)の更新通知が来ていた。 たまたまNewsを確認したところ下記のような記載があり、 sfパッケージがspに置き換わるまでしばらく更新はなさそうである。
NOTE: this will be the last version before the major update (in which sf fully replaces sp)
issuesには現時点で25項目上がっているが、これらは今後解決していくのだろうか。
