日本の海岸線データの加工

かなり国土地理院のデータを間引いて、なんとか100MBのデータを数MBまでにしましたが、ノートPCなんかで使うにはまだまだサイズが大きいです。思い切って1/100に間引いてみましたが、日本の海岸線が途切れてしまいます。
やはり少し頭を使って、元になる国土地理院の海岸線データの特徴を調べてみます。

1.海岸線データの固まり方
先にも言いましたが、元データは県別のxmlデータになっていて、行政界、島毎にタグが分かれています。それを緯度経度だけ抽出してくると以下ののようになります。

35.04924833 136.83906528
35.05145167 136.83815250
35.05175111 136.83802834
   

35.05899417 136.87838278
35.05921944 136.87846166
35.05952000 136.87857805
35.05949694 136.87866195
35.05947139 136.87874555
35.05914306 136.87861805
35.05894917 136.87855139
35.05899417 136.87838278
   

35.05730000 136.87774750
35.05771861 136.87788778

途中、２行の空白行がありますが、そこが行政界の区切り、あるいは島のデータの始まりを表します。上記サンプルを見ると、２つ目のデータの固まりのデータの始まりと終わり（6行目と13行目です）が同じ緯度経度になっています。つまりこれは島だということです。
前回のプログラムでは、以下に注意しました。

• 行政界、島の区切りを示す空白行は残す。
• 空白行を出力したところで間引きカウンタを０クリアして、空白行に続く先頭データは積極的に出力する。

しかし、特に島のデータに顕著にでましたが、あまり間引く数が大きいときはデータの固まりの最後（空白行の直前のデータ）も残してやらないと島が閉曲線になりません。

2,改良版プログラム
ということで、島の最後のデータも残してやるように改修しました。

#!/usr/bin/python
# coding: utf-8

import sys

if __name__ == "__main__":
 argv = sys.argv  # コマンドライン引数を格納したリストの取得
 argc = len(argv)  # 引数の個数

# print argv
# print argc

 fi = open(argv[1], 'r')
 fo = open(argv[2], 'w')

 # n行に1回出力するようにする
 ct = 0

 preLine = ""

 for line in fi:
  if line.strip() == "": # 島等のデータとの境界である空白行は残すようにする
   if (preLine.strip() != ""): # 空白行の直前のデータは残す（島の閉曲線を守るため）
    fo.write(preLine)
   fo.write(line)
   ct = 0 # 島のデータは少ないので１行でも残す
   preLine = line
   continue

  if ct == 0:
   fo.write(line)

  ct = ct + 1
  preLine = line

  if ct == 100:  # 100行に１回出力する
   ct = 0

 fi.close()
 fo.close()

さてこれで1/100に間引いてみましたが、どうもデータのサイズが思うように小さくなりません。

3.改良版の処理結果を見る
再度、処理結果を調べてみました。

35.08614861 136.89070723
35.08342250 136.86787278
35.05776194 136.84691778
35.06403806 136.84048584
   

35.05175111 136.83802834
35.05175111 136.83802834
   

35.05899417 136.87838278
35.05899417 136.87838278

確かに島のデータは残されるようになっていますが、思いの外国土地理院のデータには小島がしっかりはいっているようで、同じ緯度経度が２回連続する固まりが大量にありました。こいつらは改良版のプログラムではどんなに間引率を大きくしても残されるようにしてあるため、思ったようにデータのサイズが小さくならなかったようです。
でも実際日本地図を使うようなときはこんな小島、表示もされないのでデータの無駄です。

4.小島のデータを削除する
手作業で削除するかとも思ったんですが、すぐに諦めました。さすがにデータのサイズが小さくならない原因だけあって大量に存在しています。あきらめて、２行同じ緯度経度のデータが続く箇所はフィルタするようなプログラムを作りました。（行政界の他のところで偶然そういう場所がないか、という恐れはありますがまあ少しくらい抜けてもそれに気づくほど拡大して使うことはないでしょう。もし使う必要があるなら、それはもう手作業でなんとかしてやるしかありません）

#!/usr/bin/python
# coding: utf-8

import sys

if __name__ == "__main__":
 argv = sys.argv  # コマンドライン引数を格納したリストの取得
 argc = len(argv)  # 引数の個数

# print argv
# print argc

 fi = open(argv[1], 'r')
 fo = open(argv[2], 'w')

 # n行に1回出力するようにする
 ct = 0

 preLine = ""

 for line in fi:
  if line.strip() == "": # 島等のデータとの境界である空白行は残すようにする
   if (preLine.strip() != ""): # 空白行の直前のデータは残す（島の閉曲線を守るため）
    fo.write(preLine)
   fo.write(line)
   ct = 0 # 島のデータは少ないので１行でも残す
   preLine = line
   continue

  if preLine != line:
   fo.write(preLine)
  else:
   preLine = ""
   line = ""
   continue

  preLine = line


 fi.close()
 fo.close()

かなり手抜きプログラムですが、これでサイズを1/100にし、かつ２行しかないような小島のデータを削除して、データサイズが1.4MBまで縮小できました。結果を以下に示します。

先日の海岸線とほとんど見た目は変わりません。よく見れば、小さな小島は消えているんですがね。

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

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

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

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

  <gml:LineStringSegment>
   <gml:posList>
35.55061306 140.11532166
35.55081833 140.11538084
35.55206806 140.11517666

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

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

# coding: utf-8

import xml.dom.minidom

if __name__ == "__main__":

 dom = xml.dom.minidom.parse("sample.xml")

 for line in dom.getElementsByTagName("gml:posList"):
  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が扱えないといってきます。どうせそんな細かいところまで見ないので、データの間引きをします。そのために用意したプログラムを以下に示します。

#!/usr/bin/python
# coding: utf-8

import sys

if __name__ == "__main__":
 argv = sys.argv  # コマンドライン引数を格納したリストの取得
 argc = len(argv)  # 引数の個数

# print argv
# print argc

 fi = open(argv[1], 'r')
 fo = open(argv[2], 'w')

 # n行に1回出力するようにする
 ct = 0

 for line in fi:
  if line.strip() == "": # 島等のデータとの境界である空白行は残すようにする
   fo.write(line)
   ct = 0 # 島のデータは少ないので１行でも残す
   continue

  if ct == 0:
   fo.write(line)

  ct = ct + 1

  if ct == 10:  # 10行に１回出力する
   ct = 0

 fi.close()
 fo.close()

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

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

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