地形の標高の格子データを作図する例として、 [地形 (標高) データを Pythonで読み、作図する (試作版)] の ページでは、アメリカ合衆国の NCEI (旧 NGDC) で編集された世界の標高・水深データ ETOPO1 の 日本の部分を切り出したものをつかった。
ここでは、日本の国土交通省の 「国土数値情報ダウンロード」のサイト https://nlftp.mlit.go.jp/ksj/ にある 「国土数値情報」 {1. 国土 (水・土地)」 「地形」 のうちの 「標高・傾斜度3次メッシュデータ」 https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-G04-a.html (いわゆる 1 km メッシュ) の標高データを読んで作図してみる。
これは、国土地理院の「基盤地図情報」の「10 m メッシュ」の標高データをもとに 2011年度に編集されたものである。
【日本語では 1980年代あるいはそれよりまえから、この種類のデータセットを「メッシュデータ」ということが多い。 「メッシュ」の語源は英語 mesh にちがいないのだが、 英語でこのようなものをさす表現は mesh ではなく、grid あるいは raster である。】
【国土数値情報の「1 km メッシュ」の標高データは、 それよりもまえの 1981 年度につくられたものがあるが、これは、測地系 (緯度経度の基準) が当時の標準だった「日本測地系」になっている。 「基盤地図情報」やそれにもとづいて編集された国土数値情報の新版のデータセットは「世界測地系」によっている。】
国土数値情報は、日本付近の地球表面を、まず、緯度{2/3}度、経度 1度 の箱にわける。 これを「1次メッシュ」という。これはほぼ国土地理院の 20万分の1地勢図の図幅に対応している。 つぎのリンクさきに1次メッシュ番号と図幅の対応を示す。たとえば「東京」図幅は 5339 である。 (ただし、リンクさきのページの情報は2000年現在のものであり、 その後に国土数値情報が対象とする1次メッシュの箱はさらに追加されている。)
「1次メッシュ」の箱を、東西、南北それぞれ 8等分したものが「2次メッシュ」である。 これはほぼ、かつての 2万5千分の1地形図の図幅に対応している。(現在の図幅は変更されている。)
「2次メッシュ」の箱を、東西、南北それぞれ 10等分したものが「3次メッシュ」である。 緯度幅は {1/120}度、経度幅は {1/80} 度であり、長さにして 1 km に近いので、「1 km メッシュ」ともいわれる。
データは、1次メッシュごとに、たとえば「G04-a-11_5339-jgd_GML.zip」のような名まえの zip アーカイブになっている。 これを全国ぶんダウンロードした。 (ひとつひとつポチポチやるよりはてまが少ない方法をみつけたのだが、 どのような手順だったか、わすれてしまった。)
ダウンロードしたものを展開すると、 ファイル名の拡張子が gfs, xml, dbf, shp, shx であるものがふくまれている。 このうち xml, shp, dbf に、実質データが重複してふくまれているらしい。
ダウンロードサイトにある「操作マニュアル」には、QGIS というソフトウェアをつかって shp ファイルを読む方法がかかれている。 わたしの Windows パソコンでも、読みこんだデータに対応する位置に四角が表示されるところまではいったが、 それからどうすれば標高の数値がわかるのかは、QGISを勉強しないとわからない。
xml はテキストファイルであり、そのうち数値データの部分は csv の形でふくまれているようだが、 それを正しくとりだすには xml を解釈できるソフトウェアが必要であり、わたしはまだもっていない。
dbf は、8ビットパソコンの時代からつかわれていた「dBASE II」のデータ形式にちがいない。 LibreOffice Calc や Microsoft Excel でも読める。 Python には dbfread というパッケージがあり、 Anaconda Navigator の environment からインストールすることができた。
そこで、まず、dbfファイルのリストをつくった。 (Unix互換のフリーソフトウェア UnxUtils の ls コマンドで「ls -1 *.dbf > ls_dbf.txt」とした。 ただしこれは Windows の「コマンドプロンプト」で実行した。 もし Anaconda Power Shell の中でコマンド標準出力でファイルをつくると 文字コードが UTF-16 になり、コード変換が必要になる。)
そして、1次メッシュの箱ごとに、 dbf ファイルを読み、その内容を Pandas のデータフレームにつめかえた。 それを実行したプログラムは、下にしめす dbf_to_raster.py の前半部分にあたるものである。
データフレームには 10の列がある。
データ欠損のところには「unknown」という文字がはいっている。
そのうち平均標高 (G04a_002) の値を、NumPy の 32ビット (4バイト) 浮動小数点数の 80×80の配列に つめかえて、単純バイナリ形式でファイルに書きだすプログラムをつくり、 全部の dbfファイルについて実行した。 「unknown」のところは、NumPy の nan (not a number) を入れた。
1次メッシュの複数の箱のラスターデータを連結するプログラムをつくり、 ひとまず、関東地方をふくんでそれよりやや広い範囲を連結してみた。 連結範囲は、1次メッシュ番号の上2けた (南北)、下2けた (東西) で指定している。 なお、長方形領域のうちで dbfファイルがないところは、nan でうめた。
なお、このバイナリデータにつぎの ctl ファイルを追加すれば、GrADS で読むことができる。 【「UNDEF -99999」は GrADS の ctl ファイルの形式をみたすためにバイナリデータに出てこない値を書いただけであり、 バイナリデータでは欠損は -99999 ではなく nan であらわされている。 GrADS では nan は仕様の外のようだが、これまで見たところでは欠損あつかいされているようだ。】
[地形 (標高) データを Pythonで読み、作図する (試作版)] の 例と同じ空間領域について、同じ段彩で表示してみた。