世界測地系の地図上の日本測地系3次メッシュポリゴンを作る(5)

QGISでシェープファイルの表示設定

QGIS（http://qgis.org/index.php）は今日（6/19），新バージョンの1.7.0をリリースしたそうです。でもこのエントリは1.6.0を使って書いています。
前回（http://d.hatena.ne.jp/kooi/20110618）は国土地理院の基盤地図情報をダウンロードし，シェープファイルに変換しました。今回はそれをQGISで表示，印刷します。
QGISを起動し，新しい地図プロジェクトに「ベクタレイヤの追加」で基盤地図情報のシェープファイルを追加します（複数同時選択可）。都道府県の広さやパソコンの性能にもよりますが，AdmArea（行政区域［面］），AdmBdry（行政界［線］），AdmPt（市町村名［点］），RailCL（鉄道），RdEdge（道路），WL（水涯線），Cntr（等高線），Cstline（海岸線），BldA（建物［面］），ElevPt（標高点）の10個のシェープファイルを一気に読み込むことができます。

読み込んだら，表示の重なりを考えながらレイヤの順番を修正します。左側にあるレイヤのボックス内でレイヤを選択してドラッグすれば順番を変えることができます。また，それぞれのレイヤで右クリックしてプロパティを表示し，シンボル（色・線の太さなど）を設定します。

AdmPtではシンボルを非表示（シンボルの大きさを0とする）にして，ラベルを表示します。プロパティのラベルを選択して，ラベルを表示するにチェックを入れるとラベルとして表示する属性をプルダウンで選択できるようになります。AdmPtでは「名称」を，ElevPtでは「標高」を選びます。フォントサイズや色も指定します。

csv（デリミテッドファイル）の読み込み

続いて，前回ダウンロードした位置参照情報のcsvファイルを読み込みます。メニューの「プラグイン」→「デリミティッドテキスト」→「デリミティッドテキストレイヤの追加」とすると以下のような画面になります。デリミティッドテキストファイルの横のボタンをクリックしてcsvファイルを選択します。

ファイルを選択すると，csvファイルの内容が表示されます。1行目にカラムのタイトルが入っています。その中からXに経度，Yに緯度を選択します（順番に注意）。

レイヤを追加すると点が表示されます。他のPtレイヤと同様にシンボルを非表示にしてラベルに大字町丁目名を表示させます。
なお，デリミティッドテキストレイヤも右クリックで「名前を付けて保存」からシェープファイルに変換して保存できるはずですが，今回の例に使っているcsvファイルは1行目のカラムのタイトルが日本語文字で6文字以上になっているためエラーになります。シェープファイルに変換したい場合は，csvファイルの1行目に入っているカラムのタイトルを全て英数字で10文字以内にしておきます。（シェープファイルを構成するdbfの制約）

投影座標系への表示の切り替え

緯度経度表示にしていると，地図の南北方向と東西方向で縮尺が異なってしまうので，「設定」→「プロジェクトのプロパティ」で「'オンザフライ'CRS変換を有効にする」にチェックを入れ，座標系を平面直角座標系（画面の例では佐賀県なのでII系）にします。表示が消えますが適当なレイヤを選んで右クリック→「レイヤの領域にズーム」とすれば再度表示されます。

なお，デリミティッドテキストファイルを読み込む際，XとYを経度緯度で指定しますが，「'オンザフライ'CRS変換を有効にする」がチェックしてあればプロジェクトのプロパティが平面直角座標系であってもちゃんとした場所に点が表示されます。ですから，最初にQGISを起動してレイヤの追加をする前にプロジェクトのプロパティで座標系を選択しても問題ありません。
平面直角座標系については，以下を参照して，各都道府県に適した座標系を選択してください。
http://vldb.gsi.go.jp/sokuchi/patchjgd/download/Help/jpc/jpc.htm

プリントコンポーザーの設定，印刷

