ラベル GIS の投稿を表示しています。 すべての投稿を表示
ラベル GIS の投稿を表示しています。 すべての投稿を表示

2017年4月19日

気象庁XML電文を「正しく」画像化するために必要な地図の話

このエントリーをブックマークに追加 このエントリーを含むはてなブックマーク

皆さんは「地図」気にしていますか？
伊能忠敬が日本全国を歩いて測量してまわり、精巧な日本地図を作成してから約200年、今や誰もが手のひらの上で正確な日本地図・世界地図を自由な縮尺で見ることができる便利な時代です。

それにしても、久しぶりの投稿ですが...
今回は、気象庁が提供しているXML電文を画像化する、その時に使う地図の話をしようと思います。

以前、地震発生時の震源・震度情報をブラウザ上でグリグリ動かして閲覧可能なWebサービス「EVI 地震火山詳報」を作りました。
(参照: 地震とか火山噴火情報を閲覧できるWebサイトを作った)

このとき使用した地図は、国土交通省が提供している国土数値情報 行政区域というもので、行政区（雑に言うと市区町村）単位で日本を分割した地図です。
地震が発生した後に気象庁から発表される「震度速報」「震源・震度に関する情報」では、日本の行政区単位で震度が発表されます。

さて、気象庁防災情報XMLフォーマット電文では、地震以外にも多数の情報を配信しています。
例えば...
  • 気象特別警報・警報・注意報
  • 土砂災害警戒情報
  • 竜巻注意情報
  • 地方気象情報
  • 府県気象情報
  • 全般週間天気予報
  • 噴火に関する火山観測報
  • 台風解析・予報情報
などなど...
紫外線やスモッグ、生物季節観測では桜の開花や満開などの情報も配信されています。
こういった情報の多くは行政区単位ではなく、気象庁が定める4つの地域の分け方単位で発表されています。

気象庁｜予報用語　特別警報・警報・注意報や天気予報の発表区域

簡単にまとめるとこんな分類です。下になるほど細かく分割されています。
つまり気象庁は、一番細かい単位として基本的に「二次細分区域」という地域の分け方で情報を発表していることになります。
この二次細分区域ですが、実は行政区域と1対1ではありません。
日本の市区町村数は1741ですが、二次細分区域は1775の地域からなっています。(2017年4月 現在)
ということは、行政区を更に何らかの基準で分割しているということになります。

ここまで、気象庁の用いる発表区域について簡単に説明しましたが、ここで少し気象庁防災情報XMLフォーマット電文を利用するうえでの話をしましょう。
気象庁防災情報XMLフォーマット電文(以下、気象庁XML電文)で配信される情報は、普段私たちがテレビやラジオ、新聞、Yahoo!天気やウェザーニューズなどで見る情報と根本は同じです。
私たちは、気象庁XMLを受信して処理することで、テレビなどのマスコミが発信している気象情報と同じ物を世の中に発信することができます。

ですが、これらの情報は、当然のことですが勝手に改変してはいけません。
例えば、とある市に住んでる人が「気象庁の発表する警報は大げさだ。うちの市は注意報くらいに修正しておこう。」などと勝手に書き換えて、この情報をインターネットに流してはいけません。
以下は、気象庁が公開している、気象庁XMLを利用するうえでの留意事項の抜粋です。
警報の取り扱いについて
警報は重要な情報であり、万が一、誤った警報事項や錯誤を生じさせる情報が流通した場合、気象業務法第23条へ抵触する可能性があるほか、社会への影響が大きいことから、元の電文の本質を損なうような編集は認められません。
公開XML電文の編集・加工について
気象庁が発表した予報内容と異なる独自の予報を発表することは予報業務に該当し、気象業務法により許可を受けた者しか行えません。
気象庁ホームページを通じて公開するXML形式電文のご利用にあたっての留意事項より引用
http://xml.kishou.go.jp/open_trial/considerationforxml.pdf (PDF)

さてさて、そこで地図の話です。
気象庁が発表する情報は、気象庁の定める発表区域で出る訳ですが、もし気象庁の情報を国交省国土数値情報の行政区域を使って塗ったらどうなるでしょうか？

これは、とある会社がTwitterや独自のアプリ内で配信している気象警報・注意報の画像です。
そして、こちらはゲヒルン株式会社が許諾して特務機関NERVで配信している気象警報・注意報の画像です。

注目して頂きたいのはここです。


