エクセルで作る地形解析図【参考】
|
【参考 1 】 DEM間隔と座標系
【参考 2 】 tiff画像 【参考 3 】 マクロを知らない人へ 【参考 4 】 カラーホイール 【参考 5 】 3D表示ファイル ： VRML 【参考 6 】 局所積和演算とSUMPRODUCT関数
【参考 7 】 3D表示ファイル ： 3D-PDF
【参考 8 】 画像の合成 【参考 9 】 10mDEMと5mDEMの合成法 【参考 10 】 CSE（配列）数式 【参考 11 】 ステレオ投影 【参考 12 】 風向ローズダイアグラムとsumproduct関数 【 参考 1 】 -- DEM間隔と座標系 -- 10mDEMデータは公称10m間隔グリッドとなっているが、実際は緯度および経度を等分割して0.4秒間隔グリッドであるため緯度によってそのグリッド間隔距離は異なる。このため、対象とする地区のグリッド距離は、とりあえず正確なものを使ったほうが良い。10mDEMデータの1ピクセル（グリッド）の正確な大きさは、国土地理院のHPの「便利なプログラム・データ」で公開されている測量計算で、平面直角座標による距離計算ページを利用して求められる。 
なお、右画面画像は2013年7月くらいまでのもので、1点ごとにしか計算できなかった。それがその頃にリニューアルされて、2013年9月現在、複数点を一括して計算できるようになっている。ちなみに、本稿作成時に住んでいた北緯45度、東経142度付近では、東西方向：8.76m、南北方向：12.34mである。このため、例えば地上開度図で検索距離を250mとすると、東西方向は28ピクセル、南北方向は20ピクセル、45度方向は16ピクセルとなる。グリッドが正方形ではないので45度方向は実は正確に45度になっていないけど。 ついでに座標系のことも触れておく。 平面直角座標系というのは公共測量に使われている座標系だが、ようするにX-Y座標で位置を示している。なんだ、座標値が2点あれば距離も方向もピタゴラスと三角関数で簡単に計算できるじゃないか、と思うのだが、まぁ、メートル単位だとだいたい合ってるのだが、厳密には違う。その理由は、地球は丸いから。子午線の東西幅も北に行くほど狭くなる。つまり子午線は平行ではない。見通し範囲であれば違わないのだろうが、その曲面が微妙に影響してくる。東京を基準にして北海道で測量するとゆがみが大きくて座標軸が直角にならなくて不便だ、ということなんだろう。なので平面直角座標系は、日本国内をそれぞれの中では大きくはゆがまないだろうという範囲で19地区に分割している。 東京と北海道を含むような範囲で同じ座標で表したいという場合もあるかと思う。そういう時はUTM座標系を使う。UTMの説明は、ネット検索するといっぱい出てくる。とりあえず、日本地図センターのこちらのページをどうぞ。 UTM座標系だったら北海道と九州もいっしょに表せるかといったら、これまた難しい。UTM座標系もゾーンに区切られていて、北海道のほとんどは第54帯、九州のほとんどは第52帯で、基準原点が異なる。このスケールになると、緯度・経度で表現するしかないかな。 【 参考 2 】 -- tiff画像 -- XMLファイルをtiff画像に変換する方法について書いていたが、GRASS GISを使えばその面倒くさいことをしなくて良いので、差し替える。 ＜差し替え元＞ 【 参考 2 】 -- geotiff画像 -- 原稿差し替え（2014.9.2） ３章に、「GRASS GISを入手していれば、3）でXML文書にするより、書き換えたエクセルデータをそのままテキストエディタにコピペし、位置座標等の情報を加えてGRASS GISに取り込み、GeoTiff画像にしたほうが簡単だし、そのままGISソフトで表示・加工できるようになる。」と書いた。この位置座標等の情報とは、以下のヘッダ情報。 north: 45.16666667 south: 45.04832667 east: 141.6521111 west: 141.5818889 rows: 1065 cols: 632 null: -9999.00 type: float とりあえず適当な数値を入れているが、これはデータの範囲と行数、列数、データなしの値、データのタイプ。 ここで注意するのは、データの範囲はセルの外側であること。ワールドファイルの場合はセルの中心の座標であるので、セルサイズの0.5分ずれる。 これに表形式のデータをそのままコピペしてテキスト形式で保存したものが「GRASS ascii形式」となる。あとはこれをGRASS GISで読み込み、geotiffファイルとして保存すればよい。 参考文献 【 参考 3 】 -- マクロを知らない人へ -- エクセルのマクロなんてさわったことがない、という人も少なくないようだ。 全部を説明するのはできないので、ここではマクロのエディタを開く方法と、それを動作させる方法だけを書く。 まず、普通にエクセルシートを開いたら、『名前をつけて保存』のところでファイルの種類をプルダウンから『マクロ有効ブック』を選択し、ファイル名を入力して保存する。   次に、マクロ画面を開く。『alt＋F8』キーで右の画面が出てくる。  『マクロ名』に、作成するマクロの名前を入れると、『作成(C)』の文字が濃くなるので、クリックする。 下記の画面が出てくるので、 「Sub ＊＊()」 から 「End Sub」 の間に必要な数式等を入れる。 この「Sub ＊＊()」 から 「End Sub」がひとつのマクロとなる。 ＊＊がマクロ名になっている。   マクロを実行するには、『alt＋F8』キーで右の画面が出てくるので、『実行(R)』をクリックする。 エラーが出て途中で止まってしまったら、・・・ 申し訳ないがケースバイケースなので、勉強してくださいとしか言えない。4章で示したものは、基本的にはそのままコピペすれば大丈夫と思うが、「Sub ＊＊()」 から 「End Sub」にはさまれた範囲がひとつのマクロということだけはお忘れなく。 なお、マクロが有効なエクセルファイルが開いていれば、別なファイルからもそのマクロが使える。4.1に示したように、GML形式のXMLファイルをcsvファイルとして開き、そのcsvファイル上でマクロを実行できる。 【 参考 4 】 -- カラーホイール --
斜面方位図で8方位に区分する場合の方位の着色は、カラーホイールの並びにするのが一般的らしい。 カラーホイールは色の相対関係を示す基本となるものだそうで、カラー間の関係を把握し覚えるのも便利なもの。レッド、グリーン、ブルーは加算混合色で、シアン、マゼンタ、イエローは減算混合色。それぞれの正反対側にあるのがその補色（レッドとシアン、グリーンとマゼンタ、ブルーとイエロー）。 減算混合色の多様な色は、その補色以外の2つの加算混合色で生成される。したがって、画像内のある原色の量を増やすとその補色の量が減る。例えば、イエローはグリーンとレッドで生成され、ブルーを含まない。ということは、画像ソフトなどでイエローを調整するには、ブルー成分を増減させると、イエローが弱くなったり強くなったりするということになる。 色の話のついでだが、RGB立方体なるものも色の組み合わせの理解に役立つと思う。下左図のように三次元の直交座標軸にRとGとBを割り当てると、立方体の頂点が、R G B と C M Y 、それに黒と白の8つの色になる。 カラーホイールに対応するのが下右図のHSIカラーモデルなんだそうだけど、私にはこの図ではなんだかピンとこない。 
【 参考 5 】 -- 3D表示ファイル ： VRML --
（2014.4.27追記）
いまや、webの世界でVRMLは完全に廃れて、WebGLに切り替わっているらしい。国土地理院のサイトでは、地形図だけでなく空中写真までも簡単に3D表示して閲覧できるようになった。ただ、3Dプリンタの中にはまだVRMLファイルじゃないと模型作成できないものもあるらしい。これも時間の問題なんだろうな。 地形解析図とは直接関係ない。ただ、VRMLファイルで地形の3D表示が簡単にできることを知った。プラグインを別途インストールする必要があるが、ウェブブラウザでそのままくるくる回転させて上から、横から、下から、そして裏からものぞき見ることができる。
（表示例の右下図は、同じ箇所の段彩陰影図） VRMLファイルは、標高のグリッドデータをカンマ区切りあるいはタブ区切りなどで並べるだけのテキストファイル。サンプルとして4列5行のデータ例を右に示す。xDimensionで列数、zDimensionで行数、xSpacingで東西方向のグリッド間隔、zSpacingで南北方向のグリッド間隔を指定するだけ。右のサンプルでは標高データをわかりやすいように4列5行に配列したが、半角スペースや改行は意味がないので1データごとに改行してもかまわない。カンマの代わりにタブでもよい。 ということは、GML形式のxmlファイルがそのまま使えるということ。拡張子を .csv に書き換えてExcelで開き、B列の標高データを選択してコピーし、VRMLファイルの height [ の次の行に挿入し、拡張子を .wrl にして保存すればでき上がり。 これに衛星画像を貼り付けたりもできる。 （私の環境では、IEでは表示されるが、Firefoxだと真っ白画面となり表示されない）。チェックボックスによりパーツごとの表示非表示切り替えができる。そして下記は結構重いファイルなのだが、読み込みにはやや時間がかかるものの読み込んだ後は動きが軽いのにはびっくり。これはとあるところの地層と断層を表示した3D-PDFの例。 （マウスとCtrlキーで回転・拡大・縮小・移動が自在） |
|
【 参考 8 】 -- 画像の合成 -- 画像のソフトとしては、もっぱらphotoshopを利用しているのだが、知らない機能がいっぱいある。その中で、このサイトでも大きく関係するのが画像の合成。 2枚の画像を重ねて表示させるには、上のレイヤーを透過させる。これまでこれしかしていなかった。ところが、重ね方にもいろいろ方法があった。 レイヤーウィンドウは常時表示させていて、重ねる順番とか、不透明度とかで常に使っている。しかしその隣のプルダウンは使ったことがなかった。この中に 『乗算』 というのがある。 反転地下開度図に、地上開度図（50%透過）、傾斜量図（50%透過）をそのまま 『通常』 で表示したのが左図。傾斜量図だけ 『乗算』 で表示したのが右図。  色の調整をしていないのでちょっと暗いが、わかるだろうか。今まで赤色立体地図に似せた図をいろいろ試行錯誤していたが、どうしても谷底を暗くできなかった。ところが、乗算すると谷底がしっかり暗くなっている。このクリックに気づかず、ずいぶん無駄な労力をしていたことに、少々がっくりしてしまった。 【 参考 9 】 -- 10mDEMと5mDEMの合成法 -- 
国土地理院の基盤地図情報数値標高モデル（5mメッシュ） （以下、5mDEM） は、平野部については航空レーザースキャナ測量 （一部写真測量） の成果によるデータがかなり網羅されているが、山地部のデータは未整備なところが多い。一方、基盤地図情報数値標高モデル（10mメッシュ） （以下、10mDEM） は1/25,000地形図から作成されており、日本全域を網羅している。ただし地形図からの読み取りであるため、平坦部の小起伏はほとんど表現されない。 この例では、5mDEMがあるところは平野部の道路や河川堤防が表現されているが、ないところは扇状地みたいになり、谷間もまっ平らになっている。 ●ArcGISによる作成手順 1）国土地理院のホームページより、10mDEMおよび5mDEMをダウンロード 2）ArcGISの変換ツール（基盤地図情報）により、それぞれをインポート この時、『同一種別のデータは1レイヤーとして保存』 にチェックを入れる 3）それぞれを地理座標系から投影座標系 （平面直角座標など） に投影変換 （このとき、『セルサイズ10m、Bilinear法による補間』 とすると、10mメッシュとなる） 4）モザイク処理により、10mDEMに5mDEMを上書き 
===== 5mDEMのうち、写真測量によるデータの整備がかなり進んできた。 5mDEMでは、人工地形がきれいに反映されている。道路の切り盛り地形や田圃のあぜ道は当然くっきり。と思いながら見ていたら、右図のように山中に 奇妙なリニアメント があった。何だろうと思いつつ 空中写真を見ると、これも土地の境界だった。（2022.03.24記） 【 参考 10 】 -- CSE（配列）数式 -- 最近、CSE数式なるものを覚えた。正式には『配列数式』というらしいが、数式の入力確定に Ctrl + Shift + Enterキー を使うことから、その頭文字をとった呼び方だそうだ。CSEで入力すると、数式は自動的に ｛ ｝ で囲まれる。 例えば、sumproduct の項に書いたが、3×3の係数を次のようにして、それぞれ「keitate」,「keiyoko」と名前を付けたとする。
{＝SUM( A7:C7*{1,2,1} - A9:C9*{1,2,1} )} {＝SUM( A7:A9*{1;2;1} - C7:C9*{1;2;1} )} この場合は SUMPRODUCT関数の方が良さそうだが、例えば二次方程式 Y ＝ aX＾2 + bX + c の計算式を示す。係数 a，b，c がA1〜C1セルにあり、 X が A2セルにある場合、そのまま書くと、 ＝ A1*A2^2 + B1*A2 + C1 であるが、CSE数式では、 {＝SUM( A1:C1 * A2^{2,1,0} )} となる。 この理屈（仕組み）を知るだけでもいろいろ応用できそうで、今まで長々と書いてきた計算式が短くて済む場合があるように思える。 【 参考 11 】 -- ステレオ投影 -- 地形解析で使うかは知らないが、地質解析ではステレオ投影を使うことがある。このエクセルファイルは以前より極の投影とローズダイアグラムの表示に使っていたが、最近、気が向いて大円と小円の表示もできるようにした。ついでにカウントダイアグラムも作ってみようかとやりかけているが、まだできていない。これはそのうち。ステレオ投影ファイル 【 参考 12 】 -- 風向ローズダイアグラムとsumproduct関数 -- sumproduct関数では、countifやsumif関数と同じ計算ができる。なおかつ、countifsやsumifs関数の条件は、and関数だけだと思うのだが、sumproduct関数だとor関数も条件に入れられる。
ステレオ投影で方位別のローズダイアグラム（円形のヒストグラム）はだいぶ以前より使っていたが、気象データを見ていたら、風向ローズダイアグラムも作りたくなった。ステレオ投影では方位は10°刻みにしたが、風向は16方位なので、16分割（22.5°刻み）に変えるのはすぐできる。方位別の集計は、countif関数でできる。 そうすると次に、季節による風向の違いを知りたくなった。AMEDASの日データを過去数年分ダウンロードし、countifs関数を使うと、例えば1月の北風の日数はすぐに集計できる。データの例を右に示す。月データはダウンロードした後、日付からmonth関数で取り出す。 夏季の6月から9月の北風の日数であれば、条件を3つ並べれば良い。 =COUNTIFS([風向データ],"北",[月データ],">="&6,[月データ],"<="&9) 次に冬季の11月から3月の北風の日数を求めようとすると、11以上または3以下となり、countifs関数だと、それぞれを合計しなければならない。 =COUNTIFS([風向データ],"北",[月データ],"<="&3)+COUNTIFS([風向データ],"北",[月データ],">="&11) これがsumproduct関数を使うと、これですむ。 =sumproduct(([風向データ]="北")*(([月データ]<=3)+([月データ]>=11))) かっこがわかりにくいので単純化すると SUMPRODUCT( a * ( b + c ) ) で、 countifs関数の場合は COUNTIFS( a , b ) + COUNTIFS( a , c ) である。 ついでに書くと、6月から9月のsumproduct関数の場合はこうなる。 =sumproduct(([風向データ]="北")*([月データ]>=6)*([月データ]<=9)) これも単純化すると、SUMPRODUCT( a * b * c ) である。 これで、任意の月別で16方位ごとに日数を集計すると、風向ローズダイアグラムが描ける。よく遊びに行く「遠別」の日最多風向と、最大瞬間風速時の風向を描かせたのがこちら。   ちなみに、Excelに備わっているレーダーチャートで同じデータを描かせると、右のように似たような表現にはなる。しかし個人の感想ではあるが、これではやはり “風” というイメージがしない。 以下余談であるが、遠別の春（4-6月）の風向は、南南西と東南東が卓越している。初冬（10-12月）になると、西と東南東になる。一年を通じて東南東のことが多く、これは遠別川の方向が関係しているんだろうな、と思いつつ、 AMEDAS 遠別観測地点をチェックしてみた（中心の ＋ 位置）。 すると確かに遠別川の脇ではあるが、すぐ南は標高差は50m程度だが、崖の下である。どおりで南南東の風がほとんどないわけだ。こんな、明らかに地形に左右されるようなところに観測点があって、良いのだろうか？ それとも、遠別の風向はそれを加味してデータを見ろ、ということなのだろうか？ |