QGISで地図を印刷するにはプリントコンポーザーを設定する必要があります。まず，「新コンポーザマネージャ」をクリックすると以下のような画面になります。

まず地図を配置します。コンポーザマネージャの「新規地図を追加」で適当にドラッグすると長方形のフレーム（アイテム）が配置され中に地図が表示されます。配置されたフレームを選択して，右側の「アイテム」タブの中で地図フレームの大きさや縮尺を設定します。
地図の表示範囲は，「アイテムの中のコンテンツを移動」（地球を手でつかんでいるようなアイコン）で地図をドラッグすると動かすことができます。

また，スケールやノースマーク，凡例などを表示することもできます。スケールバーは「新規スケールバーを追加」でアイテムを追加し，アイテムタブの中で設定します。

ノースマークは，「イメージを追加」（カメラのアイコン）で追加するアイテムとして設定することができます。

印刷する地図のイメージが設定できれば，プリンタに印刷することができます。また，pdfやsvgへの出力も可能です。ただし，DistillerやPDFMakerが入っているなら，「印刷」でpdfファイルを作る方がきれいにできるように思います。

これで日本測地系3次メッシュをつかった自然環境調査のための地図を印刷する話は一区切りとします。ここで印刷した地図を使えば，パソコンやGPSが得意でないおじさまやおばさまでも，サンプリング位置を記録してもらえると思います。
次は野外調査データをQGISで整理する話へいきたいと思います。GPSを使う場合と使わない場合を考えています。記事タイトルがだいぶずれてきましたので，次回からタイトルを改めます。（その前に，間奏エントリをあげるかもしれません。）

世界測地系の地図上の日本測地系3次メッシュポリゴンを作る(4)

基盤地図情報の利用

前回（http://d.hatena.ne.jp/kooi/20110612）はWMS地形図の表示まで紹介しましたが，印刷はできないところで終わりました。今回は，使いたいところの地図をメッシュ入りで印刷するために必要なデータを紹介します。
さて，基盤地図情報というのは，国土地理院が公開している日本国内の地形図や都市計画図の基盤となる地理情報です。
http://www.gsi.go.jp/kiban/index.html
基盤地図情報の閲覧・ダウンロード→ダウンロードサービスを開きます。
http://fgd.gsi.go.jp/download/
6月18日現在，ID登録制へ移行する予定というアナウンスは出ていますが，自由にダウンロードすることができます。まずは下の方にある基盤地図情報閲覧コンバートソフト（6.2MB zipファイル）をダウンロードします。zipファイルを解凍するとFGDV.exeなど8つのファイルができます。FGDV.pdfがマニュアルになっています。
インストールをしなくてもexeファイルをダブルクリックすると動きますが，Start.cmdを実行して任意のフォルダにインストールしておきます。

次に，データのダウンロードです。基盤地図情報ダウンロードサービスのダウンロードファイル形式選択から「基盤地図情報縮尺レベル25000」「JPGIS形式」をクリックします。するとダウンロード項目指定として都道府県名がずらっと並びます。
ダウンロードしたい県の「+」をクリックすると，「海岸線」「行政区画の境界線及び代表点」「道路縁」「軌道の中心線」「標高点（数値標高モデルを除く）」「水涯線」「建築物の外周線」のチェックボックスがあらわれます。全部チェックして「選択して次へ」のボタンをクリックすると，ダウンロードファイルリストが表示されます。「道路縁」「標高点（数値標高モデルを除く）」「建築物の外周線」の３つは特にファイルサイズが大きく，県によっては20〜30MBごとに複数のファイルに分割されていると思います。北海道を見たら「標高点（数値標高モデルを除く）」のzipファイルは71個もありました。これ，どの地域がどのファイルに入っているかわからないのはちょっと大変ですね…
ちなみに，「標高点（数値標高モデルを除く）」とは等高線のことです。