和歌山県の田辺市は、行政区域としては画像上段のようになっていますが、気象庁の二次細分区域としては画像下段のように5つの地域に別れています。
ここで、2017年 4月18日 3:23に発表された気象警報・注意報を確認すると、田辺市龍神には大雨警報が発表されていますが、田辺市の他の地域には発表されていないことが分かります。

(画像が横に長いので、拡大してご覧ください)

別の地域でも見てみましょう。

以下は北海道の渡島・檜山(おしま・ひやま)地方に発表された気象警報・注意報です。

もう見ただけでお分かり頂けると思いますが...


北海道二海郡八雲町は、気象庁二次細分区域では八雲町八雲と八雲町熊石に分割されています。

そもそも、分割された地域は、市町村等をまとめた地域でも同一の地域に属していません。


更にその上の一次細分区域でも同一地域ではありません。


ここで、2017年 4月18日 8:13に発表された気象警報・注意報を確認すると、八雲町熊石には暴風警報が発表されていますが、八雲町八雲には発表されていないことが分かります。


さてさて...
ここまで例を挙げながら解説してきましたが、いかがでしたでしょうか。
気象警報が出ていない地域に警報が出ている、あるいはその逆という情報を、あたかも真実のように世の中に配信してしまう、そんなことはあってはいけません。

ゲヒルン株式会社では、正しい情報を迅速かつ確実に分かりやすく伝えるために、気象庁二次細分区域の地図データセットやレンダリングエンジンを開発しています。

実は、気象庁二次細分区域よりも更に細かい地域分けとして、土砂災害警戒情報で使用する地域という発表区域もありますが、これに対応する地図データセットも開発しています。


ふう、だいぶ長々と書いてしまいましたが...
こういうシステムを作っている人たちの気持ちは「一人でも多くの人に、安全で安心できる生活を送ってもらいたい」だと思います。
そのためにも、現在気象庁XMLを受信して活用している方も、これから利用して何か作ろうと思っている方も、正確な防災情報を配信するように心がけましょう！

Time: 22:42 Tag: , ,
0 Comments

2016年10月15日

e-statのshpを全国分取り込んで国交省の行政区域shpと比べてみる

このエントリーをブックマークに追加 このエントリーを含むはてなブックマーク

久々にGISネタです。

国交省の国土数値情報ダウンロードサービスでダウンロードできる行政区域のShapefileは以前の記事でも取り扱いました。
国土数値情報 行政区域データは全国の行政界を記録したものですが、これは市区町村界レベルであって、大字(〜丁目)レベルの境界ではありません。(行政界なので当然なのですが)
そこで、現状無料で手に入る大字レベルのshpとして、総務省統計局のデータを利用してみることに。

e-Stat 政府統計の総合窓口

地図で見る統計(統計GIS) → データダウンロード → 国勢調査 → 平成２２年国勢調査（小地域） 2010/10/01 → 男女別人口総数及び世帯総数 → 市区町村を選択 → 世界測地系平面直角座標系・Shape形式

これをダウンロードします。
ダウンロードするのですが... これを全国分となると、市区町村を選択してShapeを選んで...という作業をおよそ1900回繰り返さなければなりません。
一括でダウンロード出来ない理由は負荷軽減のためで、全国分欲しい人は地図屋からデータ買ってねという記載があるので、ちゃんと正攻法で買いましょう。
そもそも、このShapefileは数値情報を上に重ねたりグラフ化したりするためのもので、境界域を得る為に使用するものではありません。

...ただ、今回はスクリプト書いて全部ダウンロードさせてもらいました。(ごめんなさい)
一括ダウンロードのためのスクリプトは載せません。

結論としては、やはり地図屋から買った方がいいです。
しかし、同じことを企む人たちのために、何故そういう結論に至ったのか書き残そうと思います。

まずは、視覚的に比較するためにShapefileを取り込みます。
e-Statから手に入れたデータはzip化されてますので、まずは解凍。

unzip \*.zip

出てきたShapefileをQGISで結合します

1. ベクタ → データマネジメントツール → 複数のシェープファイルを1つに結合する
2. フォルダのレイヤによって選択する にチェック
3. 入力ファイルで .shpファイルを全て選択して、エンコードにSift_JISを選択 → Open
4. 出力シェープファイルで保存先を選んでエンコードにUTF-8を選択 → Save
5. OK で処理開始

この段階で結合されたshpができあがるのですが、CRS(測地系)がEPSG2454(JGD2000)になっていると思うので、今後扱いやすいようにEPSG4326(WGS84)で保存し直します。

