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

![04.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/33e73886-1bb7-77a8-4cb4-42c35be343a1.png)


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

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

# パッケージ読み込み

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

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

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

ダウンロードした基盤地図情報のzipファイルを解凍するとさらにzipファイルが現れます。
![01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/81670bfc-e47d-e7c6-d8bc-a8c01f1aa158.png)

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

```
unzip '*zip'
```
解凍が終わるとxmlファイルが現れます。
![02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/e15aa72a-c925-fa8d-19f5-b53391ad7919.png)


# 数値標高モデルの変換

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

```{r}
#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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/f9cf84b3-bb78-c5dd-9fc5-8b977fbe0904.png)


# 基本項目

続いて基本項目（道路や河川）の変換です。こちらもベースは[「Rで国土地理院 基盤地図情報データを扱う」](https://uribo.hatenablog.com/entry/2018/12/16/084535)を参考にしています。

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

```{r}
# 水域
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でポリゴンを融合する（実践編）](https://misaki-blog.com/190806-r-dissolve-2/)）

以上の作業で出力したtiffやgpkgファイルを読み込んで、QGISで表示できました。
![04.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/60438/0fbc8134-283e-0b09-e50b-26bbcab47fe9.png)