ダウンロードファイルリストの右端のダウンロードボタンを押して，一つずつファイルをダウンロードします。ダウンロードしたファイルは，適当なフォルダにまとめておきます。zipファイルを解凍する必要は特にありません。
基盤地図情報閲覧コンバートソフトを起動し，「新しいプロジェクト」で「追加」ボタンを押すとファイルを選択する画面になりますので，今ダウンロードしたzipファイルを選択します（複数同時選択可）。

見たいファイルが全部ボックス内のリストに表示されたら，「OK」をクリックすると読み込みがはじまります。Core2+メモリ2GBぐらいなら比較的さっと動くと思いますが，わがノートパソコン（ThinkPad T42）では佐賀県のファイルを読み込むだけで20分ぐらいかかりました。広大な都道府県では，道路や等高線は一部のzipファイルだけを読み込んでみて必要な地域が含まれているか確認したほうがよいかもしれません。

もうちょっと見たい部分を拡大すると，例えば以下のような表示になります。

次にこのデータをシェープファイルに変換します。
基盤地図情報閲覧コンバートソフトのメニューの「コンバート」→「シェープファイルへ出力」で変換項目を選択し，出力先のフォルダを指定します。今回は直角座標系へ変換はせず，佐賀県の全項目のファイルを作ります。前記ノートパソコンで10分ぐらいの処理時間でした。処理時間がかかりすぎる場合は，必要な部分を矩形や多角形で選択して，その中だけを変換することもできます。

処理が終わったら，.shp, .shx, .dbf, .prjの4種類各10個計40個のファイルができているはずです。
これで，基盤地図情報のシェープファイルができました。
なお，地域によっては，市町村ごとに1:2500（都市計画図の基図）のデータが公開されています。ただし，1:2500レベルの等高線データはありません。調査範囲や目的によっては1:2500レベルの道路や建物の地図情報をフィールド野帳用の地図に使うのもよいでしょう。

位置参照情報ダウンロードサービス

ところで，今作ったシェープファイルでは市町村名までは表示できますが，大字や町名のデータはありません。地名の表示されていない地図ではどこかわからないという方もあるかもしれませんので，もう少し細かい地名を表示することを考えてみます。
国土交通省国土計画局のGISサイトにはいろいろとGISで使えるデータが公開されています。
http://nlftp.mlit.go.jp/index.html
この中で，今は位置参照情報ダウンロードサービスを利用します。
http://nlftp.mlit.go.jp/cgi-bin/isj/dls/_choose_method.cgi
都道府県単位でダウンロード→必要な県を選択・大字町丁目レベルのみ→最新年度（平成22年度）を選択してダウンロードしてください。街区レベルというのは市街地では番地に相当するので，今回はいりません。
ダウンロードしたzipファイルを解凍するとcsvファイルと書式とメタデータのファイルができます。このcsvファイルを先ほどの基盤地図情報のシェープファイルと同じフォルダにでも保存しておいてください。
次回はQGISで基盤地図情報のシェープファイルと位置参照情報csvをひらいて，メッシュと一緒に印刷します。

GIS関係のエントリの「インテルメッツォ」