1. 結合処理が終了して、レイヤに追加されたデータを右クリックして名前を付けて保存
2. パスに保存先を選び、CRSにEPSG4326(WGS84)を選択してOK

これで準備完了。


千代田区の大手町1丁目を選択してみました。
なんだか良さそうな雰囲気。

では、すこしズームアウトして湾岸部を見ましょう。


なんだか海岸線とは思えない鋭利なラインが目立ちます。
この上に国交省 行政区域データを重ねてみます。


どうやら、海域の一部も含まれているので、陸地の境界域として利用するにはイマイチでしょう。
(こういう海域が取り込まれている部分が全国の海岸部分に沢山あるので、1つ1つ手で除去とかアホみたいに時間掛かることはしません)

それと、もう1点。
これは統計局のデータではなく国交省のデータの問題なのですが...


中央の濃い赤は統計局shpの皇居、周囲の赤は国交省shpの千代田区です。
周囲の境界をよーく見てみると...


濃い赤の選択部分は統計局shpの八重洲1丁目、周囲の赤は国交省shpの千代田区です。
八重洲1丁目が千代田区と中央区に跨がっています。
こういうポイントが全国各地に山ほどあります。

どちらのデータに問題があるのか、原典として使われている地図を調べてみないと分からないですが、おそらく国交省行政区域データの方が粗いのではないかと思われます。
ベクタなので、点の数が増えれば複雑な形状にもフィットする境界を描けて、逆に点の数が少ないと簡素な境界しか描けません。
もちろん、点を増やせばデータサイズも増えます。

ということで、結論ですが...

1. 大字レベルのShapefileが必要なら地図屋さんから買いましょう。
2. 厳密な市区町村界が欲しいなら、購入した大字レベルのShapefileから市区町村レベルで地物の結合をして生成しましょう。

以上です。

普通、業務でこういったGISデータを利用されている方にとっては当然なのでしょうが、個人であれこれやるとぶつかりそうな問題だなぁと思ったので書いておきます。

Time: 8:18 Tag:
1 Comments

2015年12月8日

地震とか火山噴火情報を閲覧できるWebサイトを作った

このエントリーをブックマークに追加 このエントリーを含むはてなブックマーク
作るのを思い立ってから半年くらい掛かった... (実質作業時間は短い)

地震火山詳報
http://evi.prioris.jp

地震と火山(桜島に限らない)に関する気象庁XML電文情報と、それに付随する情報を配信するWebサイトです。
これの立ち上げに伴い、桜島詳報を閉じました。
現時点では火山噴火情報は未実装ですが、統合するつもりで地震火山詳報とネーミングしたので、近々実装します。
と言っても、桜島詳報のコードをちょっと改造して合体させるだけなのですが、最近全国的に火山活動が落ち着いているのでモチベーションが...
→ 実装しました！
12月12日より、噴火に関する火山観測報も追加しました。
桜島詳報では桜島の噴火・爆発のみを扱っていましたが、地震火山詳報では桜島以外の火山(阿蘇山など)の情報を配信されます。

発生した地震の詳細ページでは、震源と揺れを観測した地域の震度がGoogle Mapsライクな操作感で閲覧できます。
例 : 震源・震度情報(2015年12月4日 5時40分頃に発生した地震) - 地震火山詳報

震源や震度のアイコンをクリックすると、上記画像のように地名や座標が表示されます。
地図はOpenStreetMapを使っています。
あとレスポンシブ対応なのでスマホでもPCでも良い感じに閲覧可能です。

地震火山詳報は、地震に関する情報をTwitterとPushbulletでも配信しています。
(Pushbulletって何？っていう人は、EngadgetとかLifehackerの記事を読むといいです。)

Twitter : @Prioris_EVI
Pushbullet : #earthquake_jp

上記Pushbulletのチャンネルは、日本全国のすべての地震(震度1〜)をプッシュ通知で教えてくれるので若干鬱陶しいです。
(日本全国でこんな頻度で地震が起きてるんだ〜！っと認識するのには最高ですが)
というわけで、震度3以上の地震があった時だけプッシュ通知してくれるチャンネルもあります。

Pushbullet : #earthquake_int3over

桜島、阿蘇山の噴火・爆発時プッシュ通知はYoにて配信しています。

Yo : SAKURAJIMA (鹿児島県 桜島の噴火・爆発時に配信)
Yo : ASOSAN (熊本県 阿蘇山の噴火・爆発時に配信)

技術的な話


以下の記事の集大成みたいな感じです。

