目次
はじめに
データを入手する
データの構成について
フォーマット
フォーマット例
メッシュ番号
3次メッシュデータを読みこんでみよう
緯度経度の範囲
メッシュを構成するグリッドセルの配列数
標高データ
データの並べ方
グリッドセルの開始点
GeoTiff変換しよう
変換用コード
GeoTiff変換用の関数
できたものを見てみよう
おわりに
参考
はじめに
やっほー、HerniaBabyです。
みんなはこのWebセミナーを見たことがありますか？
RoadRunner:日本の空中写真・標高データ活用法
けっこう革新的で素晴らしい動画なんだけどひとつだけ気になるところがあります。
それはこれ
え？！なんで外部ツール使うん？！
Mapping Toolboxあるやろ！
てなわけで、今回は国土地理院の数値標高データ(.xml)を自分でGeoTiff画像にします。
(このためにMapping Toolbox買いました。)
データを入手する
数値標高モデルについては利用規約に基づき基盤地図情報ダウンロードサイトから好きな場所をダウンロードしてください。
今回の題材として、メッシュ番号5841-46を選択して1mメッシュの航空レーザ測量データを使用します。
場所は以下の通りです。
出典：国土地理院の画像について
図は以下を基に承認申請しておりません。必要な場合は申請します。
その際は画像を削除する可能性があります。
Q1-11． Webサイトに地図を挿入して利用したいが申請は不要か？
ウェブサイトへの地図の挿入（貼り付け）は、書籍に地図を挿入する場合と同等と見なして、原則として申請不要とします。ただし、挿入した図をクリックした後に別窓が開き単体の地図が表示されるものや、ウェブサイトのメインコンテンツが地図である場合を除きます。(承認申請Q&Aより)
※2023年11月30日(木)より、最新の航空レーザ測量データを基に作成した１ｍメッシュデータが提供されるようになりました。
いやあ、すごいっすね！詳細リンクはこちら！
データの構成について
データの構成に関しては以下を参考にしました。
フォーマット
ファイル名とそれぞれの意味は以下のようになっています。
FG-GML-pppp-qq-rr-DEM5X-yyyymmdd.xml
フォーマット例
FG-GML-5841-46-00-DEM1A-20230929.xml
※今回は1mメッシュの航空レーザ測量なので1Aになっています。
メッシュ番号
メッシュ番号については地域メッシュ統計(総務省統計局)がわかりやすいです。
本記事では3次メッシュを1つずつGeoTiffデータにするというものです。
3次メッシュデータを読みこんでみよう
まずはひとつデータを読み込んでみましょうか。
xmlファイルは readstruct から読み込むことができます。
clear,clc;
file_name = 'FG-GML-5841-46-00-DEM1A-20230929.xml';
T = readstruct(file_name,'FileType','xml');
データ自体は T.DEM にあります。中身を見てみましょう。
緯度経度の範囲
lowerCornerで北西端、upperCornerで南東端の緯度経度をとってきます。
ここから範囲を決めることができるわけですね。
disp(T.DEM.coverage.gml_boundedBy.gml_Envelope)
srsNameAttribute: "fguuid:jgd2011.bl"
gml_lowerCorner: "39 141.75"
gml_upperCorner: "39.008333333 141.7625"
メッシュを構成するグリッドセルの配列数
GridEnvelopeに入ってます。セルを何分割して標高データを入れてるかになります。
low属性は北西端を示す番号です。常に(0,0)です。
high属性は南東端を示す番号です。5mメッシュは(224,149)、1m & 10mメッシュは(1124,749)固定です。
disp(T.DEM.coverage.gml_gridDomain.gml_Grid.gml_limits.gml_GridEnvelope)
gml_low: "0 0"
gml_high: "1124 749"
標高データ
tupleListに入っています。splitlines で改行部分で分割し一部見てみましょう。
水面や標高データがない場合は、標高の値が -9999 になります。
Surf = splitlines(T.DEM.coverage.gml_rangeSet.gml_DataBlock.gml_tupleList);
disp(Surf(1:5))
"海水面,-9999."
"海水面,-9999."
"海水面,-9999."
"海水面,-9999."
"海水面,-9999."
データの並べ方
orderAttributeに入ってます。データをどういう風に並べるか決めています。
以下引用(基盤地図情報~数値標高モデルを利用するより原文ママ)
orderの"+x"はxが正の方向(西→東)であることを、"-y"はyは負の方向(北→南)であることを示します。
北西端から南東端に向けて次のような並びで一列にしたデータがtupleListタグ内に配置されます。
だいたい+x-yなのであまり気にしてません。
国土地理院以外でやるなら見ておく必要がありそうです。
disp(T.DEM.coverage.gml_coverageFunction.gml_GridFunction.gml_sequenceRule.orderAttribute)
+x-y
グリッドセルの開始点
startPointに入ってます。
こいつが本当に厄介です。(0, 0)からスタートしない場合もあります！
実際のデータ数が(low + 1)×(high + 1)に必ずしもならないことに注意です！
disp(T.DEM.coverage.gml_coverageFunction.gml_GridFunction.gml_startPoint)
0 414
今回も(0, 414)からスタートしてますね。
つまり(0, 0) ~ (0, 413)の標高情報はここにはありません。
startPoint がゼロスタートでない場合はNaNで埋めることにします。
埋め方について、ツールによっては-9999mをNaNとして認識するように設定するものもあるらしいのですがMATLABはそうじゃないそうです(MathWorks情報)。
GeoTiff変換しよう
先の情報がわかればGeoTiff変換は簡単です。
georefcells に緯度, 経度, 標高データを入れて、geotiffwrite で書き出すだけです。
変換用コード
clear,clc,close all;
path = 'FG-GML-5841-46-DEM1A';
tmpファイルを作成しGeoTiff化したファイルを入れます。
if ~exist(fullfile(path,'tmp'))
mkdir(path, 'tmp')
end
ひとつずつ変換します。ファイルの数だけループします。
流石に重たいので parfor で並列処理します。
fnames = dir(fullfile(path,"*.xml"));
path = fnames(1).folder;
file = {fnames.name}';
fsize = length(file);
parfor ii = 1:fsize
file_name = fullfile(path,file{ii});
MyGeoTIFFConvert(file_name);
end
終わったらメンションします。
disp('嗚呼…GeoTIFF変換完了…')
↓
嗚呼…GeoTIFF変換完了…
嗚呼…GeoTIFF変換完了…
ふざけないと死んじゃうのは私の悪い癖ですね
GeoTiff変換用の関数
function [] = MyGeoTIFFConvert(file_name)
% 参考：https://qiita.com/tobira-code/items/43a23362f356198adce2
% 出力用に名前を保存
[fpath, name, ~] = fileparts(file_name);
name = string(strsplit(name,'-'));
name = strjoin(name(3:5),'');
% 緯度軽度：これをもとにGeoTiffに変換
T = readstruct(file_name,'FileType','xml');
T = T.DEM.coverage;
LatLon = struct2array(T.gml_boundedBy.gml_Envelope)';
LatLon = LatLon(2:end);
LatLon = double(split(LatLon))';
% DataBlock：各グリッドセルの標高データ
Raster = splitlines(T.gml_rangeSet.gml_DataBlock.gml_tupleList);
Raster = str2double(extractAfter(Raster,','));
% GridEnvelope：メッシュを構成するグリッドセルの配列数
GridEnvelope = struct2array(T.gml_gridDomain.gml_Grid.gml_limits.gml_GridEnvelope)';
GridEnvelope = double(split(GridEnvelope));
% startPoint：グリッドセルの開始点を指定
StartPoint = str2double(strsplit(T.gml_coverageFunction.gml_GridFunction.gml_startPoint))+1;
% サイズ
RasterSize = diff(GridEnvelope) + 1;
RasterLength = prod(RasterSize);
% startPointが0スタートでない場合、NaNで埋める
% MATLABではNaNで出力するのがいいらしい(-9999mをNaNとして認識できない)
if length(Raster) == RasterLength
idx = isnan(Raster);
% Raster(idx) = -9999;
Raster(idx) = NaN;
RasterMap = reshape(Raster,RasterSize)';
else
StartPoint = prod(StartPoint);
tmp = nan(RasterLength,1);
tmp(StartPoint:StartPoint+length(Raster)-1) = Raster;
idx = isnan(tmp);
% tmp(idx) = -9999;
tmp(idx) = NaN;
RasterMap = reshape(tmp,RasterSize)';
end
% GeoTiffファイルとして書き出し
R = georefcells(LatLon(1,:), LatLon(2,:), flip(RasterSize),"ColumnsStartFrom","north");
geotiffwrite(fullfile(fpath,'tmp',strjoin([name,'.tif'],'')),single(RasterMap),R,"CoordRefSysCode",'EPSG:6668')
end
できたものを見てみよう
さすがに RoadRunner で確認なんてできませんので、フリーソフト QGIS で確認します。
そう、これでは2次メッシュがきれいに繋がらないんですよねー
欠損があったり境界部がきれいにつながってなかったり…
おわりに
今回は国土地理院のDEMデータについて中身をだいたい理解してGeoTiffファイルにしました。
次回はモザイク結合で2次メッシュデータをきれいに繋げます。