QGISの使い方指南のような記事を書き始めたところ，思いもかけずいろいろな方に読んでいただいているようです。ありがとうございます。ちょうど個人的にも，職場にユーザを増やしたくてQGISの使い方についてまとめようとしていたところなので，今後も続けて書いていこうと思っています。
GISというとどうしても高価であつかいが難しいソフトウェアというイメージが強いですし，実際，それなりの費用をかけて導入はしたものの活用されないままバージョンが古くなって今さら使えなくなったという話も聞きます。日常的に，スキャンした地形図の上にイラストレータを使って確認位置を落としていって，報告書用の図面を作成している私の職場でもGISの活用は余り進んでいません。
GISをはじめて使ってみようとする人にとって，わかりにくいのは，まず座標系の緯度経度と投影座標系のことのようです。GISの緯度経度表示で出力した図面をイラストレータでむりやり縦横比を変えてスキャンした紙地図に重ねていたのを見たこともあります。今のQGISでは，同じ測地系のデータであれば地理座標系も投影座標系も自動的に変換して重ねてくれるので，あまり問題になりません。今回の一連のエントリでは，基本的にデータは緯度経度で，表示は平面直角座標系を使う予定です。
もう一つわかりにくいのが測地系で，平成14年の測量法改正前まで使われていた日本測地系と現在使用されている世界測地系の話が，初心者にとって最大の難関のように思います。何が難関かというと，「ベクタ→データマネージメントツール→新しい座標系へエクスポートする」では日本測地系と世界測地系の座標変換ができない点でしょう。いきなりこの一連のエントリの最初でExcelマクロやらRが出てきたのは，これが大変だからです。
私の記事の想定読者は，マクロ生物学（分類・生態）の学生さんとか環境調査関係者です。GISを使う目的も，まずはフィールドへ持参する地図をきれいに印刷すること，次にフィールドでとってきたデータを地図上で整理すること，そして地形や植生など関連情報と重ね合わせて分析することと進めていく予定です。

なお，この一連のエントリのきっかけとなったtwitterのやりとりをtogetterに保存しています。微妙にからむ別の話題も入っていますが，非公開にしていたのを公開しました。
http://togetter.com/li/147487

次回は基盤地図情報をダウンロードしてシェープファイルに変換して使用する方法を紹介しようと思っています。

世界測地系の地図上の日本測地系3次メッシュポリゴンを作る(3の続き)

WMSサービスの地図が平面直角座標でうまく表示されない件

平日も毎日更新はつらいので，続きのエントリはしばらくお待ちいただきたいのですが，昨日できなかった，基盤地図WMSを平面直角座標で表示することについて教えていただいたりしたので補足します。
http://www.finds.jp/siteinfo/news/20110607-01.html.ja
http://d.hatena.ne.jp/yellow_73/20110607

難しいことはよくわかりませんが，サーバから送る書式をWMS1.3.0という新しいバージョンに変更したため，座標の数値の順番が入れ替わってしまったということのようです。
QGISでは，サーバのurlに"VERSION=1.1.1&"を追加すれば，元通り平面直角座標でWMSを表示させることができるようになります。（基盤地図25000，地名WMSとも）

世界測地系の地図上の日本測地系3次メッシュポリゴンを作る(3)

QGIS

オープンソースのデスクトップGISソフトウェアであるQGISは，以前はpython errorとかで止まってしまって泣くことがあったのですが，現行の1.6.0はとても使いやすいと思います。
フィールド調査に行きその調査結果を報告書や論文にまとめるような仕事をしている方はぜひ自分のパソコンに入れて使ってみてください。何十万円もする商用GISを使わなくても，かなりいろんなことができます。
ダウンロードサイト（英語）
http://www.qgis.org/wiki/Download
オープンソース地理空間ソフトウェアのサポートサイト（OSGeo.JP）
http://www.osgeo.jp/
OSGeo.JP内のQGIS ユーザーズガイド（日本語訳）
http://www.osgeo.jp/user_guide/user_guide.html

QGISをインストールして起動したら，さっそく前回作成したf3meshN.shpとf3meshW.shpを表示してみます。ツールバーの「ベクタレイヤの追加」をクリックしてC:\GIS\に入っているシェープファイルを選択します。なお，このメッシュファイルは英数字だけですが，Windowsの場合，文字コードはSHIFT-JISに変更しておいた方がよいと思います。

読み込んだら左側のレイヤのところで1つを選んで右クリックし，プロパティを表示させます。塗りや線の有無や色を変更します。

日本測地系のメッシュf3meshNを塗りなし・黒線で，世界測地系のメッシュf3meshWを青線にしました。また，メッシュ１つ１つについている属性をラベルに表示することもできます。拡大してみると，日本測地系のメッシュの方が西北側へずれていることがわかります。

