R

西之島の標高(平成27年7月28日) (変更あり)

rgl , raster , rasterVis パッケージ

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

(使用するデータ)
「地理院地図」に西之島付近の噴火活動関連情報を掲載しています　－７月２８日に撮影した西之島－
標高データ 撮影日:平成27年7月28日

(準備)
データをダウンロードして解凍。作業フォルダにデータを移動。
作業フォルダは /home/user/R/work とします。

今回、読み込みにはdata.table::freadを使ってみます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
library(rgl)
library(data.table)
library(raster)
library(rasterVis)
#今回、読み込みにはdata.table::freadを使ってみます。
nishi <- fread("/home/user/R/work/nishinoshima_20150728_2.5m.txt")
head(nishi)
# 特定の列の要素を変更
#-9999をNAに変更
#data.tableでは
nishi[nishi$Z==-9999, Z:=NA] 
#data.frameでは
#nishi[,3][x[,3]==-9999] <- NA
#欠損値のあるデータを削除
#nishi<-na.omit(nishi)
#日本の平面直角座標系は x 軸が垂直方向に上の方向を正の向き
#x -> y ; y -> x
#nishi<-nishi[,c(2,1,3)]
#data.tableの場合、with=Fが必要」
nishi<-nishi[,c(2,1,3),with=F]
#標高データの仕様 西之島の平面直角座標（第14系） EPSG:2456
nishi_rt<-rasterFromXYZ(nishi)
proj4string(nishi_rt) <- CRS("+init=epsg:2456")
nishi_rt

class : RasterLayer
dimensions : 836, 840, 702240 (nrow, ncol, ncell)
resolution : 2.5, 2.5 (x, y)
extent : -111898, -109798, 137438, 139528 (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:2456 +proj=tmerc +lat_0=26 +lon_0=142 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
data source : in memory
names : Z
values : -2.8, 149.7 (min, max)

1
2
#rasterVis::levelplot
levelplot(nishi_rt, col.regions=terrain.colors)

1
2
3
4
5
6
7
8
9
10
11
#解像度を落とす(fact = 5)
nishi_rt5 <- aggregate(nishi_rt, fact = 5)
#サイズを指定
r3dDefaults$windowRect <- c(50,50,500,500)
#3D表示(R上では動きます。)
plot3D(nishi_rt5, adjust=F, col=terrain.colors)
rgl.snapshot( "plot3dimage.png", fmt="png")
#
#decorate3d()
#postscriptで保存する場合。
#rgl.postscript("plot3dimage.ps")

スクリーンショット

以下、コードのみ 

1
2
3
4
5
6
7
8
9
10
11
#解像度を落とす(fact = 8)
nishi_rt8 <- aggregate(nishi_rt, fact = 8)
plot3D(nishi_rt8, adjust=F)
decorate3d( xlab = "x", ylab = "y", zlab = "",  box = F, axes = T)
#動く HTML ファイル
#もっと解像度を落とす(fact = 20)
nishi_rt20<- aggregate(nishi_rt, fact =20)
plot3D(nishi_rt20, adjust=F)
decorate3d( xlab = "x", ylab = "y", zlab = "",  box = F, axes = F)
#動く HTML ファイルとして保存
writeWebGL(width=356, height=356)