はじめに
動機
最近、週末にいろいろなデータをベクトルタイルにしてウェブ地図にプロットすることを試しています。先日の国土数値情報に続き、今回はFITの事業計画認定情報を可能な範囲でウェブ地図化するためのベクトルタイルにしてみました。自分のためのメモとして残す意味あいが大きいですが、記録を残しておきます。
私の環境
WindowsPCですが、使ったツールは以下のようなものです。
- Windows 11 Home
- WSL2（Ubuntu）
- tippecanoe/felt
- Windows PowerShell
- R (4.3.2)
- メモ帳（Windows メモ帳 11.2405.13.0）
- Excel
- QGIS ver. 3.34.2
- WSL2（Ubuntu）
手順
Step 1: データダウンロード
再生可能エネルギー 事業計画認定情報のページから都道府県別のエクセルシートをダウンロードします。個別にクリックしてダウンロードするしかないですが47都道府県をダウンロードしても10分もかかりませんでした。
Step 2: データ加工・CSVデータ作成
エクセルファイルだと扱いにくいので、CSVファイルにします。こちらも地味にExcelで作業します。
今回は、「位置」に注目したいので二つあるシートのうち、「すべての場所所在」シートのデータを使うことにしました。（ただし認定日や稼働日がないので、それらを見たいときはもう一つのデータのほうが良いかもしれません。）
「すべての場所所在」のデータをコピーします。
新しいシートにデータだけをはりつけCSVファイル（UTF8ではない）として保存します。この辺りも省力化したいのですが、とりあえず手作業でやります。（開く→タブを変える→コピー→新しいシート→貼り付け→ラベル削除→CSVで保存 の流れを繰り返します。）
こうして47個のCSVファイルができます。
Step 3: 東大のCSVアドレスマッチングサービスを使って緯度経度を付与する
CSVファイルの住所は文字列なので、住所の文字列から緯度経度に変換する必要があります。これをジオコーディングというのですが、東京大学がCSVアドレスマッチングサービスというものを提供しています。利用規約は以下の通りです。
このCSVアドレスマッチングサービスでは、こちらで変数を調整してCSVをアップすると処理されたデータがダウンロードできるようになっています。「全国街区レベル（経緯度・世界測地系）」にして、住所のカラムは7番目など、各種パラメータを指定してファイルを選び送信を押していきます。
処理が終わったファイルを見てみると、分割できた住所・地名、緯度経度、変換の信頼度（iConf）、変換された住所レベル(iLvl)もついています。これは助かります。
こうして47の緯度経度座標を持つファイルができました。
Step 4: 47ファイルをまとめ場所の重複排除
まず47ファイルをまとめます。もう少し早いタイミングでやってもいいのですがCSVアドレスマッチングの処理を軽くするために前段までは都道府県ごとのファイルでやっていました。ファイル統合にはWindowsPowerShellを使います。各CSVの最初の行はラベルなので、1行目を飛ばしてマージします。
$files = Get-Childitem -Path 2_geocoded/*.csv
foreach($file in $files){cat $file | Select-Object -skip 1 | Add-Content 3_merged.csv}
このくらいのサイズになりました。
このデータですが、Rを使って件数などを見てみました。以下の通りでした。
ジオコーディングでプロットしてきましたが、同一地点に配置された同じIDのものは削除するようにします。GISソフトでの操作を念頭に書き出しはUTF8にしましたが、エラー等はありませんでした。
# FITの作業用
# June 23, 2024
#
#ライブラリの読み込み
library("dplyr")
#ファイルの読み込み
fit <- read.csv("3_merged.csv",fileEncoding="cp932",header=FALSE,colClasses=rep("character",13))
#件数
nrow(fit) #809554
#設備ID数
length(unique(fit[,1])) #410393
#発電事業者数
length(unique(fit[,2])) #181400
#発電代表者数
length(unique(fit[,4])) #164395
#発電施設の所在地数
length(unique(fit[,7])) #723880
#ジオコーディングした場所
length(unique(fit[,9])) #164395
#同一地点プットで同一IDのものはまとめる
fit_sel <- select(fit,1,2,4,5,6,8,9,10,11,12,13)
fit_sel <- unique(fit_sel)
nrow(fit_sel)
#書き出し（utf-8で）
write.csv(fit_sel, "4_geocoded.csv",fileEncoding="utf8")
(追記)2024.6.25
上記の方法だと、例えば町字にプロットされる点などは同じ場所に落ちてしまいウェブ地図でも重なってしまうため、少し工夫して場所をずらします。Rで以下のように処理しました。
#ライブラリの読み込み
library("dplyr")
#ファイルの読み込み
fit <- read.csv("3_merged.csv",fileEncoding="cp932",header=FALSE,colClasses=rep("character",13))
#同一地点プットで同一IDのものはまとめる
fit_sel <- select(fit,1,2,4,5,6,8,9,10,11,12,13)
fit_sel <- unique(fit_sel)
nrow(fit_sel) # 437124
###
#同一ポイントに落ちるものを半径radの円上に散らす
rad <- 0.00833333 * 15/1000 #0.00833333 eq 1km
length(unique(fit_sel[,7])) #164395
basho <- unique(fit_sel[,7])
if(exists("result")){rm(result)}
for(i in 1:length(basho)){
sameloc <- filter(fit_sel, V9 == basho[i])
for(j in 1:nrow(sameloc)){
if(exists("out")){rm(out)}
out <- sameloc[j,]
if(exists("kaku")){rm(kaku)}
kaku <- 2*j/nrow(sameloc)
if(nrow(sameloc) != 1){
out[,8] <- as.character(as.numeric(out[,8]) + rad*sinpi(kaku) )
out[,9] <- as.character(as.numeric(out[,9]) + rad*cospi(kaku) )
}
#resultに書き出し
if(!exists("result")){result <- out}else{ result <- rbind(result,out)}
}
#print(basho[i])
}
#書き出し（utf-8で）
write.csv(result, "4_geocoded-2(modified).csv",fileEncoding="utf8")
Step 5: GISデータ化
QGISに緯度経度付きのCSVファイルを読み込ませます。
GeoCodingのときに座標値が取れなかったのがいくつかあるようで、緯度経度0、0のところにもデータがみられますね。
そしてレイヤを右クリック→エキスポートからGeoJSONファイルとして出力します。SHP形式でもいいですが、私はベクトルタイル化するのでGeoJSONにしておきます。
拡大してみてもある程度の場所には来ているように思います。
（追伸）
同じところに落ちる点を周囲に散らばせる処理もまあまあうまくいっています。
点が多すぎるところはちょっとむずかしいですね。。。
Step 6: ベクトルタイル化
ベクトルタイル変換ツールfelt/tippecanoeを使うためにWSL2を使います。（変換ツールや関連ソフトもインストール済みとします。）
GeoJSONファイルは既にあるので、あとは単純に変換するだけです。ズームレベルは13～16にしておきます。
tippecanoe -o fit202405.pmtiles --no-tile-compression --no-feature-limit --no-tile-size-limit --drop-rate=1 -Z13 -z16 -L fit:5_fit.geojson
少し待ちますが、300MB超くらいのPMTilesファイルができました。このサイズ感ならズームレベルを15までにすればpbf形式でも十分対応できそうなサイズですね。
まとめ
FIT事業計画認定情報について47都道府県のエクセルファイルからCSVファイルを作り、ジオコーディングしてベクトルタイルにすることができました。
参考