ところで，座標系については後で説明しますが，この地図がどの座標系によって表示されているか確認しておきます。メニューの設定からプロジェクトのプロパティを表示させます。デフォルトの座標系はWGS84だと思います。今は特に変更する必要はありません。「'オンザフライ’座標系変換を有効にする」にチェックを入れておきます。

それから，Rで生成したシェープファイルは座標系の定義がされていないので，定義をしておきます。メニューの「ベクタ→データマネジメントツール→現在の座標系を定義する」で座標系はJGD2000にします。

基盤地図情報25000 WMS配信サービスを利用する

メッシュは表示できたけれど，地図が表示されていなければどこがどこだかわかりません。地図のシェープファイルがあればレイヤに追加すればいいのですが，それは後述するとして，まずはインターネットを使った便利な方法を紹介します。
基盤地図25000　WMS配信
http://www.finds.jp/wsdocs/kibanwms/index.html.ja
「ご使用方法」にQGISでの使い方があります。
接続したらAdmAreaBdr, Cntr10, BldA, WL, RdEdg, RailCLを選択して追加します。

表示できましたか。
え，地名表示がほしい？　では地名WMS配信サービスにも接続しましょう。
http://www.finds.jp/wsdocs/pnwms/index.html

こんな感じで表示できますか。以下は福岡県大牟田市の南部。日本測地系メッシュでは493034の2次メッシュにわずかに福岡県域が含まれています（49304490の黒枠の上の方）。実は世界測地系のメッシュ（青）では493044のメッシュにおさまります。

ところで，なんかおかしい。上下につぶれていると感じませんか。GISの緯度経度表示では，緯度の1秒（約31m）と経度の1秒（約26m，北へ行くほど短くなる）が同じ長さで表示されてしまうため，東西へつぶれた形になります。
普通の紙の地図と同じように東西も南北も距離に会わせた座標にするためには，緯度経度ではなくある決められた原点からのX軸とY軸の距離座標できめられた投影座標系に変更する必要があります。ここではとりあえず福岡県あたりの25000地形図と同じUTM52系にします。*1
再び，設定→プロジェクトのプロパティで，参照座標系を今度はJGD2000/UTM 52Nに変更します。画面から地図が消えますが，メッシュのレイヤを選択して「右クリック→レイヤの領域にズーム」とすればシェープファイルのレイヤは再表示できます。右下に表示されていた座標の値が「度」の表示からかなり大きな数字にかわったことがわかると思います。
残念ながら，WMSレイヤは座標系を変更すると表示できなくなります。いったんレイヤを削除して，再度WMSレイヤを追加してください。その際，投影法がUTM52Nになっていることを確認してください。
以下のような表示ができましたか？　今回は日本測地系のメッシュだけを赤で表示しました。

九州大学伊都キャンパスはまだ造成中ですが，道路などは平成15年ぐらいの状況が反映されているように思います。これで，20年も昔の1:50000地形図に印刷された環境庁の3次メッシュマップを使わなくても，新しい地図の上にメッシュが表示できます。
ところが，この地図を印刷しようとすると，実は問題があるのです。印刷のためにはコンポーザマネージャで新しいコンポーザを作って，地図やスケール，凡例などを設定してから印刷するのですが，WMSレイヤは基本的に印刷できないのです。すごく狭い範囲なら印刷もできるようですが，とてもフィールド調査の現地野帳用に印刷するのは無理でした。（スクリーンキャプチャーを印刷する手はありますが）
ということで今回はここまで。
次回は，インターネットで公開されている地理情報をダウンロードして利用する方法について書きたいと思います。

（続く）

*1:平面直角座標系にしようと思ったが，WMS配信地図が位置ずれを起こしてしまってよくわからなかったため

世界測地系の地図上の日本測地系3次メッシュポリゴンを作る(2)

日本測地系から世界測地系への座標変換

