【R】geodiv

1. はじめに

geodivは、ラスターやマトリックスなどの空間データの傾斜等を計算してくれるパッケージです。

2. インストール

CRANからインストールできます。

3. つかってみる

多摩地区のあるデータを国交省のページからダウンロードしてみてみます。

library(geodiv)
library(tidyverse)
library(xml2)
library(jpmesh)
library(sf)
library(raster)

tgt_file <- "FG-GML-5339-33-46-DEM5A-20161001.xml"
x <- read_xml(tgt_file)

xy_lc <-  x %>% 
  xml_find_all("/*/*/*/gml:boundedBy/gml:Envelope/gml:lowerCorner") %>% 
  xml_contents() %>% 
  as.character() %>% 
  strsplit(" ") %>% 
  unlist() %>% 
  map_dbl(as.numeric)
xy_lc

mesh <- jpmesh::coords_to_mesh(xy_lc[2], xy_lc[1])
mesh

x %>% xml_find_first("/*/*/*/gml:gridDomain/gml:Grid/gml:limits/gml:GridEnvelope")

df_dem <- x %>% 
  xml_find_first("/*/*/*/gml:rangeSet/gml:DataBlock/gml:tupleList") %>% 
  xml_contents() %>% 
  as.character() %>% 
  readr::read_delim(delim = ",", col_names = c("type", "value"))

raster_dem <- df_dem$value %>% 
  matrix(nrow = 225, ncol = 150) %>% 
  t() %>% 
  raster()
raster_dem

mesh <- 
  jpmesh::export_mesh(mesh)

bb <- mesh %>%
  st_bbox() %>%
  as.numeric()

extent(raster_dem) <- extent(bb[1], bb[3], bb[2], bb[4])

crs(raster_dem) <- sp::proj4string(as_Spatial(mesh))

plot(raster_dem)

こんな地形データです。

average surface roughness, average surface roughness, mean summit curvature を計算してみます。

sa(raster_dem)
svk(raster_dem)
ssc(raster_dem)

> sa(raster_dem)
[1] 7.006029
> svk(raster_dem)
[1] 13.28818
> ssc(raster_dem)
[1] -0.004826036

4. さいごに

地表面の情報も簡単に計算できるといいですよね。