PostgreSQLとPostGISで国土数値情報(行政区域)を扱ってみる (2015/06)
PostGISで気象庁の細分区域&市町村等に対応する重心を求める (2015/06)
OpenLayers 3で気象庁発表の震度をマッピングする (2015/07)
OL3で扱える地図タイルをShapeから生成する (2015/12)

あと、ベースとして以下も。
桜島が噴火したYoの裏側 (2014/08)

上記の記事中にも書いていますが、気象庁XML電文から震度マップ(GeoJSON)を作るのに必要なデータやGeoJSON生成コードはオープンソースとして公開しています。

9SQ/jma-eqarea-centroid (気象庁が用いる地域区分の重心座標テーブル)
9SQ/jma-eqxml2geojson (気象庁XMLから上記テーブルを利用してGeoJSONを生成する)
9SQ/seismic-intensity-map (震度GeoJSONをOpenStreetMapにオーバーレイ表示する)

...以下詳細...

気象庁から送られてくる震源・震度に関する情報のXML電文には、発生時刻や地震の規模を表すマグニチュード、震源地の座標、津波に関するコメント、そして揺れを観測した地域の地域コードと地名、震度が含まれています。
このXML電文から、震源地の座標を地図にプロットするのは非常に簡単ですが、各地の震度を、その地域の中心となる場所にプロットするというのが結構面倒です。
というのも、気象庁は地域コードや地名のテーブルは公開していますが、その中心(厳密に言えば重心)座標は提供していません。
また、気象庁の区域割りはちょっと特殊で、例えば「石狩地方北部」みたいに、ある程度の市区町村をまとめていたりします。
そこで、まずこれらの区域の中心座標を求めよう、ということで...

1. PostgreSQLとPostGISで国土数値情報(行政区域)を扱ってみる (2015/06)
国土交通省から提供されている行政区域(群市区町村みたいな区域割り)のShapefileをPostgreSQLに読み込んで

2. PostGISで気象庁の細分区域&市町村等に対応する重心を求める (2015/06)
PostGISを使って気象庁が提供している区域とすり合わせて、重心を出力し、テーブルを作り

3. OpenLayers 3で気象庁発表の震度をマッピングする (2015/07)
作成した重心テーブルと気象庁XML電文からOpenLayersなどの地図システムで扱えるGeoJSON形式のデータを生成する

という流れで、震度を地図上のその地域の上に表示させることが出来ました。
ちなみに、表示される位置が中心ではなく重心なのは、地域の形が収まるバウンディングボックスから中心を求めてしまうと、突起やへこみ、極端に長細い岬が出ているような地形ではそれらに引っ張られて視覚的(感覚的?実際的?)な中心から外れてしまうからです。
重心をSQLを叩いてパッと計算(全国分計算するのには数分掛かりましたが)することができるのは、PostGIS最強としか言えません。

この生成した重心座標のテーブルは、上にもリンク張っていますが、GitHubに置いているので、これから気象庁による区域割りの地図上に何か情報を表示させるようなことをしたい方は使ってみてください。

裏側の話


昨年作った桜島の噴火情報配信サイト「桜島詳報」と若干設計が異なります。
以下は、桜島詳報とプッシュ通知システムの構成。


さらに詳しいスライド版はこちら

桜島詳報ではSubscriber自身に色々と他の処理もさせていました。
そして、詳細ページにアクセスする度にXML本体をパースしてサーバ側でページを組み立ていました。

今回のシステムでは、SubscriberはXMLの保存と、電文の種類別に行う動作を振り分けるBridgeに概要を渡すところまでさせています。
それと、概要とUUIDを一緒にMySQLへ入れていた処理を廃止して、XMLをJSONに変換してMongoDBに入れるようにしています。
(桜島詳報のシステムを合体させるために、元のXMLも一緒に保存しています)


今回のシステムでは電文が届く度にjsonが生成され、詳細ページにアクセスすると生成されたjsonを取得してクライアントサイドでページを組み立てていく仕組みです。

届いた電文が震源・震度情報の場合の処理は以下のとおり。


将来的にAPIとか公開するかもなぁーということで、複数のサーバに分けています。
(現時点では性能に余裕があるので1つのVPSの中に上記4つが詰まっている)
地震火山詳報は上記のEVIの部分になります。
今回新たにEDSという地震に関する情報を保存＆提供するサーバと、AZMAPという日本の地図タイルを持っているサーバを作りました。

