image_print

Table of Contents

9 JGD2024の取扱とGrassのLocation作成

 ここまでトライしてきて,大きな問題にぶつかった。JGD2024の座標参照系が未定なのである。国土地理院はEPSG

(Google AI回答: EPSG(European Petroleum Survey Group)データセットの登録組織は、現在、国際石油・ガス生産者協会(IOGP)のジオマティクス委員会によって管理されています。 元々は、ヨーロッパの石油会社に所属する測量専門家による非公式な科学団体であるEuropean Petroleum Survey Group (EPSG) がデータセットを作成・維持していましたが、2005年に同組織は解散し、その活動はIOGPに引き継がれました。 現在もデータセットの名称として「EPSG」というブランド名が残されていますが、管理運営はIOGPが行っており、測地系、座標参照系、変換パラメーターなどの情報を含むこのデータセットを維持・公開しています。データセットは、EPSG Registryのウェブサイトで無償で公開されています。 )

に申請してはいるだろうけど承認が遅いのか。不思議な気がする。GrassGISはEPSGが決まらないと使えない(,と思ってきた)。

 JGD2011とJGD2024の違いは,日本測地系2024(JGD2024)への移行で何が変わる?建設・測量現場への影響と対策 2025年4月18日 掲載

に記されているように,平面直角座標系の(easting, northing)座標値はそのままで,高度情報が最大 +/- 60cmほど変化する。当方の使用目的からすると,気にしなくて良い,というのが,Perplexityさんだが,これは嘘のたぐいだ。

 なお,平面直角座標系座標値の表現の順序は,(easting, northing),,か,(northing, easting)か。Google AIに聞いてみた。

日本の平面直角座標系における座標値の表現の順序は、原則として(northing, easting) 、すなわち (X, Y) の順です。 一般的な数学の座標系(右手系)では,𝑋 が東西(easting)、𝑌 が南北(northing)を表しますが、日本の測量における平面直角座標系では以下の通り定義されています: 
X軸(northing): 南北方向。原点を通る子午線の方向を正(北)とする。
Y軸(easting): 東西方向。原点においてX軸に直交する方向を正(東)とする。 
 したがって、座標値を記述する際は、まず南北を示す X値(northing)、次に東西を示す Y値(easting) の順で並べるのが標準的です。 
 なお、国際的な地理情報システム(GIS)やUTM座標系では、数学的な慣習に合わせて (easting, northing) の順で扱う場合もあるため、利用するソフトウェアやデータの仕様を確認することが重要です。

 今後,ぼくの時間は少ないのは確かであるが,わざわざ,現実とずれた座標系を使いたくはない。そこで,EPSGが決まるまでは,GrassGISは使わない,ことにする。

 QGISは,on the flyという機能があって,QGISに質問した。その回答を次に。

——— 引用〜

「不明なCRS」のままGSI数値地図を表示するのは実用的問題ありません。

問題なしの理由

✅ QGIS on the fly → 形状位置自動調整(内部座標一致)
✅ GRASS GIS → DEM形状自動認識(CRS無視で解析)
✅ KyotoOsaka_extent → 形状正常重畳(位置正確)

実用確認ポイント

■ マップ表示 → GSIデータ位置正確 ✓
■ 縮尺・座標値 → 正常表示 ✓
■ レイヤ重畳 → 凸包・等高線位置一致 ✓

推奨運用

プロジェクトCRS → 「不明なCRS」維持
各レイヤCRS → GSIデータ自動認識
GRASS Location → DEM自動設定(完璧)

形状表示・解析精度に影響なし![1]

「不明なCRS」そのままQuickDEM4JP → GRASS解析へ直行!

GSIデータ位置正常なら完璧進行!

引用〜おわり ———

 Grassに関する記述はやってみなとわからない。図11はプロジェクトのプロパティの座標参照系の情報だ。JGD2024の投影法などの情報がQGISによって正しく認識されている。GRASS利用なしで,QGISの機能の範囲で使用したいと思う。

図11 CRSなし(または未知)

