経度緯度からその位置の流域を取得したい
タイトルの通り、そういう仕事があった。流域というのは、ある地点に降り注いだ雨がどの河川に流入するのかを示した概念で、たとえばどこかの山に降り注いだ雨が地面を伝い、あるいは染みこみ、最終的にどこの川に流れ込むのかということを示している。そういう河川流域やその水系データは国土交通省がオープンデータとして提供していて、国土をだいたい100メートル四方に区切り（地球は扁球なので実際の辺の長さはちょっとちがう）、その区域がどの流域に所属するのかを網羅している。いわゆるメッシュデータといわれるものだ。しかしこのままではデータ量も膨大で大変使いにくい。データ形式はGMLやShapefileといった形式なのだが、日本の国土全体で10ギガバイトをかるく超える。必要なのは流域だけなので単純なポリゴンであればいい。今回はデータがメッシュである必要はないのだ。
まずはすべてのメッシュデータを取得する。「流域メッシュデータ」などのキーワードで検索するとすぐに国交省のページがヒットするので、ここからZipで圧縮されたメッシュデータをダウンロードしたい。しかし使用目的についてのアンケートが必要だったり、ダウンロードしたいメッシュコードをいちいちクリックして選択させたり、ダウンロードボタンをひとつひとつクリックして、さらに「○○MBあります。本当にダウンロードしますか」的なダイアログが毎回出たりする。いったいどういう発想でこんな使いにくいサイトが用意できるのだろう。オープンデータとして無料で公開してくれるのは本当にありがたいが、手間をかけて不便なサイトを作るくらいなら全てのデータをまとめてどこかの大学でサーブしてもらえば良いんじゃないかと思う…。
面倒なので一度だけアンケートに回答して普通にファイルをダウンロードする。そのときのファイルのURLから他のファイルのURLを予想して適当なスクリプトでまとめて落とした。並列にやりたいのでただのCLIコマンドをたたくだけだがRubyを使った。全部で154の書庫ファイルが手に入る。よく見てみるとファイルによって中身のディレクトリの構造が違っている。面倒だ。
ns = "3036 3623 3624 3641 3653 3741 3841 3927 3928 3942 4027 4028 4040 4042 4142 \
4440 4530 4531 4540 4630 4631 4729 4730 4731 4739 4740 4828 4829 4830 4831 4839 \
4928 4929 4930 4931 4932 4933 4934 4939 5029 5030 5031 5032 5033 5034 5035 5036 \
5038 5039 5129 5130 5131 5132 5133 5134 5135 5136 5137 5138 5139 5229 5231 5232 \
5233 5234 5235 5236 5237 5238 5239 5240 5332 5333 5334 5335 5336 5337 5338 5339 \
5340 5433 5436 5437 5438 5439 5440 5536 5537 5538 5539 5540 5541 5636 5637 5638 \
5639 5640 5641 5738 5739 5740 5741 5839 5840 5841 5939 5940 5941 5942 6039 6040 \
6041 6139 6140 6141 6239 6240 6241 6243 6339 6340 6341 6342 6343 6439 6440 6441 \
6442 6443 6444 6445 6540 6541 6542 6543 6544 6545 6546 6641 6642 6643 6644 6645 \
6646 6647 6741 6742 6747 6748 6840 6841 6842 6847 6848"
Parallel.each(ns.split, in_process: 4) do |n|
url = "http://nlftp.mlit.go.jp/ksj/gml/data/W07/W07-09/W07-09_#{n}-jgd_GML.zip"
puts `wget #{url} -nv -O shape/#{n}.zip`
puts `unzip -q -d shape/ shape/#{n}.zip`
puts `rm shape/#{n}.zip`
end
国がオープンデータとして提供してくれる情報は毎度のことながら大変使いづらい。普通に国土を等分割していたら現れないようなこまかな隙間が存在していたり、日本語の表記がSJISだったりなどのうっとうしい部分はあるが、どれも修正するためのプログラムさえ書ければだいたいなんとかなる。一時期話題になったsyukujitsu.csvや悪名高いken_all.csvに比べればはるかに簡単な前処理ですんだ。
Shapefileはバイナリデータなので、Rubyで加工しやすくするためGeoJSON形式に変換する。gdalというソフトウェアをインストールするとogr2ogrというCLIツールがインストールできるので、ここからの作業はこれを多用する。GISツールはけっこうたくさんの種類があるが、無料で使えるのはツールは限られている。QGISは非常に高機能で処理性能も高いため広く使われており、黒い画面を使いたくない場合はQGISの綺麗なGUIをつかっても問題ない。
Parallel.each(Dir['./shape/**/*.shp'], in_process: 12) do |file|
puts file
puts `ogr2ogr -f geoJSON shape/#{file[/\d{4}/, 0]}.geojson #{file}`
end
指定したディレクトリ以下のすべてのShapefileをGeoJSONに変換して一つのディレクトリにまとめる。本来はここで各メッシュデータのうち、流域コードが同じものをマージしていきたいが、そううまくいかない。メッシュデータの座標処理が甘いのか、各値に若干の誤差がある。するとポリゴンとしては接続されていないと判断されてしまったようでマージされない。変換したGeoJSONファイルの中身をテキストエディタでのぞいてみると、どうやら経度緯度の少数第7位くらいから循環小数を普通に丸めたものと異なる値になっていることがわかる。とりあえず適当に下7桁目以降を丸めてあげるとうまくいった。
ちなみに、各キーの意味は以下のように定義されている。
- W07_001 細分メッシュコード
- W07_002 水系域コード
- W07_003 河川コード
- W07_004 水系名
- W07_005 河川名
- W07_006 単位流域コード
Parallel.each(Dir['./shape/*.geojson'], in_process: 12) do |file|
json = JSON.parse(File.read(file))
json['features'].each do |feature|
feature['geometry']['coordinates'][0].each do |coordinate|
coordinate.map! { |v| v.round(6) }
end
end
File.write("./rounded/#{file[/\d{4}/, 0]}.geojson", JSON.dump(json))
end
これでやっとマージ処理ができるようになる。そしてここが最も時間のかかる部分だ。
Parallel.each(Dir['./shape/*.geojson'], in_process: 12) do |file|
puts file
output = "./dissolved/#{file[/\d{4}/, 0]}.shp"
sql = "SELECT ST_Union(geometry),* FROM 'OGRGeoJSON' GROUP BY W07_003"
puts `ogr2ogr #{output} #{file} OGRGeoJSON -dialect sqlite -sql "#{sql}"`
end
ogr2ogrでST_Unionを使うと同じ辺を共有するポリゴン同士を融合することができる。とにかくこの処理は重たいので一昼夜はかかることを覚悟しよう。QGISでも同じことはできるし、理由はわからないがそっちのほうがだいぶ速い。私は結局QGISのバッチ処理ツールをつかってやった。並列処理はできないが、それでも如実に早く処理が終わる。
QGISを使う場合、154個分のファイルをGUIでいじるのはかなり大変なのでバッチファイルの形式で書き出してあげる。中身はふつうのJSONなので以下のようにすると、GUIをぽちぽちしなくてすむ。デフォルトのJSONエクスポート機能はぶっこわれていて、エクスポートしてもそのファイルをインポートできないのだ。FalseをNoと解釈しないとか、勝手に文字をエスケープしたりだとかするので自分でつくってしまったほうが早い。
list = Dir['./shape/*.geojson'].map do |j|
num = j[/\d{4}/, 0]
({
OUTPUTS: {
OUTPUT: "/Users/xxx/yyy/shape/#{num}.shp"
},
PARAMETERS: {
INPUT: "/Users/xxx/yyy/shape/#{num}.geojson",
DISSOLVE_ALL: 'No',
FIELD: 'W07_003'
}
})
end
File.write('list.json', JSON.dump(list))
メッシュ同士を融合させたGeoJSONファイル154個を一つのファイルにまとめてもう一回融合させる。最初にめっちゃでかいGeoJSONファイルひとつにconcatした状態で融合してもいいんだけど、試してみたらメモリが枯渇してしまった。まぁそんなわけでバラバラのまま加工したが、これも分割統治的なアプローチなんじゃなかろうか。
完成した日本全国の流域データは、各流域が一つのポリゴンとして扱えるので元データに比べると非常に軽量なものとなる。とはいえ数万個の流域があると100MB弱くらいはあるのでブラウザから使うようなものではない。なので汎用的でファイルサイズの小さいShapefileに戻してあげてもいいだろう。
さて、当初の目的は「経度緯度からその位置の流域を取得したい」というものだった。ここまで適当にRubyのスクリプトで元データを処理してきたが、ここからはうまいこと加工後の流域ポリゴンにアクセスできないといけない。具体的に言えば地図上の一点と流域ポリゴンの内外判定を行えばいい。でもそのポリゴンの内外判定とかがめんどくさいのでできるだけあるもので間に合わせたい。そういえばMongoDBにGeoSpatial型を保存できたなと思い出した。調べてみると経度緯度の配列を保存することができ、その内外判定をクエリとして投げられるようだ。早速ためしてみる。
require 'json'
require 'mongo'
client = Mongo::Client.new('mongodb://localhost:27017/basin')
collection = client[:record]
records = []
JSON.parse(File.read('./japan-basin.geojson'))['features'].each do |f|
next if f['geometry']['coordinates'].empty?
att = f['properties']
records.push({
mesh_code: att['W07_001'].to_i,
basin_code: att['W07_002'].to_i,
river_code: att['W07_003'].to_i,
basin_name: att['W07_004'],
river_name: att['W07_005'],
loc: f['geometry']
})
end
collection.insert_many(records)
適当にバルクインサートしてインデックスを張って経度緯度を含めたクエリを投げると数ミリ秒でちゃんと返ってくる。クリックした点から50km以内の流域を表示する、とかのクエリも投げられるのでだいぶ柔軟につかえる。いいかんじ。あとはそれっぽくHTTPリクエストで取得できるようにしておけばとりあえず完成かしら。あ、ここまでぜんぶRubyで書いてるんだからSinatraにすればよかった。まぁ中身は一緒だ。
const app = require('express')();
const MongoClient = require('mongodb').MongoClient;
MongoClient.connect('mongodb://localhost:27017/basin', (err, db) => {
app.get('/', (req, res) => {
const query = {$geoIntersects: {$geometry: {
type: 'Point', coordinates: [parseFloat(req.query.lon), parseFloat(req.query.lat)]
}}};
db.collection('record').find({loc: query}).toArray((err, result) => res.json(result));
});
app.listen(3000, () => console.log('started.'));
});