地震火山詳報の震源・震度情報一覧から任意の1ページを開くと、info.json(発生日時や最大震度、マグニチュードなどが入っている)、smallScalePoints.json(広域ポイント)、largeScalePoints.json(詳細ポイント)の3つがEDSからダウンロードされて、JavaScript(jQuery, OpenLayers3)で処理されます。

地震火山詳報のページ内では提供していませんが、情報を各種SNSに流す時に付属させるデータとして、震度マップをWebKitに描画させてキャプチャした画像も生成しています。
こんな画像がTwitterに震源・震度情報の概要と共に投稿されます。


似たような震源・各地の震度画像出力エンジンを作っている会社もありますが、こちらは普通のHTML/CSS/JSで記述されたWebページから画像を生成しているので、簡単に出力する画像のデザインが可能というわけです。

この画像地図は詳細ページの地図と違って、OpenStreetMapの地図タイルを使用していません。
国土交通省の国土数値情報(行政区分)Shapefileから地図タイルを生成して、使用しています。(この地図タイルを持っているのがAZMAPサーバ)

関連する話


これも公開していないのですが、気象庁XML電文が届く度にheadlineをJSONに変換して流しているWebSocketサーバも作って使っています。
これは、グローバルIPを持つサーバが無いと受け取れない気象庁XML電文を、どこでも受信できるように...という考えで立てているのですが、キャパシティがそんなに無いので(今のところ)自分用です。
これをRaspberry Piで受けとって64x16ドットのフルカラーLEDマトリクスに表示させたりしている話は、また改めて記事にしようと思います。

感想と告知


  • 自分の想像(妄想)を形にするものづくりは楽しい。生きている、生きていける理由。
  • Web系欲が満たされたので、またしばらく組み込み系に戻ります。
  • 3か月くらい前から、はぐれエンジニアになって実家療養中です。

Time: 4:02 Tag: , , , ,
0 Comments

2015年12月5日

OL3で扱える地図タイルをShapeから生成する

このエントリーをブックマークに追加 このエントリーを含むはてなブックマーク
OpenLayers3で扱える地図タイルを国土交通省 国土数値情報(行政区域)から生成する方法のメモ。
この方法で出力すればOpenLayers3以外(Leafletなど)でも利用可能です。

使用するのは

1. Shapefileからmbtilesを生成する


LearnOSMに詳しい手順が載っているので、おおまかな手順だけ。

1. TileMillを起動し、New projectから新しいプロジェクトを作成する
2. 適当なFilenameを入力し、Image formatを選択し、Addをクリック
この時、デフォルトの世界地図を利用しない(読み込むShapefile以外必要ない)場合は、Default dataのチェックを外す
3. Editor画面になったら、左下のLayersからAdd layerをクリック
4. Datasourceで読み込みたいShapefileを選択し、Save & Styleをクリック
5. Layerが追加されるので、適宜デザインを調整する(デザインはCartoCSSで記述可能)
6. 右上のExportからMBTilesを選択する
7. 書き出す範囲やズームレベルを以下のような形で選択し、Export
おおまかな容量がZoomの下に表示されるので、参考程度に。
(100GB+とかの場合、出力に結構な時間を要する＆ファイルサイズが大きくなるので注意)

2. mbtilesから地図タイルを生成する


今回はTileMillで出力したjapan.mbtilesからjapan/tiles下に"Z/X/Y.png"の形で出力した。
git clone git://github.com/mapbox/mbutil.git
cd mbutil
python setup.py install
mb-util japan.mbtiles japan/tiles

tilesディレクト下の構造を保ったまま、適宜設置する。

3. OL3で開く


ol.min.jsなどを呼び出してあげてから、以下のような感じで。
var map = new ol.Map({
    target: 'map',
    renderer: 'canvas',
    layers: [
        new ol.layer.Tile({
          source: new ol.source.XYZ({
            url: '/tiles/{z}/{x}/{y}.png'
          })
        }),
        new ol.layer.Tile({
          source: new ol.source.TileWMS({
            attributions: [new ol.Attribution({
              html: "地図データ © 国土交通省 国土数値情報(行政区域)"
            })]
          })
        })
    ],
    view: new ol.View({
        //中心座標(仮指定)
        center: ol.proj.transform([134.15, 35.27], 'EPSG:4326', 'EPSG:3857'),
        zoom: 6
    })
});

(OpenLayers 3.11.2で動作確認済み)

WGS84あたりに測地系を揃えてあげれば、QGISなどで作ったShapefileも同じ手法で地図タイルにできるので、色々なデータをOL3で見ると面白い(かも)
Time: 21:39 Tag:
0 Comments

