raster パッケージ
国土地理院標高データのprojectionをランドサットのprojectionに変更する。
(使用するデータ)
「地理院地図」に西之島付近の噴火活動関連情報を掲載しています
標高データ 撮影日 平成27年7月28日
(準備)
データをダウンロードして解凍。作業フォルダにデータを移動。
読み込みにはdata.table::freadを使う。
1 2 3 4 5 6 7 8 9 10 11 12 13
| library(raster) #library(rasterVis) library(data.table) nishi<-fread("nishinoshima_20150728_2.5m.txt") #欠損値 -9999 -> NA nishi[nishi$Z==-9999, Z:=NA] #日本の平面直角座標系は x 軸が垂直方向に上の方向を正の向き coordinates(nishi) <- c( "Y","X") proj4string(nishi) <- CRS("+init=epsg:2456") #projectionの変更 #ランドサット：西之島 測地系(datum):WGS84, 投影法:UTM、zone:54 res <- spTransform(nishi, CRS("+proj=utm +zone=54 +ellps=WGS84")) res
|
class : SpatialPointsDataFrame
features : 702240
extent : 487131.6, 489244.8, 3012584, 3014687 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=54 +ellps=WGS84
variables : 1
names : Z
min values : -2.8
max values : 149.7
1 2 3 4 5 6 7 8 9 10
| #SpatialPointsDataFrameをrasterへ x <- data.frame(res@coords[,1],res@coords[,2],res@data) colnames(x) <- c('x', 'y', 'z') # set up an 'empty' raster, here via an extent object derived from your data e <- extent(x[,1:2]) #解像度を落とす。元はdimensions : 836, 840, 702240 (nrow, ncol, ncell) r <- raster(e, ncol=400, nrow=400) nishi_rt <- rasterize(x[,1:2], r, x[,3], fun=mean) proj4string(nishi_rt) <- CRS("+proj=utm +zone=54 +ellps=WGS84") nishi_rt
|
class : RasterLayer
dimensions : 400, 400, 160000 (nrow, ncol, ncell)
resolution : 5.282908, 5.25812 (x, y)
extent : 487131.6, 489244.8, 3012584, 3014687 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=54 +ellps=WGS84
data source : in memory
names : layer
values : -2.8, 149.6 (min, max)
nishi_rtを保存します。(nishi20150728と名前を変える)
1 2
| nishi20150728<-nishi_rt save("nishi20150728",file="nishi20150728.Rdata")
|