R

西之島の標高(陰影起伏図)

raster , パッケージ

日本の平面直角座標系は x 軸が垂直方向に上の方向を正の向き

(使用するデータ)
「地理院地図」に西之島付近の噴火活動関連情報を掲載しています
公開中のデータ一覧 標高データ
撮影日 平成25年12月4日～平成27年7月28日

(準備)
データをダウンロードして解凍。作業フォルダにデータを移動。

読み込みにはdata.table::freadを使う。

平成27年7月28日のデータのextentは
-111898, -109798, 137438, 139528
（表示範囲はこれに合わせる。）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
library(raster)
#library(rasterVis)
library(data.table)
data<-c("nishinoshima_20131204_2.5m.txt","nishinoshima_20131217_2.5m.txt",
"nishinoshima_20140216_2.5m.txt","nishinoshima_20140322_2.5m.txt",
"nishinoshima_20141204_2.5m.txt","nishinoshima_20150301_2.5m.txt",
"nishinoshima_20150728_2.5m.txt","nishinoshima_20140704_1.0m_mm.txt")
for (i in 1:8){
nishi<-fread(data[i])
ymd<-paste0(substr(data[i],14,17),"年",substr(data[i],18,19),"月",substr(data[i],20,21),"日")
#欠損値 -9999 -> NA
nishi[nishi$Z==-9999, Z:=NA] 
#日本の平面直角座標系は x 軸が垂直方向に上の方向を正の向き
#x -> y ; y -> x
#data.tableの場合、with=Fが必要」
nishi<-nishi[,c(2,1,3),with=F]
#XYZをラスターに変換
nishi_rt<-rasterFromXYZ(nishi)
#平面直角座標（第14系）
#EPSG:2456
proj4string(nishi_rt) <- CRS("+init=epsg:2456")
#表示範囲を2015年07月28日にあわせる。
e <- extent(-111898, -109798,137438, 139528)
nishi_rt<-extend(nishi_rt,e)
#陰影起伏図を描く
slope <- terrain(nishi_rt,opt='slope')
aspect <- terrain(nishi_rt,opt='aspect')
hill <- hillShade(slope, aspect, 40, 270)
png(paste0(substr(data[i],1,21),".png"),width=800,height=800)
plot(hill, col=grey(0:100/100),legend=FALSE,axes=F,box=F)
plot(nishi_rt, col=terrain.colors(25, alpha=0.6),add=TRUE,axes=F,box=F)
title(paste0("西之島 ",ymd))
dev.off()
}

スライドにした。

(参考)
サムネイルもスライドするスライドショーをbxsliderで作ってみた

スライド