国土地理院のサイトから座標変換プログラムTKY2JGDと座標変換用パラメータファイルをダウンロードして解凍します。最近のパソコンなら処理時間はずっと短縮されているのでパラメータファイルは全国版で問題ありません。
http://vldb.gsi.go.jp/sokuchi/tky2jgd/download/agreement.html
TKY2JGDに一括変換させるために，m3_dsのワークシートのD，E列をテキストエディタに貼り付けて保存します。この際，貼り付けたままだとtab区切りになっているのでtabをspaceへ置換します。また，1行目はタイトルなのではじめに#をつけてコメントアウトします。
ファイル名はmesh3.inとかがいいでしょう。

#Lat Lon
341000 1300000
341030 1300000
341030 1300045
341000 1300045
...

TKY2JGDを起動し，一括変換のボタンを押します。緯度・経度→緯度・経度を選択し，入力ファイルに今保存したmesh3.inを指定します。出力ファイルはmesh3.outでいいでしょう。
変換を実行すると「1万行に1時間ぐらいかかります」と表示されますが，私のノートパソコンでも4万行に10分ぐらいでしたから，そんなにはかからないと思います。
出力ファイルは最初に12行ほどコメントがついて，変換されたデータが続きます。

#Lat Lon
 341011.47758   1295951.64632   3parameters
 341041.47411   1295951.64550   3parameters
 341041.47496   1300036.64139   3parameters
 341011.47843   1300036.64221   3parameters
...

この3parametersというのは海上など地域ごとの変換パラメータがなかった点です。
#Lat Lonの行から最後までをコピーし，エクセルファイルに貼り付けます。うまくセル区切りをしてくれなかったら，データ→区切り位置でスペースを区切り文字にして実行してみてください。
先ほどのExcelファイルのm3_dsワークシートのF列に貼り付けた場合は以下のようになります。（1行目をつけなおしています）

ID	X	Y	Lat	Lon		#Lat	#Lon	Comment
1	130.000	34.167	341000	1300000		341011.478	1295951.646	3parameters
…

シェープファイルへの変換

この緯度経度の値は，(d)ddmmss.sssssという書式ですから，これをGISに持って行くには度（ddd.dddddddd）表示にする必要があります。
そのため，J2のセルに以下の式を入れます。

=INT(G2/10000)+(INT(G2/100)-INT(G2/10000)*100)/60+(G2-INT(G2/100)*100)/3600

J2セルをJ列とK列のセルにコピーします。またこの2列のセルの書式を「数値」にして小数点以下の桁数を8にします。すると以下のようになります。GISとはXYの順が逆になることに注意してください。

最後に，このm3_dsワークシートのA列とK列とJ列をコピーし（順番注意）新しいファイルのABC列に数値と書式をペーストしてテキスト形式で保存します。ファイル名はm3N_ds.txtでいいでしょう。

ID	Xw	Yw
1	129.99767953 	34.16985488 
1	129.99767931 	34.17818725 
1	130.01017816 	34.17818749 
1	130.01017839 	34.16985512 
...

次は，統計解析ツールである「Ｒ」を使ってこのファイルをGISで使うシェープファイルに変換します。Rをインストールし，さらにshapefilesというパッケージをダウンロードしてインストールしておいてください。
作業フォルダはC:\GIS\とします。先ほど作ったm3_dd.txt, m3W_ds.txt, m3N_ds.txtをこのフォルダに入れておいてください。
Rを立ち上げ，「新しいスクリプト」を開き，以下のコードを入力します。

library(shapefiles)
dd <- read.delim("C:/GIS/m3_dd.txt")
ds <- read.delim("C:/GIS/m3W_ds.txt")
ddShapefile <- convert.to.shapefile(ds, dd, "ID", 5)
write.shapefile(ddShapefile, "C:/GIS/f3meshW", arcgis=T)
ds <- read.delim("C:/GIS/m3N_ds.txt")
ddShapefile <- convert.to.shapefile(ds, dd, "ID", 5)
write.shapefile(ddShapefile, "C:/GIS/f3meshN", arcgis=T)

