メモ：sfのオブジェクトをtidyr::gather()でtidyにするとサイズはでかくなる

R sf kokudosuuchi 国土数値情報 GIS

キャプテンアメリカがこんな記事を上げていたので、ちょっと気になってたことをメモ。

対象のデータ

先日、kokudosuuchiのアップデートを告知しましたが、

国土数値情報のデータの中には、

各年度別最深積雪（A22_01XXXX）XXXX：西暦
(http://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-A22-m.html)

みたいなフィールドがあります。これは、例えば、列名がA22_012000なら「2000年度の最深積雪」という意味になります。ただし、A22_010001は西暦1年のデータかといえばそんなことはなく全データの平均値、みたいな例外もあるので注意が必要です…

実際のデータはこんな感じです。

library(kokudosuuchi)
library(dplyr)

# 67MBあるので注意
d <- getKSJData("http://nlftp.mlit.go.jp/ksj/gml/data/A22-m/A22-m-14/A22-m-14_01_GML.zip")

names(d)
#> [1] "A22-m-14_MaxSnowDepth_01"      "A22-m-14_TotalSnowfall_01"    
#> [3] "A22-m-14_MinTemperature_01"    "A22-m-14_MeanWindSpeed_01"    
#> [5] "A22-m-14_SnowDisaster_01"      "A22-m-14_SnowCoolingSystem_01"

今回はA22-m-14_MaxSnowDepth_01を見てみましょう。上に挙げたA22_01XXXXの列がたくさんありますね。

names(d$`A22-m-14_MaxSnowDepth_01`)
#>   [1] "A22_000001" "A22_000002" "A22_000003" "A22_010001" "A22_010002"
#>   [6] "A22_012013" "A22_012012" "A22_012011" "A22_012010" "A22_012009"
#>  [11] "A22_012008" "A22_012007" "A22_012006" "A22_012005" "A22_012004"
#>  [16] "A22_012003" "A22_012002" "A22_012001" "A22_012000" "A22_011999"
#>  [21] "A22_011998" "A22_011997" "A22_011996" "A22_011995" "A22_011994"
#>  [26] "A22_011993" "A22_011992" "A22_011991" "A22_011990" "A22_011989"
#>  [31] "A22_011988" "A22_011987" "A22_011986" "A22_011985" "A22_011984"
#>  [36] "A22_011983" "A22_011982" "A22_011981" "A22_011980" "A22_011979"
#> ...

translateKSJData()はこういう年度もいい感じに変換してくれます。

max_snow <- translateKSJData(d$`A22-m-14_MaxSnowDepth_01`)

glimpse(max_snow)
#> Observations: 505
#> Variables: 149
#> $ 観測点名                     <chr> "東中山", "当別", "無意根", "月形除雪ST", "厚田", "八幡", ...
#> $ 観測点の所在地               <chr> "虻田郡喜茂別町字川上中山峠", "石狩郡当別町対雁43の１", "札幌市南区定山渓国有林...
#> $ 観測点の管理者               <chr> "北海道開発局", "北海道開発局", "北海道開発局", "北海道開発局", "北海道開...
#> $ `各年度別最深積雪（平均値）` <dbl> 289.50000, 109.28571, 268.92857, 181.00000, 164.428...
#> $ `各年度別最深積雪（最大値）` <dbl> 387, 197, 344, 255, 284, 172, 136, 148, 275, 109, 2...
#> $ 各年度別最深積雪2013年度     <dbl> 286, 166, 265, 177, 164, 159, 136, 137, 219, 65,...
#> $ 各年度別最深積雪2012年度     <dbl> 288, 187, 294, 202, 192, 172, 121, 141, 264, 94,...
#> $ 各年度別最深積雪2011年度     <dbl> 293, 197, 238, 255, 204, 165, 119, 119, 275, 91,...
#> $ 各年度別最深積雪2010年度     <dbl> 247, 93, 193, 120, 116, 133, 118, 140, 144, 99, ...
#> $ 各年度別最深積雪2009年度     <dbl> 228, 47, 215, 175, 226, 120, 118, 107, 182, 56, ...
#> $ 各年度別最深積雪2008年度     <dbl> 233, 32, 266, 87, 88, 104, 89, 80, 149, 38, 136,...
#> $ 各年度別最深積雪2007年度     <dbl> 235, 43, 230, 143, 169, 123, 129, 95, 195, 90, 1...
#> $ 各年度別最深積雪2006年度     <dbl> 382, 99999999, 338, 203, 115, 83, 93, 83, 198, 4...
#> $ 各年度別最深積雪2005年度     <dbl> 387, 99999999, 344, 254, 284, 109, 104, 148, 236...
#> $ 各年度別最深積雪2004年度     <dbl> 273, 99999999, 280, 193, 249, 142, 132, 146, 109...
#> ...

さて、ではこれをtidyなデータ（整然データ）にするために、tidy::gather()で「各年度別最深積雪XXXX年度」という列を集めます。

max_snow_tidy <- tidyr::gather(max_snow,
                               key = "fiscal_year",
                               value = "max_depth",
                               dplyr::matches("各年度別最深積雪\\d+年度")) %>%
  mutate(fiscal_year = readr::parse_number(fiscal_year))

できたデータはこんな感じ。縦長になっていますね。

glimpse(max_snow_tidy)
#> Observations: 71,710
#> Variables: 9
#> $ 観測点名                     <chr> "東中山", "当別", "無意根", "月形除雪ST", "厚田", "八幡", ...
#> $ 観測点の所在地               <chr> "虻田郡喜茂別町字川上中山峠", "石狩郡当別町対雁43の１", "札幌市南区定山渓国有林...
#> $ 観測点の管理者               <chr> "北海道開発局", "北海道開発局", "北海道開発局", "北海道開発局", "北海道開...
#> $ `各年度別最深積雪（平均値）` <dbl> 289.50000, 109.28571, 268.92857, 181.00000, 164.428...
#> $ `各年度別最深積雪（最大値）` <dbl> 387, 197, 344, 255, 284, 172, 136, 148, 275, 109, 2...
#> $ A22_050001                   <int> 1, 1, 1, 1, 3, 3, 1, 1, 3, 1, 3, 1, 3,...
#> $ geometry                     <S3: sfc_POINT> POINT (141.096111111111 42.....
#> $ fiscal_year                  <dbl> 2013, 2013, 2013, 2013, 2013, 2013, 20...
#> $ max_depth                    <dbl> 286, 166, 265, 177, 164, 159, 136, 137...

実際にデータを可視化したり分析したりするには、もうちょっとデータの前処理が必要そうです（欠損値が99999999で表されていたりとか）が、今回はそれがメインではないのでここまでにしておきます。可視化するには以下のような感じにすればいいはず。

ggplot(max_snow_tidy) +
  geom_sf(aes(colour = max_depth)) +
  facet_wrap(~fiscal_year)

で、何が気になっているかというと…オブジェクトのサイズです。オブジェクトのサイズは、pryr::object_size()で調べられます。

pryr::object_size(max_snow_tidy)
#> 794 kB

pryr::object_size(max_snow)
#> 5.04 MB

明らかに膨れ上がっています。これは、縦長になった分だけgeometryが複製されているからです。これはまだ点データなのでよかったですが、ポリゴンデータとかだともっと大きくなるでしょう。

tidyなのはいいことですが、こういう時空間データのよい持ち方がいまだによくわかりません。たぶん↓が解決してくれるのかな…？