2015年7月5日

OpenLayers 3で気象庁発表の震度をマッピングする

このエントリーをブックマークに追加 このエントリーを含むはてなブックマーク
色々忙しくて、前回から丸1ヶ月空きました...
気象庁の防災情報XML電文で発表される震源地と各地の震度情報をOpenLayers3で表示しようプロジェクト、その3です。

その1 : PostgreSQLとPostGISで国土数値情報(行政区域)を扱ってみる
その2 : PostGISで気象庁の細分区域&市町村等に対応する重心を求める

こんな感じになりました。


ココマデ長イ道ノリダッタ...

1. XML電文からGeoJSONを作る

まず、PostGISで気象庁の細分区域&市町村等に対応する重心を求める で生成した jma_area_centroid.csv と jma_city_centroid.csv を使って、気象庁防災情報XML電文の震度速報、震源・震度に関する情報から、区域別の震度(数値)情報と重心(緯度経度)を持つGeoJSONを出力します。

コードはまるっと以下に。
9SQ/jma-eqxml2geojson

#相変わらずのウンコードですが、動けば良いの精神で。

動作は
1. jma_area_centroid.csv と jma_city_centroid.csv をメモリ上のSQLiteデータベースにそれぞれ展開
2. 電文XMLを開いてBeautiful Stone Soupでパース
3. 震源地(epicenter)はそのまま経度緯度を取り出して辞書に入れる
4. 震度毎にリストを作り、各区域の区域コードに対応する重心(経度緯度)をデータベースから探してきてリストに追加
5. 上の3と4を一緒にしてjsonで出力
という感じ。

2. OpenLayers 3でGeoJSONのデータを描画

できたGeoJSONをOpenLayers 3でOpenStreetMap上に描画します。

コードはまるっと以下に。
9SQ/seismic-intensity-map

データが揃えばOpenLayers 3のパワーでサクッと表示できてしまいます。

XML電文の震度情報は、細分区域と市町村等区域でやってくるので、縮尺に応じて表示を切り替えするようにしました。
ズームして、ある一定以上になると市町村等レベルの震度表示に切り替わります。

初期のズームレベルと中央は、細分区域レベルでのポイントから、だいたい全体が入るように自動調整されます。

これが100行(スタイル定義を除けば実質60行)くらいで実現できるので、OpenLayers 3 凄い。

デモを以下に設置しましたので、ご自由にご覧ください。
http://www.quitsq.com/demo/seismic_intensity_map/
(2015年5月30日 20:24分頃に小笠原諸島西方沖で発生した最大震度5強の地震)

以上で終わりです。

もう少し暇になったら、地震発生毎に各地の震度情報を上記デモのように表示できるサービスを作ろうと思います。
(ついでに、Yoとかtwitterへのツイート機能も含めて。)

追記 2015/12：作りました→ 地震とか火山噴火情報を閲覧できるWebサイトを作った
Time: 23:03 Tag: , , , , ,
0 Comments

2015年6月5日

PostGISで気象庁の細分区域&市町村等に対応する重心を求める

このエントリーをブックマークに追加 このエントリーを含むはてなブックマーク

震源地と各地の震度情報をOpenLayers3で表示するために、気象庁が独自に分けた区域の重心をPostgreSQLとPostGISで求めようというお話。
今回は、国土数値情報(行政区域データ)をPostgreSQLに入れてQGISで閲覧した前回の続きです。

1. 気象庁の資料からテーブル作成


1. 気象庁防災情報XMLフォーマット 技術資料から個別コード表(jmaxml_20150430_Code.zip)をダウンロード
2. 地震火山関連コード表.xlsのシート24の内容をコピー、新しくファイルを作る
3. AreaInformationCityのコードを隣の列にコピー、下2桁を削って6桁にし、5桁(先頭0埋め)にする
4. CSVで書き出し