コードを選択してCtrl+Rを押すと実行します。数分でshapefileができあがると思います。

QGISで表示させるところは次回に
（続く）

世界測地系の地図上の日本測地系3次メッシュポリゴンを作る

作業の流れ

1. 緯度経度による定義通りのメッシュポリゴンをつくる（Excel）
2. 日本測地系の緯度経度を世界測地系に変換する（TKY2JGD http://vldb.gsi.go.jp/sokuchi/program.html ）
3. 世界測地系の座標を十進度数表示に戻し，Shapefileに変換（Excel，R http://www.r-project.org/ ）
4. ShapefileをフリーのGISソフトで表示（Quantum GIS http://www.qgis.org/ ，基盤地図情報閲覧コンバートソフト http://fgd.gsi.go.jp/download/ ）
5. 生物記録の地理データベース構築

Excel以外はフリーウェアです。それぞれダウンロードしてインストールが必要です。
かなりちからわざですので，もっとスマートなやり方があれば教えてください。

メッシュポリゴンデータの作成（Excel2003使用）

ファイルの作成

ワークシート3枚のExcelファイルを作ります。
マクロで使う関係で，ワークシート名を次のように変えます。
Sheet1 → m2_sel
Sheet2 → m3_ds
Sheet3 → m3_dd

3次メッシュポリゴンを作る範囲の2次メッシュを選択

m2_selのシートにA1のセルから下へ順に3次メッシュを作りたい場所の2次メッシュコード（6桁数字）を入れます。Excel 2003の場合，入力する2次メッシュの数は131以下にしてください。2次メッシュ1つで500行のデータができるので，65536行を越えてしまいます。
福岡県全域をカバーしたい場合はこんな感じ。環境庁のメッシュマップで2次メッシュ87個あります。

513020
513005
513006
503075
503076
503077
503170
503060
503063
…
493042
493043
493044
493033
493034

ExcelVBAでマクロコードを入れる

Excelのツール→マクロ→Visual Basic Editorで新しいモジュールを追加して以下のコードを入れます。

Sub mesh3_generate()

Dim m2_sel As Worksheet
Dim m3_ds As Worksheet
Dim m3_dd As Worksheet

Set m2_sel = ActiveWorkbook.Worksheets("m2_sel")
Set m3_ds = ActiveWorkbook.Worksheets("m3_ds")
Set m3_dd = ActiveWorkbook.Worksheets("m3_dd")

cur_row = 2
cur_obj = 1
cur_m2 = 1

m3_ds.Cells(1, 1) = "ID"
m3_ds.Cells(1, 2) = "X"
m3_ds.Cells(1, 3) = "Y"
m3_ds.Cells(1, 4) = "Lat"
m3_ds.Cells(1, 5) = "Lon"

m3_dd.Cells(1, 1) = "ID"
m3_dd.Cells(1, 2) = "Mesh3Code"