9’追記 2026年1月10日

 流域の面積を求めるには経緯度座標系のままではできない。平面直角座標系を選択せざるを得ない。QiitaのEPSG を紐解く (2) では,「日本の場合、一つ前の改定が東日本大震災の影響で同時に測地系と鉛直基準系(水準原点の標高値)両方であったため、どちらも「日本測地系2011」とみなされており、さほど問題にはなりませんでした。個人的には、 Vertical Datum に対し「日本測地系2024 / JGD2024」という名称を使ったり、改定していない Geodetic Datum に対し「日本測地系2024 / JGD2024」とするのはやめて欲しかったなあ。と思います。」,とあり,ぼくは”nodding yes”であった。

 QGISのon the fly機能はインポートしたコンテンツの来歴に変更を加えるものでもないし,こだわりをおさえて,EPSG:6674(JGD2011 / Japan Plane Rectangular CS VI),を使おうと思った。ところがEPSG:10167(JGD2011 / Japan Plane Rectangular CS VI + JGD2011 (vertical) height)というのもある。この実質的な違いがわからない。

 PerplexityはEPSG:10167がJGD2024対応だというが嘘だ。

図11′ EPSG6674と10167比較

 Google Chromeに聞いてみた。

——— 引用〜

QGISでJGD2024の平面直角座標系Ⅵを扱う場合、既存のプロジェクトのCRS(座標参照系)をJGD2011のEPSG:6674に設定し直すことで対応できます。QGISには「オンザフライCRS変換」機能があるため、プロジェクトのCRSを変更しても、インポート済みのレイヤーは自動的に新しいCRSに合わせて表示されます。 (ぼくのメモ:嘘)

JGD2024はJGD2011とほぼ同等の測地系であるため、現状はEPSG:6674で問題なく扱えます。 

QGISでの設定方法

1. プロジェクトのCRSを変更する手順

既存のプロジェクトの座標系をデフォルトの経緯度座標系(おそらくEPSG:4326など)から変更する手順です。

  1. QGISを開き、該当のプロジェクトをロードします。
  2. メニューバーから[プロジェクト] > [プロパティ] を選択します。
  3. [プロジェクトプロパティ]ウィンドウの左側パネルで [CRS] タブを選択します。
  4. 上部にあるフィルタ(検索)ボックスに 6674 または JGD2011 Japan Plane Rectangular CS VI と入力して検索します。
  5. 検索結果から JGD2011 / Japan Plane Rectangular CS VI (EPSG:6674) を選択します。
  6. [適用] ボタンをクリックし、次に [OK] ボタンをクリックしてウィンドウを閉じます。 

プロジェクトの座標系が変更され、インポート済みのレイヤーが新しい座標系に自動的に再投影(オンザフライ変換)されて表示されます。 

2. 新しいプロジェクトのデフォルトCRSを設定する(推奨) 

今後新しいプロジェクトを作成する際に、最初からJGD2024(EPSG:6674)を使いたい場合は、QGISのデフォルト設定を変更しておくと便利です。

  1. メニューバーから [設定] > [オプション] を選択します。
  2. [オプション]ウィンドウの左側パネルで [CRS設定] を選択します。
  3. [新しいプロジェクト作成時のデフォルトCRS] の項目で、JGD2011 / Japan Plane Rectangular CS VI (EPSG:6674) を検索して選択します。
  4. [新しいレイヤ作成時またはCRS未定義レイヤ読み込み時のデフォルトCRS] の項目も同様に設定することをお勧めします。
  5. [OK] をクリックして設定を保存します。 

これにより、新しいプロジェクトやCRSが未定義のデータを読み込む際に、自動的にEPSG:6674が適用されるようになります。

注意点

  • データソースは変更されない: この設定はQGIS上での表示方法(プロジェクトの解釈)を変更するものであり、元のデータファイル自体の座標系が書き換わるわけではありません。
  • JGD2024とJGD2011: JGD2024は最新の測地成果ですが、GISソフトウェアのEPSGレジストリではまだ正式な専用コードとして広く普及していない場合があります(2025年4月時点)。現状では、JGD2011のコード(EPSG:6674)が実質的な標準として利用されています。 

