Rでおっぱい山を分析する

これはFOSS4G Advent Calendar 2017 24日目の記事です。

「おっぱい山」とは、おっぱいの形をした山（英：breast-shaped hill）のことです。（出典： Wikipedia）
今回は統計ソフト「R」を使って各地のおっぱい山を分析してみたいと思います。

なぜ、そんなことするかだって？
そこにおっぱい山があるからさ
（ジョージ・マロリー）

Rと地理院タイルで作成した北アルプスCS立体図

まずは、おっぱい山を探す所からです。地理院地図の検索機能を利用してみましょう。

検索窓に"おっぱい山"と入力してみます。
しかしながら、0件、見つからず。東大の協力も力及ばずと言った結果となりました...

仕方ないのでググってデータを集めました。潜在的なおっぱい山はもっとあるはずですが、なかなか表に出てこないようです。
データ：https://gist.github.com/tmizu23/29c24df1662bbf976e2407f843f43c15

area	name	lat	lon
北海道上士幌町	西クマネシリダケ、ピリベツ岳	43.524422	143.229961
岩手県八幡平市	七時雨山	40.070009	141.107197
福島県郡山市	安達太良山	37.621130	140.287864
福島県天栄村	二岐山	37.247941	139.965477
東京都小笠原村母島	乳房山	26.659560	142.161412
東京都小笠原村父島	乳頭山	27.096348	142.210046
埼玉県小川町	笠山	36.018509	139.189873
山梨県甲府市	要害山	35.703142	138.599333
広島県安芸太田町	574ピーク	34.583811	132.211576
沖縄県宮古島市	69.1ピーク	24.750397	125.405974

おっぱい山の可視化

データ分析の基本は可視化です。
これらのおっぱい山をRと地理院タイルを利用して、可視化してみましょう。

地理院タイルをRで読み込めるようにベクトルタイル Advent Calendar 2017 16日目の記事で作成したコードを改造しました。指定した範囲の画像タイル、PNG標高タイルをダウンロードし、Rasterクラスに格納するので、表示・分析に利用できます。

コード：https://github.com/tmizu23/gsitiles.R

使い方

最初に、以下のRライブラリをインストールしてください。

install.packages(c("raster", "curl", "png", "jpeg", "spatialEco"),dependencies=T)

関数をsource("gsitiles.R")で読み込んでから、こんな感じで使用します。

#使用例
source("gsitiles.R")

#標高データを取得・表示・投影変換・書き出し
dem<-getTile(137.665215,36.294582,137.7,36.314582,15,type="dem5a")
plot(dem)
demUTM53<-projectRaster(dem, crs=CRS("+init=epsg:3099"))
writeRaster(demUTM53,"demUTM53.tif")

#標準地図とCS立体図を乗算合成・表示・geoTIFF書き出し
cs<-makeCS(dem,shaded=TRUE)
stdmap<-getTile(137.665215,36.294582,137.7,36.314582,15,type="std")
stdmap_cs<-overlayColor(stdmap,cs,type="multiply")
plotRGB(stdmap_cs)
writeRaster(stdmap_cs,"stdmap_cs.tif",datatype="INT1U") #datatypeを指定

#陰影起伏図の作成・表示
slope<-terrain(dem,"slope")
aspect<-terrain(dem,"aspect")
shade<-hillShade(slope,aspect)
plot(shade,col=grey(0:100/100),legend=F)

#航空写真と陰影起伏を透過表示
photo<-getTile(137.665215,36.294582,137.7,36.314582,15,type="photo")
plot(shade,col=grey(0:100/100),legend=F)
plotRGB(photo,alpha=180,add=T)

• getTile関数で取得したい範囲の左下と右上の座標、ズームレベル、タイルの種類を指定します。（左下と右上の座標が同じ場合は、その座標を含む1タイルを取得します。）
• タイルの種類は、標準地図："std"、シームレスフォト："photo"、標高PNG："dem10b","dem5a","dem5b"を選択できます。提供されているズームレベルと範囲は、地理院タイルのページを参照ください。
• crop=TRUEを指定すると、指定範囲で切り抜きます。デフォルト(FALSE)では、指定範囲を含むタイルの範囲になります。

この他に以下の関数を用意しました。

• makeCS(dem,shaded=TRUE)：

　　　標高データからCS立体図を作成。shaded=TRUEで陰影起伏図とオーバレイ合成したCS立体図を作成

• overlayColor(stackA,stackB,type="overlay | multiply")：

　　　RGB値の入ったstackクラス同士を画像合成。種類は、オーバレイ合成："overlay"、乗算合成："multiply"。（stackAとstackBは同じ範囲、ズームレベルの必要があります）

CS立体図については、こちらの記事を参考にさせて頂きました。

http://waigani.hatenablog.jp/entry/2017/12/03/225834
https://qiita.com/frogcat/items/7e91d3070a7a8d3e2c94

おっぱい山の可視化

前置きが長くなりましたが、可視化の結果はこちらです。

日本を代表するおっぱい山です。両ピークに名前が付いているのが特徴です。

なかなかの良型です。おっぱい山として文句はありません。

日本百名山からのエントリーです。片乳ながら、乳頭のリアリティを追求しています。

デコルテが激しく何かを主張しています。もっとリスペクトされても良いかもしれません。

母島のおっぱい山です。あるべくしてある。必然の結果です。

父島のおっぱい山です。父にも乳はあるんです。

地元では名が知られているが垢抜けない。伸びしろに期待です。

史跡と一体となった妖麗なおっぱい山です。誘われたら断れません。

一見、変哲もなくつまらない。それでいいんです。

言わんとすることは分かるんです。でも危険と裏腹な恐ろしさがあります。

まとめ

地理院タイルとRで、おっぱい山の可視化が簡単にできるようになりました。
しかしながら、我々が為すべきことは山積みです。

どのようにして、おっぱい山ができたのか？
なぜ、おっぱい山はそこにあるのか？
これから、おっぱい山はどうなってゆくのか？

おっぱい山の頂き目指し、これからも一歩一歩登っていきたいと思います。
それでは、良いクリスマスを！