目的
日本のいろんなところにあるポリゴンの面積を計算したいのだけど、場所によって適した投影座標系に合わせなければならないが、いちいちESPGコード変える作業がめんどくさいので、R (とQGIS)を使ってできるだけ自動的に処理する。
(だけどちょっと小さくなった、ナンデ...)
準備
- 日本各地にポリゴンがあるデータ。
- 今回は国土数値情報 世界自然遺産からA28-10_GMLを使っています。
- QGIS 3.4以上
- R version 3.6以上
平面直角座標系とは
ざっくりいうと、平面なら面積を計算しやすいが、地球は丸いので、地球の狭い範囲を部分的に平面に直してから、面積を計算します。この平面に変えた座標系が平面直角座標系で、日本では19の地域に分けられていて、それぞれの地域で当てはめる平面直角座標系が違う。
手順
1. 各ポリゴンにあったESPGコードを追記する。ここではQGISで属性テーブルを編集しました。1,2
-
どの地域がどのEPSGコードに対応するかは、国土地理院のマップや系番号とESPGコードの対応表をみながら入力.
-
国土地理院：１９の平面直角座標系のマップ
-
検索すると、既に作っているサイトもあるよ。ありがたい
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で計算(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 |
なんかぁ微妙ーーーに小さいのね..