目的

日本のいろんなところにあるポリゴンの面積を計算したいのだけど、場所によって適した投影座標系に合わせなければならないが、いちいちESPGコード変える作業がめんどくさいので、R (とQGIS)を使ってできるだけ自動的に処理する。
(だけどちょっと小さくなった、ナンデ...)

準備

日本各地にポリゴンがあるデータ。
今回は国土数値情報世界自然遺産からA28-10_GMLを使っています。

QGIS 3.4以上
R version 3.6以上

平面直角座標系とは

　ざっくりいうと、平面なら面積を計算しやすいが、地球は丸いので、地球の狭い範囲を部分的に平面に直してから、面積を計算します。この平面に変えた座標系が平面直角座標系で、日本では19の地域に分けられていて、それぞれの地域で当てはめる平面直角座標系が違う。

手順

　1. 各ポリゴンにあったESPGコードを追記する。ここではQGISで属性テーブルを編集しました。¹,²

どの地域がどのEPSGコードに対応するかは、国土地理院のマップや系番号とESPGコードの対応表をみながら入力.
国土地理院：１９の平面直角座標系のマップ
検索すると、既に作っているサイトもあるよ。ありがたい
- http://blog.godo-tys.jp/2012/11/21/999/
- https://qiita.com/kudarisenmon/items/2003b982247448affcf6

　2. Rで面積計算をする。

epsgMENSEKI.r

# パッケージ読み込み
library(sf) #RでGISをするときは必須
library(dplyr)  #%>%(バイプ)を扱うため

setwd("SHPファイルが保存されたフォルダ") 
Shp1 <- read_sf("世界自然遺産.shp") #shpファイル読み込み
Df1 <- Shp1 %>% as.data.frame() #shpファイルのデータをデータフレームに変換
Df1$area <- 0 #データフレームに面積の列を追加しておく[中身はすべて0]

# i番目のポリゴン→座標系変換→面積計算→データフレームに追記
for(i in 1:nrow(Df1)){ #for文で一括処理
   tmp1 <- Shp1[i,] %>% st_set_crs(4612) %>% st_transform(Df1$epsg[i])
           #元のshpファイルはJSD2000で作成されているのでepsg:4612
   Df1$area[i] <- as.numeric(st_area(tmp1)) #面積計算してデータフレーム書き換え
 }
# 結果をCSVなどに出力すると良い

## 結果表示##
Df1
# 面積の単位はm^2
  A28_001    sekai keibango epsg                      geometry      area
1      01    知床       13 2455 POLYGON ((145.0242 44.05915... 710356000
2      02    白神       10 2452 POLYGON ((140.1686 40.496, ... 168888939
3      03  屋久島        2 2444 POLYGON ((130.5675 30.2487,... 109263115

計算結果の確認

QGISで各ポリゴンごとに計算した結果

Rの計算結果を、QGISで面積計算した結果³と地理院地図で計測した結果⁴で比較してみる。
QGISと地理院地図は同じ値。それらより、Rで計算した結果は0.01%未満で小さい...

	Rで計算(m^2)	QGISで計算(m^2)	地理院地図計測(km^2)	R/QGIS
知床	710356000	710397478	710.397	0.9999416
白神	168888939	168907446	168.907	0.9998904
屋久島	109263115	109278285	109.278	0.9998612

なんかぁ微妙ーーーに小さいのね..

dbfファイルを開いて直接書き込んでも良い。 ↩
数が多いときは地名や系番号とEPSGコードの対応一覧表を別ファイルで作って、R上で呼び出しながら使う方が便利。 ↩
QGISで各ポリゴンを空間参照系を変えて保存して1つ1つ面積計算しました。 ↩
kmlで保存してから地理院地図に読み込んで、ツール>計測で面積を選ぶ。 ↩