テーブルの作成
CREATE TABLE jma_areacode
(
  areaforecastlocale_code integer,
  areaforecastlocale_name character varying(10),
  areaforecastlocale_kana character varying(50),
  areainformationcity_code integer,
  areainformationcity_code_s character varying(5),
  areainformationcity_name character varying(20),
  areainformationcity_kana character varying(50),
  pointseismicintensity_code integer NOT NULL,
  pointseismicintensity_name character varying(30),
  pointseismicintensity_kana character varying(100),
  CONSTRAINT point_pkey PRIMARY KEY (pointseismicintensity_code);

データベースへCSVを読み込む
COPY jma_areacode FROM '/Users/keiichiro/jmaxml_code_s.csv' WITH ENCODING 'sjis' HEADER CSV;

試しに市町村等区域でグループ化してみる
SELECT
  min(AreaForecastLocalE_code) AS area_code, 
  min(AreaForecastLocalE_name) AS area_name, 
  AreaInformationCity_code_s AS city_code, 
  min(AreaInformationCity_name) AS city_name 
FROM
  jma_areacode
GROUP BY
  AreaInformationCity_code_s
ORDER BY
  city_code;


2. 市町村等区域で重心を求める


行政区域(国土数値情報)を気象庁の市町村等区域で結合
CREATE TABLE jma_city_shapes AS
SELECT 
  min(areainformationcity_code) AS city_code,
  n03_007 AS city_code_s,
  min(areainformationcity_name) AS city_name,
  ST_Union(geom) AS geom
FROM
  jma_areacode,
  shapes
WHERE
  jma_areacode.areainformationcity_code_s=shapes.n03_007
GROUP BY
  n03_007
ORDER BY
  n03_007 ASC;
ALTER TABLE jma_city_shapes ADD PRIMARY KEY (city_code);

Core i7-3820QM 2.7GHz, 16GB RAMなMacBook Proで537121ミリ秒、だいたい9分掛かりました。

ここで問題点が2つ。

その1. 五島と甑島の区分
気象庁区分では、
・長崎県佐世保市と長崎県佐世保市宇久島
・鹿児島県薩摩川内市と鹿児島県薩摩川内市甑島
は別域になっているので、それぞれ分割してarea_code,nameとcity_code,nameを修正、追加

730, 長崎県北部, 4220201, 佐世保市
737, 長崎県五島, 4220202, 佐世保市宇久島
770, 鹿児島県薩摩, 4621501, 薩摩川内市
775, 鹿児島県甑島, 4621502, 薩摩川内市甑島

その2. 気象庁と国交省のデータ更新日の差
栃木県下都賀郡岩舟町が栃木市に編入合併された
→栃木市にST_Unionする

気象庁の市町村等区域で重心を求める
ALTER TABLE jma_city_shapes
ADD centroid geometry;
UPDATE jma_city_shapes
  SET centroid = t.c
  FROM (SELECT city_code, ST_Centroid(geom) AS c FROM jma_city_shapes) AS t
WHERE jma_city_shapes.city_code = t.city_code;



3. 細分区域で重心を求める


市町村等を気象庁の細分区域(188区域)で結合
CREATE TABLE jma_area_shapes AS
SELECT 
  areaforecastlocale_code AS area_code,
  max(areaforecastlocale_name) AS area_name,
  ST_Union(geom) AS geom
FROM
  jma_areacode,
  jma_city_shapes
WHERE
  jma_areacode.areainformationcity_code=jma_city_shapes.city_code
GROUP BY
  areaforecastlocale_code
ORDER BY
  areaforecastlocale_code ASC;
ALTER TABLE jma_area_shapes ADD PRIMARY KEY (area_code);
これは819513ミリ秒(だいたい14分)で終了。

気象庁の細分区域(188区域)で重心を求める
ALTER TABLE jma_area_shapes
ADD centroid geometry;
UPDATE jma_area_shapes
  SET centroid = t.c
  FROM (SELECT area_code, ST_Centroid(geom) AS c FROM jma_area_shapes) AS t
WHERE jma_area_shapes.area_code = t.area_code;



テーブルを出力
COPY (SELECT area_code, area_name, ST_Y(centroid) AS lat, ST_X(centroid) AS lon FROM jma_area_shapes ORDER BY area_code ASC) TO '/temp/jma_area_centroid.csv' DELIMITER ',';
COPY (SELECT city_code, city_name, ST_Y(centroid) AS lat, ST_X(centroid) AS lon FROM jma_city_shapes ORDER BY city_code ASC) TO '/temp/jma_city_centroid.csv' DELIMITER ',';


生成したデータはGitHubに置いてます。
https://github.com/9SQ/jma-eqarea-centroid

次回は、OpenLayers3で実際に震源地・震度の表示をします。

NEXT→ OpenLayers 3で気象庁発表の震度をマッピングする
Time: 12:36 Tag: , , ,
0 Comments

2015年6月4日

PostgreSQLとPostGISで国土数値情報(行政区域)を扱ってみる

このエントリーをブックマークに追加 このエントリーを含むはてなブックマーク

先日(5/30)、小笠原諸島西方沖を震源とする地震がありましたが、そういや気象庁防災情報XMLで震源地や各地の震度情報なども配信していたなぁと思いだし...
震源地と各地の震度情報をOpenLayers3で表示するのを作ろう！
と思い立ったのが4日前。
早速作ろうと色々データを集め始めたのですが、表示の際に一番必要な各地の座標(緯度,経度)は気象庁から公開されていません。
地震が発生したときに気象庁から送信されるXML電文は、震度速報、震源に関する情報、震源・震度に関する情報などがあり、例えば震源・震度に関する情報では...
(前略)
<Pref><Name>東京都</Name><Code>13</Code><MaxInt>5+</MaxInt>
 <Area><Name>小笠原</Name><Code>359</Code><MaxInt>5+</MaxInt>
  <City><Name>小笠原村</Name><Code>1342100</Code><MaxInt>5+</MaxInt>
   <IntensityStation><Name>小笠原村母島</Name><Code>1342103</Code><Int>5+</Int></IntensityStation>
   <IntensityStation><Name>小笠原村父島西町</Name><Code>1342100</Code><Int>4</Int></IntensityStation>
   <IntensityStation><Name>小笠原村父島三日月山</Name><Code>1342101</Code><Int>4</Int></IntensityStation>
  </City>
  </Area>
 <Area><Name>東京都23区</Name><Code>350</Code><MaxInt>4</MaxInt>
  <City><Name>東京千代田区</Name><Code>1310100</Code><MaxInt>4</MaxInt>
(後略)
このようになっています。
ここで、Areaは気象庁が独自に分けた188区域を、Cityは総務省の市区町村コードをベースにした1898区域が指定されます。

参考：気象庁 | 震度情報や緊急地震速報で用いる区域の名称

このXMLには区域コードと区域名だけで、気象庁が公開している技術資料などにも各区域のコードと名称のリストはありますが、座標は無い...

というわけで、無いなら作りましょう。

今回はPostgreSQLと、その拡張でGISオブジェクを格納することができるPostGIS、GISの閲覧、編集、分析機能を持ったQGISを使って、国交省が公開している国土数値情報(行政区域データ)を扱ってみます。

1. PostgreSQL, PostGIS, QGISなどをインストール (Mac/homebrew)


まずはQGISをインストール
brew install pyqt
brew install caskroom/cask/brew-cask
brew cask install qgis

次に、PostgreSQLとPostGIS、pgAdmin3のインストール
brew install postgresql
brew install postgis
brew cask install pgadmin3

2. データベースの作成、データの読み込み


データベースを作ってログインします
postgres -D /usr/local/var/postgres
createdb -E utf8 gis_dataset
psql -U keiichiro -d gis_dataset

エクステンションを設定してPostGISを使えるようにします
CREATE EXTENSION postgis;
CREATE EXTENSION postgis_topology;

PostGISが読み込めているか確認してみる
SELECT PostGIS_Version();

読めてたらpsqlを抜ける(\q)

国土数値情報 行政区域データ(全国)をダウンロードする
http://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03.html

N03-140401_GML.zipを展開して、SQLに変換、UTF-8に変換する
unzip N03-140401_GML.zip
cd N03-140401_GML
shp2pgsql -W SHIFT-JIS -c N03-14_140401.shp shapes > N03-14_140401.sql
nkf -W N03-14_140401.sql > N03-14_140401_UTF8.sql

先ほど作ったデータベースに読み込む
psql gis_dataset < N03-14_140401_UTF8.sql

3. QGISで見てみる


QGISを開き、レイヤ->レイヤの追加->PostGISレイヤの追加


[新規]をクリックして、ホスト(localhost)、先ほど作ったデータベース名、ユーザ名を入力して[OK]


先ほど取り込んだデータは、publicスキーマのshapesテーブルに入っているので、これを選択して[追加]


測地系を選択(WGS 84がデフォルトで選ばれているはず)して[OK]


これで、取り込んだ国土数値情報 行政区域データがレイヤとして追加され、画面に表示されるはず。


元となる行政区域データの準備が整ったので、次回は気象庁のデータを取り込んで、重心を求める計算をします。

NEXT→ PostGISで気象庁の細分区域&市町村等に対応する重心を求める

参考
http://qiita.com/nakamods/items/7f25aaaba950c8b8a458
http://morphocode.com/how-to-install-postgis-on-mac-os-x/
http://qiita.com/yellow_73/items/0845451b792f4bc33e90
Time: 16:06 Tag: , ,
0 Comments
登録: 投稿 (Atom)