引用〜おわり ———

 Google Chrome (フリー)の検索結果は,Perplexity Pro よりも的確だ。Perplexity Education Proへの移行がずっと手間取っている。verifyのメールを送ったというが,何度もクレームをつけているが,届かない。そして,Perplexityの回答の不正確さはかなり致命的なものである。

 さて,上記のようにEPSGを変更した。ほぼ全レイヤーの右端に,変なアイコン(赤色のメガネの上に横一)が付いている。この上にマウス先をおくと,「EPSG:4326 – WGS84は,動的CRSですが,基準時(Epoch)が設定されていません。定義が曖昧で精度が限られます」,というメッセージが出る。それで,たとえば10mメッシュDEMレイヤーのプロパティで,EPSG:6674に変更すると,10mメッシュDEMレイヤー,を選択しているにもかかわらず,DEM雲が見えなくなる。そこでこのレイヤーを削除して,あらたにインポートすると,見えるようになった。

 そして,QGISのファイルを保存,終了して,新たに立ち上げると,また,見えなくなった。最下行のステータスバーは,EPSG:6674と表示されている。

 これに対処するには,根本的にこれまでを見直す必要がでてきたので,このページの現在の最新の章に移動して,書き込むことにしたい。「17 経緯度座標系から平面直角座標系への移行について

10 DEMマージ(QuickDEM4JP)

 リージョン region を設定するのに, 8.2のように,等高線シェープを使おうとしたが,QGISを再起動するとそのレイヤーが機能しないというメッセージが出て,結局削除して,GSIの等高線シェープファイルだけが残った。これは線分から構成されるもので,ぼくは,もともと,これをただ作業の際に目安として利用するためにダウンロードしたのであって,線分でも何も問題はない。

 どうもQGISはリージョンという意識をさせる必要がないようなのであるが,DEMはベクトルとは違って,リージョンの杭打ちのような効果があるように思う。

 QuickDEM4JPは上記のように,すでにインストールした。→25枚XML選択→出力「merged_dem.tif」(自動JGD2024認識)。Processing→「GDAL」→「ラスタをマージ」→入力複数GeoTIFF→範囲.shp指定→単一「kameoka_dem_merged.tif」出力,と,Perplexityは言っているが,実行できなかった。

 種々実行したが,成功したのは次のものである。

図12 5m メッシュDEMのマージ完成

 どのように成功したのか。QuickDEM4JPの開発者のWebページがあるがそこでは単一メッシュしか実行されていない。いわく,「zipファイルを解凍してxmlファイルにしてそれを入力すれば,実行できる。複数のxmlファイルを一つのフォルダに入れて,そのフォルダごとインプットすることもできる」,とあるが,現在流布しているこのプラグインにはその機能はない。

 GSIから選択したDEMファイルはまとめて一つのzipファイルとしてダウンロードできる,というか,普通そうするだろう。そのため,ダウンロードした一つのzipファイルには,さらにzipファイル群があって,このまとまった一つのzipファイル(Packed)に対しては,このプラグインは実行できなかった。

 QuickDEM4JPは,1回解凍した場合にxmlファイルで構成されるzipファイルか,xmlファイルそのものに対して。実行できるのである。

 QuickDEM4JPでの実行の際に,EPSGを指定する必要があり,デフォルトのEPSG4326(-WGS84),を選択せざるを得ない。これは経緯度座標系であるが,出力された座標値は幸い,平面直角座標系Ⅵのものとなっている。

 QuickDEM4JPで実行できるファイルは,一回圧縮のzipか,xmそのものである,まあ当然,zipを使うことになるだろう。

 10m meshのまとまった一つのzipファイルは25個のzipファイルから構成される。この一つ一つのzipファイルの中には,それぞれ一つのxmlファイルしかない。10mメッシュではxmlを選んでも,zipを選んでもQuickDEM4JPの実行回数は25回であり,zipファイルを選ぶメリットはない。

 5m meshは,6zipからなる,それぞれのzipを解凍すると99個ものxmlファイルがあるので,もちろん,QuickDEM4JPで実行する対象はzipということになる。

 QuickDEM4JPで実行するとgeotiffが作成される。複数のgeotiffを結合mergeするとGISで利用が可能になる。 

  1. プロセシングツールボックスの検索窓で「結合」と入力検索して,[ラスタ] > [その他] > [結合] アルゴリズム(GDALのgdal_mergeに相当)を見つける。
  2. [結合]ウィンドウを開く。
    • 画面左の[入力レイヤ] に表示されている個々のgeotiffをすべて選ぶ。5235-43_5m, ………,5235-55_5mの6枚。
    • [出力ファイル] で、結合後の最終的なGeoTIFFファイルの保存先と名前を指定する。43to55merged.tif
    • [実行] ボタンをクリックして結合する。

