1 目的

  • Hillshade effects(アクセス日:2024-04-26)を参考に光源の位置と方向を活用して,地形の陰影を生成し,山や谷の輪郭を浮かび上がらせる(陰影起伏図).
  • 標高図と陰影起伏図を組み合わせ,色付き陰影起伏図(陰影段彩図)を作成する.
    • ラスターデータ:四角いドット(ピクセル)のグリッドで表現されるデータ形式.1
library(tidyverse)
#地図(静岡県)ダウンロード
library(rnaturalearth) 
#標高データダウンロード
library(elevatr) 
#標高データ操作
library(terra) 
#標高データの色付け
library(tidyterra)
#スケールの複数使用 
library(ggnewscale)
#地図(伊豆地域,横浜市)の読込
library(sf)

2 標高図の作成

2.1 標高データのダウンロード

  • 標高を示す地域として静岡県を選定.Natural Earthに収録されている静岡県ポリゴンをne_states()で読込.2
    • 静岡県は日本の都道府県のため国はJapanを指定.
    • sfクラスを指定.
    • Japan_mapと名付ける.
Japan_map <- 
  ne_states("Japan", returnclass="sf")

Shizuoka <- Japan_map %>% 
  filter(name=="Shizuoka")
  • get_elev_raster():指定した位置の標高データ(ラスターデータ)をダウンロード.
    • locations:位置情報を指定.今回は静岡県(Shizuoka
    • z:解像度(ズームレベル)の指定.数を大きくすると解像度が増すが,時間がかかるため注意
    • cliplocationsを設定した場合,指定した位置のデータを切り取る(bboxを設定した場合,境界ボックスに合わせてデータを切り取る).
elev <- 
  get_elev_raster(locations=Shizuoka, 
                  z=9, clip="locations")
  • ggplot()で可視するため,elevをデータフレーム(as.data.frame())に変更.
    • xy=TRUEを指定することで,座標情報をデータフレームの列として取り込む.
elev_df <- elev %>% 
  as.data.frame(xy=TRUE)

2.2 データ処理 

  • elev_dfの3列目が標高を示すため,列名をelvationに変更.
  • データフレームの欠損値(NA)削除(drop_na(everything())).
elev_df <- elev_df %>%
  rename("elevation"=3)

elev_df <- elev_df %>% 
  drop_na(everything())

2.3 標高図

  • geom_raster():標高(elev_dfelevation列)の可視化.詳細
    • scale_fill_hypso_c():標高の色付け.パレットはc3t1を選択.
      • 利用可能パレット
      • 標高の単位はメートル(m).
ggplot() +
  geom_raster(data=elev_df,
              aes(x, y, fill=elevation))+
  scale_fill_hypso_c(palette="c3t1")+
  labs(fill="m",
       caption="出典:Natural Earth,
       AWS Terrain Tiles,OpenTopography API")+
  coord_sf(crs="EPSG:6668")+
  theme_void() 

3 陰影起伏図の作成

3.1 陰影計算

  • terrain():標高データ(elv)から地形特性を計算する関数.
    • 傾斜(slope)をラジアン単位(radians)で計算するように指示.
    • 斜面方位(aspect)をラジアン単位(radians)で計算するように指示.
sl <- elev %>%
  terrain("slope", unit="radians")

asp <- elev %>% 
  terrain("aspect", unit="radians")
  • shade():陰影計算.
    • angle:光源高度45度に指示.
    • direction:光源方向315度に指示.

注意:このままshade()を指示すると,

inherits(slope, "SpatRaster") は TRUEではありません

と出るためrasterSpatRasterに変更.

sl <- sl %>% 
  as("SpatRaster")
asp <- asp %>% 
  as("SpatRaster")

hill_shade <- 
  shade(sl, asp, 
      angle=45, 
      direction=315,
      normalize=TRUE)
  • hill_shadeもデータフレーム(as.data.frame())に変更.
hill_shade_df <- hill_shade %>% 
  as.data.frame(xy=TRUE)

3.2 陰影起伏図

  • 陰影起伏(hill_shade_dfhillshade列)を紫色(Purples)で表現.
    • show.legend=FALSE:凡例非表示.
ggplot() +
  geom_raster(data=hill_shade_df,
              aes(x, y, fill=hillshade),
              show.legend=FALSE)+
  scale_fill_distiller(palette="Purples")+
  labs(caption="出典:Natural Earth,
       AWS Terrain Tiles,OpenTopography API")+
  coord_sf(crs="EPSG:6668")+
  theme_void()

4 陰影段彩図

  • 標高図と陰影起伏図を組み合わせ,色付き陰影起伏図(「陰影段彩図」)を作成.
  • 地図の順番.
    • 最下層:陰影起伏図(hill_shade_df).
    • new_scale_fill():重ねる標高地図の色付けを活かすため.
    • 最上層:標高図(elev_df).alphaを利用して透過表示.
ggplot() +
  geom_raster(data=hill_shade_df,
              aes(x, y, fill=hillshade),
              show.legend=FALSE)+
  scale_fill_distiller(palette="Purples")+
  new_scale_fill() +
  geom_raster(data=elev_df,
              aes(x, y, fill=elevation),
              alpha=0.7) +
  scale_fill_hypso_c(palette="c3t1")+
  labs(fill="m",
       caption="出典:Natural Earth,
       AWS Terrain Tiles,OpenTopography API") +
  coord_sf(crs="EPSG:6668") +
  theme_void() 

5 応用(等高線の追加) 

  • 同じ標高を線で表すことで等高線を作成.

5.1 伊豆地域

上の静岡県との違い

  1. Natural Earthから伊豆半島のポリゴンを正確に切り抜くのが難しいこと.
  1. 伊豆地域のポリゴンをIzuと名付ける.
  • z:解像度を8に指定して,標高データ(ラスターデータ)をダウンロード.
  • 後は原則,静岡県のケースと同じであるためRコードは省略.
Izu_elev <- 
  get_elev_raster(locations=Izu,
                  z=8, clip="locations")

5.2 伊豆地域の陰影段彩図

  1. 陰影起伏図(データフレーム)をhillshade_Izu_dfと名付ける.静岡県同様に標高図(Izu_elev_dfと名付けている)と組み合わせ陰影段彩図を作成しようとしたところ,geom_raster()ではなくgeom_tile()を勧められたため,変更.
  • geom_contour():等高線の作成.2変数(x軸とy軸)に対するz軸の値を表現.ここで等高線の高さ.内側がより高い値.
    • 線の色(color)と太さ(linewidth)の指定.
ggplot() +
  geom_tile(data=hillshade_Izu_df,
              aes(x, y, fill=hillshade),
              show.legend=FALSE)+
  scale_fill_distiller(palette="Purples")+
  new_scale_fill() +
  geom_tile(data=Izu_elev_df,
              aes(x, y, fill=elevation),
              alpha=0.7)+
  scale_fill_hypso_c(palette="c3t1")+
  geom_contour(data=Izu_elev_df, 
               aes(x, y, z=elevation), 
               color="gray70", linewidth=0.2)+
  labs(fill="m",  
       caption="出典:国土数値情報行政区域データ,
       AWS Terrain Tiles,OpenTopography API") +
  coord_sf(crs="EPSG:6668") +
  theme_void() 

6 横浜市

Kanagawa_map<-
  read_sf("N03-20240101_14_GML/N03-20240101_14.shp")

Kanagawa_map %>% 
  filter(N03_004=="横浜市") ->
  Yokohama
  • 標高データのダウンロード.解像度を10に指定.
  • 後は原則,静岡県のケースと同じであるためRコードは省略.
Yokohama_elev <- 
  get_elev_raster(locations=Yokohama, 
                  z=10, clip="locations")

6.1 横浜市の陰影段彩図

7 参考ウェブサイト

RによるGISデータの可視化

アイコン

  • Font Awesome 6.5.1

  1. 標高・傾斜度メッシュデータの可視化では国土数値情報標高・傾斜度メッシュデータを用いて標高を可視化.このメッシュデータは基盤地図情報数値標高モデルのラスターデータをシェープファイル(ベクターデータ:ポイント,ライン,ポリゴンで表現されるデータ形式)に変換されている.↩︎

  2. Natural Earthの利用方法の詳細はNatural Earthの地理空間データの可視化参照.↩︎