LoginSignup
11
14

More than 3 years have passed since last update.

国土地理院基盤地図情報をRを使ってGISデータに変換する

Last updated at Posted at 2021-04-11

04.png

国土地理院の基盤地図情報は独自のxmlファイルで記述されています。
基本項目は国土地理院から「基盤地図情報ビューア」というコンバーターが提供されており、数値標高モデルはエコリスさんから「基盤地図情報標高DEMデータ変換ツール」というものが提供されています。普通はこの2種のコンバーターを使用して、基盤知事情報のデータをQGISなどで使用できる一般的なGISデータ(tiffやshp形式)に変換します。

しかし、これらのコンバーターはいずれもWindows専用で、MacOSやLinuxでは使えません(Wineを使用する方法もありますがうまくいかないことも多いです)。そこで、Rのfgdrパッケージを使用して基盤地図情報のデータの変換を行います。

パッケージ読み込み

最低限必要なパッケージは次のとおりです。

# パッケージ読み込み
library(tidyverse)
library(sf)
library(raster)
library(fgdr)

ダウンロードしたデータの処理

ダウンロードした基盤地図情報のzipファイルを解凍するとさらにzipファイルが現れます。
01.png

一つ一つ解凍するのは大変なのでunzipコマンドで一括解凍します。unzipコマンドのワイルドカードはちょっと変わっていて、クォーテーションマークで囲みます。

unzip '*zip'

解凍が終わるとxmlファイルが現れます。
02.png

数値標高モデルの変換

まずは、標高データをGeoTiffに変換します。
すでに国土地理院の数値標高モデルデータをラスタとしてRで扱うで紹介されている手法です。
あらかじめデータを格納するmapディレクトリを作成し、解凍したxmlファイルの入ったディレクトリをその直下に置いています。

#1 ダウンロードして解答したxmlファイルのリストを作成します。
xml <- 
  list.files(path="map/PackDLMap_DEM/", 
             pattern=".*\\.xml", 
             full.names=T) 
#2 最初のラスタを読み込む
dem <- read_fgd_dem(xml[1], resolution = 10 , return_class = "raster")  
#3 for文でxmlファイルをラスタオブジェクトに変換して結合
for(i in 2:length(xml) ) {
  dem1 <- read_fgd_dem(xml[i], 
                       resolution = 10 ,
                       return_class = "raster" )
  dem <- merge(dem , dem1,tolerance = 0.5)
}
#4  海域 = -9999に0を代入
dem[dem<0] <- 0
# 地図を表示
plot(dem)
#5 座標系を変換
dem_12k <- projectRaster(dem , crs = CRS("+init=epsg:2454") )
#6 tiffに出力
writeRaster(dem_12k, "map/FranoDEM_JGD2000_12k.tif")

03.png

基本項目

続いて基本項目(道路や河川)の変換です。こちらもベースは「Rで国土地理院 基盤地図情報データを扱う」を参考にしています。

こちらはpurrr::map関数によるループ処理を利用しています。

# 水域
list.files(path = "map/PackDLMap_kihon/",  
           pattern = "WA",
           full.names = T) %>%  # 基盤地図情報のxmlファイルをリスト化
  map(read_fgd) %>%  # xmlファイルにread_fgdを適用
  bind_rows() %>%  # xmlファイルを結合
  st_transform(2454) %>%  # JGD2000平面直角座標系12系に変換
  st_write(driver = "GPKG" ,
           append = FALSE , "map/WA_JGD2000_12k.gpkg")  # ファイルを書き出し

# 河川
list.files(path = "map/PackDLMap_kihon/",  
           pattern = "WL",
           full.names = T) %>%  # 基盤地図情報のxmlファイルをリスト化
  map(read_fgd) %>%  # xmlファイルにread_fgdを適用
  bind_rows() %>%  # xmlファイルを結合
  st_transform(2454) %>%  # JGD2000平面直角座標系12系に変換
  st_write(driver = "GPKG" ,
           append = FALSE , "map/WL_JGD2000_12k.gpkg")  # ファイルを書き出し


# 等高線
list.files(path = "map/PackDLMap_kihon/",  
           pattern = "Contr",
           full.names = T) %>%  # 基盤地図情報のxmlファイルをリスト化
  map(read_fgd) %>%  # xmlファイルにread_fgdを適用
  bind_rows() %>%  # xmlファイルを結合
  st_transform(2454) %>%  # JGD2000平面直角座標系12系に変換
  st_write(driver = "GPKG" ,
           append = FALSE , "map/Contr_JGD2000_12k.gpkg") 

#建物
list.files(path = "map/PackDLMap_kihon/",  
                  pattern = "BldA",
                  full.names = T) %>%  # 基盤地図情報のxmlファイルをリスト化
  map(read_fgd) %>%  # xmlファイルにread_fgdを適用
  bind_rows() %>%  # xmlファイルを結合
  st_transform(2454) %>%  # JGD2000平面直角座標系12系に変換
  st_write(driver = "GPKG" ,
           append = FALSE , "map/BldA_JGD2000_12k.gpkg") 

# 鉄道
list.files(path = "map/PackDLMap_kihon/",  
           pattern = "RailCL",
           full.names = T) %>%  # 基盤地図情報のxmlファイルをリスト化
  map(read_fgd) %>%  # xmlファイルにread_fgdを適用
  bind_rows() %>%  # xmlファイルを結合
  st_transform(2454) %>%  # JGD2000平面直角座標系12系に変換
  st_write(driver = "GPKG" ,
           append = FALSE , "map/RailCL_JGD2000_12k.gpkg") 

# 市町村域_polygon
list.files(path = "map/PackDLMap_kihon/",  
                  pattern = "AdmArea",
                  full.names = T) %>%  # 基盤地図情報のxmlファイルをリスト化
  map(read_fgd) %>%  # xmlファイルにread_fgdを適用
  bind_rows() %>%  # xmlファイルを結合
  st_transform(2454) %>%  # JGD2000平面直角座標系12系に変換
  group_by(name) %>%  # 市町村名でグループ化
  summarize() %>%  # 市町村名で結合
  st_write(driver = "GPKG" ,
           append = FALSE , "map/AdmArea_JGD2000_12k.gpkg") 

# 市町村域_Line
list.files(path = "map/PackDLMap_kihon/",  
           pattern = "AdmBdry",
           full.names = T) %>%  # 基盤地図情報のxmlファイルをリスト化
  map(read_fgd) %>%  # xmlファイルにread_fgdを適用
  bind_rows() %>%  # xmlファイルを結合
  st_transform(2454) %>%  # JGD2000平面直角座標系12系に変換
  st_write(driver = "GPKG" ,
           append = FALSE , "map/AdmBdry_JGD2000_12k.gpkg") 

ポリゴンを属性で結合する方法でつまずいたのですが、なんのことはなく、summarize()で良いということがわかりました。こちらを参考にしました。(Rでポリゴンを融合する(実践編)

以上の作業で出力したtiffやgpkgファイルを読み込んで、QGISで表示できました。
04.png

11
14
3

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
11
14