cucumber flesh

Rを中心としたデータ分析・統計解析らへんの話題をしていくだけ

kuniezu: 日本の国土地理を扱いやすくするRパッケージをCRANに登録しました

development kuniezu

https://repository-images.githubusercontent.com/259226015/c936a500-8e9e-11ea-81e1-8c7635268123

はじめに

kuniezuパッケージ (v0.1.0) をCRANにリリースしました。

github.com

このパッケージは、私が業務や趣味で日本国内の地理空間データを扱う時に作っていた関数を一つのパッケージに整理したものです。 空間的に世界規模のデータを扱うのではなく、日本国内に限った話であれば、日本に即した仕様や座標参照系を利用した方が良いことがあります。 そうした日本の地理空間データを処理する際に利用することがある機能や、あると便利なデータセットを提供できるように努めています。

ゆるゆると開発するつもりでいましたが、Twitterで開発中であることをつぶやいたところ想定以上に反響が大きかったのでCRAN登録を急ぐことにしました。

この記事では、そんなkuniezuパッケージの機能紹介を行います。

まずはみなさんパッケージがインストールされていないと思うので、利用されたい方は以下のコマンドを実行してインストールを済ませてください。

install.packages("kuniezu")

library(kuniezu)
# 使い方を説明するために以下のパッケージも読み込んでおきます
library(sf)
library(ggplot2)
library(leaflet)

使い方

最初のリリースでは以下の機能を実装しました。

  • DMS表記を十進数表記に変換
  • 日本測地系2011における平面直角座標系の特定
  • 南西諸島・小笠原諸島を移動した日本地図の描画
  • 地理院タイルをleafletで簡単に利用できるように
  • 国土地理に関するいくつかのデータセット

順に解説します。

parse_*_dohunbyo(): DMS表記を十進数表記に変換

地球上の任意の座標を示す際に使われる緯度と経度の表記には、十進数で表される場合と「度(Degree)・分(Minute)・秒 (Second)」を用いて表記される場合があります。 後者はDMS表記と呼ばれ、十進数で139.7414、35.6581と表される座標に対してE139°44′28.8869″、N35°39′29.1572″のように表現します。 更に日本では度分秒の記号に対して漢字使って東経139度44分28.8869秒、北緯35度39分29.1572秒等の表記が用いられることがしばしばあります。