この手順により、複数の基盤高度メッシュデータから一つの統合されたGeoTIFFファイルを作成できるのである。このようにしてできたのが,図12のマージした5mメッシュDEMである。

 前述の10mメッシュも実行したい。25個のzipファイルにはそれぞれ一つのxmlファイルがあり,QuickDEM4JPが一つのフォルダーを選ぶことができるというフェイクですでに解凍しているのでこれを一つのzipファイルに圧縮して,QuickDEM4JPを使ってgeotif化を実行してみよう。

 GIS_DataのInput_Data > GSI_5235-325o76_10mDEMフォルダー内の25個のxmlファイルを,新たに作成した一つのフォルダー 25xmls_10mDEMにまとめて(360MB),7-Zipで, ”25xmls_10mDEM.zip”(42MB)に圧縮。360MBが42MBに。10mDEM32to76merged,と名付けた。

 図13には方位とスケールを追加した。桂川全流域が表示されている。結構,隣接の流域と谷中分水界をなしている。

図13 10m-mesh DEM

11 流域解析 r.slope.aspect

 さて,10mメッシュで,流域区分図(ラスター)と水系図(ベクトル)を作成したい。ぼくの教材を次に。

地形と流域の分析 Nov. 30, 2010。15年前のGrassGIS教材である。

11.1 g.regionの取り扱い

Perplexityさんのご指導

——— 引用〜g.region

QGISのProcessing ToolboxでGRASSコマンドを実行する場合、g.regionは独立したツールとして存在せず、各GRASS r.*モジュール(r.slope.aspectなど)の入力パラメータで自動的にregionが管理されます。

region管理方法

「GRASS Region」または「GRASS_REGION_PARAMETER」パラメータから入力レイヤを選択すると、g.regionが内部的に呼び出され、そのラスタ範囲に自動調整されます。
「GRASSの範囲を標高ラスタに合わせない」オプションでこの自動調整を無効化し、カレントregionを使用できます。メモ:換言すれば,DEMに合わせることで g.region rasterDEMの効果があるということ。
QGIS GRASSプロバイダ(Processing > Toolbox > GRASS > Raster)内で直接g.regionを探さず、各ツールのパラメータを確認してください。

引用〜おわり ———

11.2 Warning: Concurrent mapset locking is not supported on Windows.

Perplexityさんのご指導

——— 引用〜

GRASSがマップセットの同時ロック機構を使用しようとするが、Windowsファイルシステムが対応していないため発生します。
QGIS Processingで一時GRASS LOCATIONを作成する際に頻出で、他のGRASSツール(r.watershed, r.fill.dir)でも同様。

引用〜おわり ———

11.3 zscale

Perplexityさんのご指導

——— 引用〜

QGIS Processing Toolboxのr.slope.aspectで「標高をメートルに換算する係数」がzscaleパラメータに対応します。この欄に「1.0」を入力してエラーを解消してください。
設定方法
「標高をメートルに換算する係数」テキストボックスに1.0を明示的に入力します。空欄だと自動で0と認識されエラーになります。
DEMがメートル単位なら1.0、

引用〜おわり ———

11.4 SetColor Table()

Perplexityさんのご指導

——— 引用〜

“ERROR SetColor Table() only supported for Byte or Unit 16 bands in TIFF format”は、GRASSのr.slope.aspectが浮動小数点(FCELL)出力をGeoTIFF保存時にカラーテーブル適用を試み、GDALがByte/UInt16のみ対応のため発生する警告です。
警告の原因
r.slope.aspectのaspect/slope出力はFCELL(Float32)で、カラーテーブルは整数型(Byte/UInt16)専用です。TIFF出力時に自動カラー割り当てが失敗。
処理自体は成功し、出力ラスタは正常に生成されます。QGISレイヤパネルで確認可能。
対処方法
• 無視して使用: 警告をログで確認後、レイヤを通常通り表示・解析
• Properties > Symbology > Singleband pseudocolorで手動グラデーション設定
• 出力形式変更: Processing > r.out.gdal > output_format=AAIGrid(ASCII)で回避
実務では無視してOKで、数値解析に影響しません。

