日本の標高データをMATLABで活用 - GeoTiff変換 -



## 目次
[はじめに](#はじめに)
[データを入手する](#データを入手する)
[データの構成について](#データの構成について)
&emsp;[フォーマット](#フォーマット)
&emsp;&emsp;[フォーマット例](#フォーマット例)
&emsp;[メッシュ番号](#メッシュ番号)
[3次メッシュデータを読みこんでみよう](#3次メッシュデータを読みこんでみよう)
&emsp;[緯度経度の範囲](#緯度経度の範囲)
&emsp;[メッシュを構成するグリッドセルの配列数](#メッシュを構成するグリッドセルの配列数)
&emsp;[標高データ](#標高データ)
&emsp;[データの並べ方](#データの並べ方)
&emsp;[グリッドセルの開始点](#グリッドセルの開始点)
[GeoTiff変換しよう](#GeoTiff変換しよう)
&emsp;[変換用コード](#変換用コード)
&emsp;[GeoTiff変換用の関数](#GeoTiff変換用の関数)
[できたものを見てみよう](#できたものを見てみよう)
[おわりに](#おわりに)
[参考](#参考)
 
<a name="endToc"></a>
<a name="H_763E9CC4"></a>
# はじめに

やっほー、HerniaBabyです。
みんなはこのWebセミナーを見たことがありますか？


[RoadRunner:日本の空中写真・標高データ活用法](https://jp.mathworks.com/videos/roadrunner-how-to-import-japanese-aerial-and-elevation-images-1629895701730.html)
けっこう革新的で素晴らしい動画なんだけどひとつだけ気になるところがあります。


それはこれ
![image_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/781521/3d112e90-1f22-cb39-6e39-eb8d1606de3a.png)




**え？！なんで外部ツール使うん？！**


**Mapping Toolboxあるやろ！**


てなわけで、今回は国土地理院の数値標高データ(.xml)を自分でGeoTiff画像にします。
(このためにMapping Toolbox買いました。)

<a name="H_1826D6CA"></a>
# データを入手する

数値標高モデルについては利用規約に基づき[基盤地図情報ダウンロードサイト](https://fgd.gsi.go.jp/download/menu.php)から好きな場所をダウンロードしてください。


今回の題材として、メッシュ番号5841-46を選択して1mメッシュの航空レーザ測量データを使用します。
場所は以下の通りです。
![image_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/781521/cd446d7c-2a51-cf48-79f6-998d54137044.png)

::: note warn
出典：国土地理院の画像について
図は以下を基に承認申請しておりません。必要な場合は申請します。
その際は画像を削除する可能性があります。

***Q1-11．　Webサイトに地図を挿入して利用したいが申請は不要か？***
ウェブサイトへの地図の挿入（貼り付け）は、書籍に地図を挿入する場合と同等と見なして、原則として申請不要とします。ただし、挿入した図をクリックした後に別窓が開き単体の地図が表示されるものや、ウェブサイトのメインコンテンツが地図である場合を除きます。[(承認申請Q&Aより)](https://www.gsi.go.jp/LAW/2930-qa.html)
:::

※2023年11月30日(木)より、最新の航空レーザ測量データを基に作成した**１ｍメッシュデータ**が提供されるようになりました。


いやあ、すごいっすね！詳細リンクは[**こちら**](https://www.gsi.go.jp/gazochosa/gazochosa61002.html)！

<a name="H_062B9C3B"></a>
# データの構成について

データの構成に関しては以下を参考にしました。


[基盤地図情報~数値標高モデルを利用する](https://qiita.com/tobira-code/items/43a23362f356198adce2)

<a name="H_F9A861EB"></a>
## フォーマット

ファイル名とそれぞれの意味は以下のようになっています。

FG-GML-pppp-qq-rr-DEM5X-yyyymmdd.xml

-  pppp:1次メッシュ番号 
-  qq: 2次メッシュ番号 
-  rr: 3次メッシュ番号 
-  X: 測量方法(A = レーザー測量, B = 写真測量) 
-  yyyymmdd:  年月日 
<a name="H_FC99A7FE"></a>
### フォーマット例

FG-GML-5841-46-00-DEM1A-20230929.xml
※今回は1mメッシュの航空レーザ測量なので1Aになっています。

<a name="H_7A815EA1"></a>
## メッシュ番号

メッシュ番号については[地域メッシュ統計(総務省統計局)](https://www.stat.go.jp/data/mesh/index.htm)がわかりやすいです。


[第１章 地域メッシュ統計の特質・沿革 (9頁目)](https://www.stat.go.jp/data/mesh/pdf/gaiyo1.pdf#page=9) 
![image_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/781521/6378b886-c90e-e663-0f84-b71b3f3db636.png)




本記事では**3次メッシュを1つずつGeoTiffデータにする**というものです。

<a name="H_36AC24A2"></a>
# 3次メッシュデータを読みこんでみよう

まずはひとつデータを読み込んでみましょうか。


xmlファイルは [readstruct](https://jp.mathworks.com/help/matlab/ref/readstruct.html) から読み込むことができます。

```matlab
clear,clc;
file_name = 'FG-GML-5841-46-00-DEM1A-20230929.xml';
T = readstruct(file_name,'FileType','xml');
```

データ自体は T.DEM にあります。中身を見てみましょう。

<a name="H_FEFCED12"></a>
## 緯度経度の範囲

<samp>lowerCornerで北西端、upperCornerで南東端の緯度経度をとってきます。</samp>


<samp>ここから範囲を決めることができるわけですね。</samp>

```matlab
disp(T.DEM.coverage.gml_boundedBy.gml_Envelope)
```

```TextOutput
    srsNameAttribute: "fguuid:jgd2011.bl"
     gml_lowerCorner: "39 141.75"
     gml_upperCorner: "39.008333333 141.7625"
```
<a name="H_913CFDEE"></a>
## メッシュを構成するグリッドセルの配列数

GridEnvelopeに入ってます。セルを何分割して標高データを入れてるかになります。


low属性は北西端を示す番号です。常に(0,0)です。


high属性は南東端を示す番号です。5mメッシュは(224,149)、1m & 10mメッシュは(1124,749)固定です。

```matlab
disp(T.DEM.coverage.gml_gridDomain.gml_Grid.gml_limits.gml_GridEnvelope)
```

```TextOutput
     gml_low: "0 0"
    gml_high: "1124 749"
```
<a name="H_FF798698"></a>
## 標高データ

tupleListに入っています。[splitlines](https://jp.mathworks.com/help/matlab/ref/splitlines.html) で改行部分で分割し一部見てみましょう。


水面や標高データがない場合は、標高の値が -9999 になります。

```matlab
Surf = splitlines(T.DEM.coverage.gml_rangeSet.gml_DataBlock.gml_tupleList);
disp(Surf(1:5))
```

```TextOutput
    "海水面,-9999."
    "海水面,-9999."
    "海水面,-9999."
    "海水面,-9999."
    "海水面,-9999."
```
<a name="H_9D2E0B5C"></a>
## データの並べ方

orderAttributeに入ってます。データをどういう風に並べるか決めています。


以下引用([基盤地図情報~数値標高モデルを利用する](https://qiita.com/tobira-code/items/43a23362f356198adce2)より原文ママ)


>*orderの"+x"はxが正の方向(西→東)であることを、"-y"はyは負の方向(北→南)であることを示します。* 
 *北西端から南東端に向けて次のような並びで一列にしたデータがtupleListタグ内に配置されます。* 

![image_3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/781521/99906134-13e2-6b64-1f4c-4238528d609a.png)



だいたい+x-yなのであまり気にしてません。


国土地理院以外でやるなら見ておく必要がありそうです。

```matlab
disp(T.DEM.coverage.gml_coverageFunction.gml_GridFunction.gml_sequenceRule.orderAttribute)
```

```TextOutput
+x-y
```
<a name="H_BD65A91F"></a>
## グリッドセルの開始点

startPointに入ってます。

:::note warn
**こいつが本当に厄介です。(0, 0)からスタートしない場合もあります！**


**実際のデータ数が(low + 1)×(high + 1)に必ずしもならないことに注意です！**
:::
```matlab
disp(T.DEM.coverage.gml_coverageFunction.gml_GridFunction.gml_startPoint)
```

```TextOutput
0 414
```

今回も(0, 414)からスタートしてますね。

つまり(0, 0) ~ (0, 413)の標高情報はここにはありません。

startPoint がゼロスタートでない場合はNaNで埋めることにします。

:::note warn
埋め方について、ツールによっては-9999mをNaNとして認識するように設定するものもあるらしいのですがMATLABはそうじゃないそうです(MathWorks情報)。
:::
<a name="H_84B250C1"></a>
# GeoTiff変換しよう

先の情報がわかればGeoTiff変換は簡単です。


[georefcells](https://jp.mathworks.com/help/map/ref/georefcells.html) に緯度, 経度, 標高データを入れて、[geotiffwrite](https://jp.mathworks.com/help/map/ref/geotiffwrite.html) で書き出すだけです。

<a name="H_E38DAD00"></a>
## 変換用コード
```matlab
clear,clc,close all;
path = 'FG-GML-5841-46-DEM1A';
```

tmpファイルを作成しGeoTiff化したファイルを入れます。

```matlab
if ~exist(fullfile(path,'tmp'))
    mkdir(path, 'tmp')
end
```

ひとつずつ変換します。ファイルの数だけループします。


流石に重たいので [parfor](https://jp.mathworks.com/help/matlab/ref/parfor.html) で並列処理します。

```matlab
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
```

終わったらメンションします。

```matlab
disp('嗚呼…GeoTIFF変換完了…')
```
↓
```TextOutput
嗚呼…GeoTIFF変換完了…
```

嗚呼…GeoTIFF変換完了…


![4494295i.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/781521/8ead5eea-1a26-3b81-258e-c26822f55a67.jpeg)


ふざけないと死んじゃうのは私の悪い癖ですね

<a name="H_28B1810D"></a>
## GeoTiff変換用の関数
```matlab
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
```
<a name="H_6CF1DCF0"></a>
# できたものを見てみよう

さすがに [RoadRunner](https://jp.mathworks.com/products/roadrunner.html) で確認なんてできませんので、フリーソフト [QGIS](https://qgis.org/ja/site/) で確認します。

![image_5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/781521/76bc0890-5522-2358-46da-f87f58593cbe.png)


そう、これでは2次メッシュがきれいに繋がらないんですよねー


欠損があったり境界部がきれいにつながってなかったり…

<a name="H_92B32A6B"></a>
# おわりに

今回は国土地理院のDEMデータについて中身をだいたい理解してGeoTiffファイルにしました。


次回はモザイク結合で2次メッシュデータをきれいに繋げます。

<a name="H_59C2ECC7"></a>
# 参考
-  [RoadRunner:日本の空中写真・標高データ活用法](https://jp.mathworks.com/videos/roadrunner-how-to-import-japanese-aerial-and-elevation-images-1629895701730.html) 
-  [総務省統計局：地域メッシュ統計](https://www.stat.go.jp/data/mesh/index.htm) 
-  [基盤地図情報~数値標高モデルを利用する](https://qiita.com/tobira-code/items/43a23362f356198adce2) 
-  [pythonで地域メッシュごとのポリゴンデータからGeoTIFFを作成する：国土数値情報平年メッシュデータを例に](https://computational-sediment-hyd.hatenablog.jp/entry/2022/11/19/000000)