この日本独自のDMS表記を十進数に変換する関数がparse_*_dohunbyo()です（アスタリスクの部分にはlon、latがそれぞれ指定できます）。 lonは東経、西経を示す(180~-180）経度、latが北緯、南緯（90~-90）の緯度を扱うために使い分けます。

# 日本経緯度原点（東京都港区麻布台二丁目）の座標
parse_lon_dohunbyo("東経139度44分28.8869秒")

## [1] 139.7414

parse_lat_dohunbyo("北緯35度39分29.1572秒")

## [1] 35.6581

これは過去のブログ記事に書いた関数をパッケージに含めたものです。 加えて、十進数に変換する前段階として漢字を記号に処理するreplace_dohunbyo_kanji()も用意しています。

replace_dohunbyo_kanji("東経139度44分28.8869秒")

## [1] "E139°44’28.8869."

日本測地系2011における平面直角座標系の特定

日本の現行の測地系である日本測地系2011には、19の平面直角座標系の区分がなされています。

平面直角座標系は狭い範囲の測量や距離の算出に用いるのに適しています。 任意の経緯度を平面直角座標に変換するには、まず座標が平面直角座標系のどの系番号に含まれるのかを知っておく必要が生じます。 これを求めるのがst_detect_jgd2011()です。

sfパッケージでおなじみのst_*()が関数名につくように、関数にはsfcオブジェクトを与えて実行します。

# 日本経緯度原点の座標
st_detect_jgd2011(st_sfc(sf::st_point(c(139.7414, 35.6581)), crs = 4326))

## [1] 6677

move_jpn_rs(): 南西諸島の一部・小笠原諸島を移動した日本地図の描画

日本だけを地図上に描画する際、表示する画面の都合からか、南西諸島の一部や小笠原諸島の位置を海上に移動することがあります（例えばNHKの気象情報のページや白地図）。 これを実現するには、予め移動対象の地物に対する操作が必要となります。

move_jpn_rs()はこれらの島々の位置を変更するための関数です。 地図描画のために奄美群島と沖縄県、小笠原諸島の位置を変更する処理が実行されます。 また、この地図化もパッケージの関数で行えるようにggplot2ベースのgeom_jpsegment()を用意しました。

使い方は次のように、まず都道府県や市区町村レベルでポリゴンが分かれた地理空間データに対してmove_jpn_rs()を適用、その後マッピングを行います。

move_jpn_rs(jgd2011_bbox) %>%
  ggplot() +
  geom_sf() +
  geom_jpsegment()

デフォルトの設定で移動した地物がわかるように境界線が引かれます。

この関数の元も、かつて記事に書いたものです。

少し前に @shoei05 さんがブラッシュアップしてくれたのをきっかけにパッケージへ含めました。

一方で、まだこの関数はさらなる磨き上げが必要で、解決策がわかる人はぜひ教えてください。

github.com

地理院タイルをleafletで簡単に利用できるように

インタラクティブなウェブ地図を利用したことがある方ならわかるかと思いますが、leafletでは任意の地図タイルを背景画像として扱うことができます。 国土地理院が発行する地理院タイルもまた、この背景画像として利用可能で、 一覧ページに示されるように多種多様な地図タイルが整備されています。

kuniezuパッケージではこの地理院タイルをleafletで簡単に導入できるようになっています。 これらはgsi_tilesオブジェクトを通して呼び出します。

現在のバージョンでは「近年の災害」をはじめとして、すべての地理院タイルに対応しているわけではありません。 kuniezuで利用可能な地理院タイルの一覧は以下のコードで確認できます。

names(gsi_tiles)

##  [1] "standard"             "pale"                 "english"             
##  [4] "lcm25k_2012"          "lcm25k"               "ccm1"                
##  [7] "ccm2"                 "vbm"                  "vbmd_bm"             
## [10] "vbmd_colorrel"        "vbmd_pm"              "vlcd"                
## [13] "lum200k"              "lake1"                "lakedata"            
## [16] "blank"                "seamlessphoto"        "ortho"               
## [19] "airphoto"             "gazo4"                "gazo3"               
## [22] "gazo2"                "gazo1"                "ort_old10"           
## [25] "ort_USA10"            "ort_riku10"           "relief"              
## [28] "anaglyphmap_color"    "anaglyphmap_gray"     "hillshademap"        
## [31] "earthhillshade"       "slopemap"             "slopezone1map"       
## [34] "afm"                  "lcmfc2"               "lcmfc1"              
## [37] "swale"                "cp"                   "jikizu2015_chijiki_d"
## [40] "jikizu2015_chijiki_i" "jikizu2015_chijiki_f" "jikizu2015_chijiki_h"
## [43] "jikizu2015_chijiki_z" "jikizu_chijikid"      "jikizu_chijikii"     
## [46] "jikizu_chijikif"      "jikizu_chijikih"      "jikizu_chijikiz"

使い方は次のようになります。leafletおよびmapviewでの例を紹介します。

# 出力結果は省略します
gsi_tiles$standard %>%
  addCircles(
    data = sf::st_transform(extreme_points %>%
      purrr::reduce(c),
      crs = 4326))

mapview::mapview(st_sfc(sf::st_point(c(139.7414, 35.6581)), crs = 4326), 
                 map = gsi_tiles$pale)

なお、これらの地理院タイル利用する際はタイルごとに設定されている利用規約に従ってください。

国土地理に関するデータセット

これまでの例で出てきたものがありますが、日本の国土地理を扱う上であると便利なデータ等をパッケージで提供するようにしました。 現在、次のようなデータがあります。

  • extreme_points… 東西南北端点の経度緯度データ
  • jgd2011_bbox… 日本測地系2011の平面直角座標系の区域を示すポリゴンデータ
  • jp47prefectural_offices… 47都道府県の庁舎の位置データ

こんなデータや機能が欲しいという方はぜひ issue に書き込んでください。 バグ報告もありがたいです。 GitHub以外ではTwitterでコメントいただいても対応します。

Enjoy!