引用〜おわり ———

11.5 r.slope.aspect

 これを実行する過程で,11.1〜11.4の課題が出てきたのであるが,QGISではaspect生成が難しい。aspectの作成にはGrassGISそのもので実行した方がよい。このコマンド処理についての覚書を次に。

FCELL: Float 32はデフォルトでこのまま。1階偏微分E-W, N-Sの計算結果だけは使える。

12. r.watershed

 引き続き,地形と流域の分析 Nov. 30, 2010。15年前のGrassGIS教材参照。

 Grassの実行名称が十五年前とはかなり変更されていて,理解しずらい。

 本コマンドで出力された各レイヤーの対応は以下の通りになる。accumulation(集水量)、drainage(流向)、basin(流域ID)はそれぞれ,英語名:accumulationdrainagebasin,に対応する。

12.1 出力対応表

出力名(日本語)英語名説明
流域ラベルbasin各流域にユニークID(偶数)を付与。NULLは不完全流域。
傾斜長勾配係数length_slopeRUSLEのLS因子(傾斜長×勾配)。侵食予測用。
排水先セルの数accumulation上流集水セル数(負値=外部流入)。流出量指標。
河川セグメントstreambasin ID対応の河川線。r.thin後ベクトル化推奨。
傾斜勾配係数slope_steepnessRUSLEのS因子(傾斜勾配)。
半流域half_basinbasinを左右半分(右偶数、左奇数)。
Stream Power IndexspiSPI = α × tan(β)。河川浸食力。
排水方向drainage流向(1-8、NEから反時計回り)。東=1、北東=2など。
Topographic indextciTCI = ln(α / tan(β))。土壌水分指標。

12.2 「外部流域(exterior basin)の最小サイズ」の指定

 このコマンドでは,何らかの専門性がなければ,ラスターファイルは,標高の10mDEMだけを指定する。

 そして他では,唯一,「外部流域(exterior basin)の最小サイズ」を指定する必要がある。これは,従来の「最小流域サイズ(Minimum size of exterior watershed basin)」に対応する。

 一般的アドバイスとしては,10m DEMの場合、5000(0.5km²相当)が推奨。値が小さいと計算時間激増・細かすぎる流域になるとされる。外部流域(exterior basin)の欄に,5000と入力せよということだ。

 √0.5km² = 0.707kmとなる。5000はセルの個数を意味し,一辺だと,√5000 = 70.7 個となる。長さにすると,70.7 x 10m = 707mだ。ラスター1個,オプション値一つを入力した結果,いくつかのラスターレイヤーが出力される。図14はdrainageだ。aspectに類似している。排水方向というよりも,ぼくには斜面方向がよくわかっていい。

排水方向drainage流向(1-8、NEから反時計回り)。東=1、北東=2など。
図14 排水方向の出力 = drainage

 次の図15は,streamを示したものである。実は,図14にもこのstreamは表示しており,drainageに見える谷地形によく一致している。

河川セグメントstreambasin ID対応の河川線。r.thin後ベクトル化推奨。
図15 河川セグメント=stream

 

 r.watershed実行で出力されたラスターはこの二つ,つまり,「河川セグメント」stream と「排水方向」drainage だけ残すことにした。ただ,「排水先セルの数」accumulation は桂川流域の限定などのために使う可能性があって,残すことにした。

 r.watershedを実行してから,必要なレイヤーをexportすると色も失われ,解像度もかなり低下する(ような?)ので,実行時に出力した方がいいように思うので,再度,再び実行してみよう。外部流域の最小サイズを教材でやったように10,000としてみたが変わらないようなので,5000を維持することにした。実行時に,ファイルを保存した。排出先 accumulation,河川セグメントstream,排水方向drainage,流域ラベルbasin。

image_print