2015年7月21日火曜日

日本の海岸線データを手にいれる

レポートを書くとき、何かと使うことの多いのが地図データです。しかも緯度経度で海岸線や県境のデータがあるとベストです。(緯度経度で地点をマークして一緒にEXCELで散布図にしてやれば正確な表現ができます)
実は国土地理院では昔から公開していますが、数年前からやたら複雑な書式のxmlなので尻込みしていました。少し長い休みがあるので、重い腰をあげてデータ整理でもしてみましょう。

1.海岸線データを入手する
ぐぐれば簡単に見つかります。DLするとき何に使う?とか簡単なアンケートに答えないといけませんが、「趣味」の項目もありますから特に気張らずに。(ただ、ライセンス的にフリーというわけではないので注意はしてください。個人目的、研究用等に使うならご自由にという程度です)
データは県ごとに分かれています。

2.xmlフォーマット
各県ごとにデータが分かれています。しかも、県の中でも行政界ごとにxmlのタグが分けられているようです。以下に例を示します。

  1. <gml:LineStringSegment>
  2. <gml:posList>
  3. 35.55061306 140.11532166
  4. 35.55081833 140.11538084
  5. 35.55206806 140.11517666

xmlのタグ、<gml:posList>で識別されています。数字は一目で緯度経度とわかります。(国土地理院なので、今は当然世界測地系です)

3.海岸線データを取り出す
テキストエディタで取り出してもいいんですが、なにせ数があります。こちらは単に日本地図が描きたいだけなんですが、きちんと県別、さらに行政界毎にデータがタグで分離されているため、テキストエディタで手作業でやってたら死にます。
とりあえず、xmlなので、以下のプログラムでこのタグ部分だけ抽出してやります。

  1. # coding: utf-8
  2.  
  3. import xml.dom.minidom
  4.  
  5. if __name__ == "__main__":
  6.  
  7. dom = xml.dom.minidom.parse("sample.xml")
  8.  
  9. for line in dom.getElementsByTagName("gml:posList"):
  10. print line.firstChild.data

以下の様な感じで海岸線のデータだけ標準出力にでてきますので、それをファイルにリダイレクトしてやります。

35.55061306 140.11532166
35.55081833 140.11538084
35.55206806 140.11517666
35.55243944 140.11512389
35.55278833 140.11507416

4.データの間引き
さてこうしてできた海岸線のデータですが、単に数値の羅列ですがなんと100MB強あります!県毎にsheetを分ければなんとか日本の海岸線を散布図で描けますが、一つのファイルにまとめてしまうと、サイズが(行数が)大きくEXCELが扱えないといってきます。どうせそんな細かいところまで見ないので、データの間引きをします。そのために用意したプログラムを以下に示します。

  1. #!/usr/bin/python
  2. # coding: utf-8
  3.  
  4. import sys
  5.  
  6. if __name__ == "__main__":
  7. argv = sys.argv # コマンドライン引数を格納したリストの取得
  8. argc = len(argv) # 引数の個数
  9.  
  10. # print argv
  11. # print argc
  12.  
  13. fi = open(argv[1], 'r')
  14. fo = open(argv[2], 'w')
  15.  
  16. # n行に1回出力するようにする
  17. ct = 0
  18.  
  19. for line in fi:
  20. if line.strip() == "": # 島等のデータとの境界である空白行は残すようにする
  21. fo.write(line)
  22. ct = 0 # 島のデータは少ないので1行でも残す
  23. continue
  24.  
  25. if ct == 0:
  26. fo.write(line)
  27.  
  28. ct = ct + 1
  29.  
  30. if ct == 10: # 10行に1回出力する
  31. ct = 0
  32.  
  33. fi.close()
  34. fo.close()

今回は起動時に入力ファイル、出力ファイルを読むようにしました。あと間引く時、元のデータでは行政界毎、及び島毎に空行が入っていました。ここを無視して間引いいてしまうと、小さな島だとデータがなくなってしまいかねません。ちょっとだけ工夫しました。間引く数は実際に試して実験です。結局、10分の1(約10MB)にしないとまともにEXCELは動きませんでした。使うPC環境にもよりますが、もう少し間引いた方がいい気がします。またどうせかなりの縮尺で使うつもりなので、こだわって島のでデータを消さないようにしました。小さな島は消したいんですが、属性として地名は入っているんですが、それが島かどうかは人間が判断してやらないといけません。(GISデータは最後は人力で編集しないと、いいデータにはなりませんね(>_<))

5.完成
以下に実際に作成した日本(白)地図を示します。


使うときは例で示したように、別途示したい位置の緯度経度のデータを用意して、データ系列を追加してやります。注意が必要なのは、普通にやると追加したデータのグラフも散布図の折れ線になってしまいます。単にX-Yの座標だけの散布図にしたい場合が多いでしょうから、その時はデータ系列を追加してから、追加したデータを選択し、グラフの種類を変えてやります。(「複合グラフ」といいます)


0 件のコメント:

コメントを投稿