Do While m2_sel.Cells(cur_m2, 1) > 0

    m2_c = Format(m2_sel.Cells(cur_m2, 1), "000000")
    m1_lat_c = Val(Left(m2_c, 2))
    m1_lon_c = Val(Mid(m2_c, 3, 2))
    m2_lat_c = Val(Mid(m2_c, 5, 1))
    m2_lon_c = Val(Right(m2_c, 1))
    m2_lat = m1_lat_c * 2 / 3 + m2_lat_c / 12
    m2_lon = 100 + m1_lon_c + m2_lon_c / 8
    
    For i = 0 To 9
    For j = 0 To 9
    
        m3_lat0 = m2_lat + i / 120
        m3_lon0 = m2_lon + j / 80
        m3_lat1 = m3_lat0 + 1 / 120
        m3_lon1 = m3_lon0 + 1 / 80
        m3_c = m2_c & Format(i, "0") & Format(j, "0")
        m3_lat0_dms = Int(m3_lat0) * 10000 + (Int(m3_lat0 * 60) - Int(m3_lat0) * 60) * 100 _
        + m3_lat0 * 3600 - Int(m3_lat0 * 60) * 60
        m3_lon0_dms = Int(m3_lon0) * 10000 + (Int(m3_lon0 * 60) - Int(m3_lon0) * 60) * 100 _
        + m3_lon0 * 3600 - Int(m3_lon0 * 60) * 60
        m3_lat1_dms = Int(m3_lat1) * 10000 + (Int(m3_lat1 * 60) - Int(m3_lat1) * 60) * 100 _
        + m3_lat1 * 3600 - Int(m3_lat1 * 60) * 60
        m3_lon1_dms = Int(m3_lon1) * 10000 + (Int(m3_lon1 * 60) - Int(m3_lon1) * 60) * 100 _
        + m3_lon1 * 3600 - Int(m3_lon1 * 60) * 60
        
        'メッシュの4隅の座標値。シェープファイル用に緯度経度の順序入れ替わりに注意
        m3_ds.Range(m3_ds.Cells(cur_row, 1), m3_ds.Cells(cur_row + 4, 1)) = cur_obj
        m3_ds.Cells(cur_row, 2) = m3_lon0
        m3_ds.Cells(cur_row, 3) = m3_lat0
        m3_ds.Cells(cur_row + 1, 2) = m3_lon0
        m3_ds.Cells(cur_row + 1, 3) = m3_lat1
        m3_ds.Cells(cur_row + 2, 2) = m3_lon1
        m3_ds.Cells(cur_row + 2, 3) = m3_lat1
        m3_ds.Cells(cur_row + 3, 2) = m3_lon1
        m3_ds.Cells(cur_row + 3, 3) = m3_lat0
        m3_ds.Cells(cur_row + 4, 2) = m3_lon0
        m3_ds.Cells(cur_row + 4, 3) = m3_lat0
        
        'TKY2JGD用
        m3_ds.Cells(cur_row, 4) = m3_lat0_dms
        m3_ds.Cells(cur_row, 5) = m3_lon0_dms
        m3_ds.Cells(cur_row + 1, 4) = m3_lat1_dms
        m3_ds.Cells(cur_row + 1, 5) = m3_lon0_dms
        m3_ds.Cells(cur_row + 2, 4) = m3_lat1_dms
        m3_ds.Cells(cur_row + 2, 5) = m3_lon1_dms
        m3_ds.Cells(cur_row + 3, 4) = m3_lat0_dms
        m3_ds.Cells(cur_row + 3, 5) = m3_lon1_dms
        m3_ds.Cells(cur_row + 4, 4) = m3_lat0_dms
        m3_ds.Cells(cur_row + 4, 5) = m3_lon0_dms
        
        'メッシュコード
        m3_dd.Cells(cur_obj + 1, 1) = cur_obj
        m3_dd.Cells(cur_obj + 1, 2) = m3_c
        
        cur_row = cur_row + 5
        cur_obj = cur_obj + 1
    
    Next j
    Next i
    
    cur_m2 = cur_m2 + 1

Loop

m3_ds.Range("B:C").NumberFormatLocal = "0.00000000_ "

End Sub

さて，このマクロを実行すると，m3_dsとm3_ddのワークシートにずらっと数値が並びます。
m3_ds

ID	X	Y	Lat	Lon
1	130.00000000	34.16666667	341000	1300000
…

m3_dd

ID	mesh3code
1	51302000
…

とりあえずエクセルファイルをm3_generate.xlsとか名前を付けて保存しておきます。
この時点のデータを後述のシェープファイルへの変換を行うと世界測地系のメッシュができます。m3_dsのワークシートのA〜Cの3列をコピーしてテキストエディタに貼り付け，m3W_ds.txtとして，m3_ddのデータの入っている2列を同様にしてm3_dd.txtとして保存しておきます。