日本の標高データをMATLABで活用 - モザイク結合 -



## 目次
[はじめに](#はじめに)
 &emsp;[今回やること](#今回やること)
[フォルダを読み込む](#フォルダを読み込む)
[2次メッシュに当てはめるための格子を作成する](#2次メッシュに当てはめるための格子を作成する)
&emsp;[3次メッシュ番号を抜き出す](#3次メッシュ番号を抜き出す)
&emsp;[3次メッシュデータを決まった場所に格納する](#3次メッシュデータを決まった場所に格納する)
&emsp;[データのない箇所を探す](#データのない箇所を探す)
&emsp;[データを読み込む](#データを読み込む)
&emsp;[データの欠損部分をNaNで埋める](#データの欠損部分をnanで埋める)
&emsp;[画像データを上下逆にする](#画像データを上下逆にする)
&emsp;[モザイク結合](#モザイク結合)
[出力](#出力)
[閑話休題：EPSG](#閑話休題epsg)
[出力結果](#出力結果)
[おまけ](#おまけ)
[参考](#参考)

 
<a name="endToc"></a>
<a name="H_05CF231C"></a>
# はじめに

<details><summary>暇な人向け挨拶はこちら</summary>

---
### やあみんな！俺だ！
#### シンボリックは数式。俺はこじき！それじゃあ！！

![やあみんな俺だ.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/781521/0f6e4042-2f4e-55de-fcd2-1442b8a3f087.png)
※この画像は載せるのをやめる予定でしたが強い要望のため掲載しました。

---

</details>

本件は[日本の標高データをMATLABで活用 - GeoTiff変換 -](https://qiita.com/Hernia_Baby/private/aa8f86f38e5f12f51c5b)の続きです。


先回はDEMデータをGeoTiffファイルにしました。


しかしきれいに繋がらなかったり、そもそも3次メッシュが一部欠損してたり…


これらを(0, 0) ~ (9,9)までの区画に当てはめて2次メッシュをきれいに結合します。


モザイク結合については[ラスターの統合 - モザイク](https://fme-casestudy.blogspot.com/2015/07/blog-post_11.html)を見てください。

<a name="H_E959D13E"></a>
## 今回やること

3次メッシュを100枚集めて結合し2次メッシュにします。


イメージとしては以下のように3次メッシュ5841-46-43を(4, 3)に当てはめていく感じです。

![image_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/781521/1edc5665-e6b1-1eba-ce6f-1e46ca35a9d1.png)



<a name="H_D23BC5B9"></a>
# フォルダを読み込む

先回作った5841-46-XXのメッシュでXXは(0, 0) ~ (9, 9)にそれぞれ当てはまる合計100このtifファイルがあるはずです。


まずは先回作ったtmpファイル内の中を見ていきましょう。

```matlab
clc,clear,close all;
path = 'FG-GML-5841-46-DEM1A\tmp';
fnames = dir(fullfile(path,"*.tif"));

path = fnames(1).folder;
file = {fnames.name}';
fsize =  length(file);
```

後は100枚のファイルを並べてモザイク結合を…

```matlab
fprintf("tifファイルは%0iこ",fsize)
```
ん…？
```TextOutput
tifファイルは82こ
```

そうです。3次メッシュが100こない場合もあります。


今回のファイルだと00の次が10になっていることがわかります。

![image_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/781521/2e66465e-2601-e33a-9081-23045531af77.png)


これは後で全てNaNのダミーデータを作って無理やり当てはめます。

<a name="H_0B5CDCFD"></a>
# 2次メッシュに当てはめるための格子を作成する

先程の画像のような格子をcell配列で作り、各cellに3次メッシュ情報をぶち込めばOKです。

この時にGeoTiffファイルがない部分には同じサイズのNaN行列を入れ込みます。


<a name="H_9EED0F88"></a>
## 3次メッシュ番号を抜き出す

まずは番地に当てはめるために番号を抽出しましょう。

```matlab
[~,name,~] = fileparts(file);
name = cell2mat(name);
rr  = str2double(string(name(:,end-1:end)));
disp(rr(1:4))
```

```TextOutput
     0
    10
    11
    14
```

<a name="H_BE0C7946"></a>
## 3次メッシュデータを決まった場所に格納する

次に全部Falseにした10×10の配列を作成します。


これは後でどこに3次メッシュデータがあるかを探索するように使います。

```matlab
IDX1 = false(10);
IDX2 = IDX1;
```
次に1桁目と2桁目を分離させ、行番号と列番号を抽出します。

10a + b ⇒ (a, b) にする感じです。

またMATLABなので(1, 1)からスタートする必要があるため1を各桁に加算します。

10の位と1の位を分割し、1からスタートさせるよう加算します。

```matlab
digits = [floor(rr./10), rem(rr,10)] + 1;
disp(digits(1:4,:))
```

```TextOutput
     1     1
     2     1
     2     2
     2     5
```

GeoTiffファイルを入れるため10×10のcell配列を作ります。

```matlab
Rasters = cell(10);
```

最初のcell部分はTrueにしておきます。
今回は(0, 0)スタートなので以下の部分がTrueになる感じです。

![image_3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/781521/71ad649f-f9c2-0d7f-1b91-f955e9c8b4df.png)

```matlab
IDX1(digits(:,1), digits(:,2)) = true;
```
<a name="H_EB32637C"></a>
## データのない箇所を探す

データのあるファイルの出力先をIDX2であぶりだします。

```matlab
idx1 = sub2ind(size(IDX1),digits(:,1),digits(:,2));
IDX2(idx1) = true;
```

わかりやすいように図示します。


データがない箇所を黒、データがある場所が白です。

```matlab
figure
image(IDX2)
axis xy
colormap([0, 0, 0; 1, 1, 1])
colorbar('Ticks',[0,1],'TickLabels',{'データ無','データ有'})
```

![figure_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/781521/595ddc0c-77b5-42b5-cdb8-8f47d753194f.png)



右下がデータないんですね。多分海だからでしょうね。

<a name="H_7DC8D76B"></a>
## データを読み込む

各GeoTIFFファイルを読み込み、緯度経度高度を各配列に格納します。

```matlab
Lats = zeros(fsize,2);
Lons = Lats;
RSzs = Lats;

for ii = 1:length(idx1)
    [Z,Ref] = readgeoraster(fullfile(path,file{ii}),'OutputType','double');
    Rasters(idx1(ii)) = {Z};
    Lats(ii,:) = Ref.LatitudeLimits;
    Lons(ii,:) = Ref.LongitudeLimits;
    RSzs(ii,:) = Ref.RasterSize;
end
```
<a name="H_B3C3D4F7"></a>
## データの欠損部分をNaNで埋める

排他的論理和でFalseとTrueを反転させます。

```matlab
IDX3 = xor(IDX1,IDX2);
```

Raster内で埋める必要があるところを750×1125のNaN配列で埋めます。

```matlab
idx2 = find(IDX3);
Rasters(idx2) = {nan(RSzs(1,:))};
```

残りは画像として排除していいのでNAN埋めします。

```matlab
idx3 = find(~IDX1);
Rasters(idx3) = {nan(RSzs(1,:))};
```
<a name="H_BB1314C0"></a>
## 画像データを上下逆にする

上下を反対にして一枚のdoubleデータにします。

:::note warn
MATLABでは左上が(0, 0)ですが、2次メッシュの当てはめ方は左下が(0, 0)です。

なので [flip関数](https://jp.mathworks.com/help/matlab/ref/flip.html) で順序を入れ替える必要があります。
:::

```matlab
check = flip(Rasters);
Rasters = cell2mat(flip(Rasters));
Rasters = Rasters(:,~all(isnan(Rasters),1));
Rasters = Rasters(~all(isnan(Rasters),2),:);
```
<a name="H_ED29F8EE"></a>
## モザイク結合

緯度・経度・高度情報(750×1125)をcell配列の決められた番地に格納していきます。


結合は [georefcells](https://jp.mathworks.com/help/map/ref/georefcells.html) で出来ます。

```matlab
LatsLim = [min(Lats(:,1)), max(Lats(:,2))];
LonsLim = [min(Lons(:,1)), max(Lons(:,2))];
R = georefcells(LatsLim, LonsLim,size(Rasters),"ColumnsStartFrom","north");
```

MATLAB上で問題ないか見てみましょう。

```matlab
geoshow(uint8(Rasters),R)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/781521/d880142b-ddfd-6f96-beb8-0b279a2fa2fb.png)
できてそうですね。

:::note alert
1m間隔のメッシュをgeoshowで確認しようとするとめっちゃメモリ使います。
特にライブスクリプトはクラッシュもあります。
気を付けましょう。
:::

<a name="H_0BE76304"></a>
# 出力

さて出力します。

```matlab
geotiffwrite('GeoTiff_Merge.tif',single(Rasters),R,"CoordRefSysCode",'EPSG:6668')
disp('嗚呼…モザイク結合完了…')
```

```TextOutput
嗚呼…モザイク結合完了…
```

嗚呼…モザイク結合完了…(画像略)

<a name="H_1F5E2925"></a>
# 閑話休題：EPSG

EPSGとはGISで使用される様々な要素に必要なパラメータを1つにまとめたものです。

それらの集合体同士を区別するためにEPSGコードで割り振ってるんですね。


実は地球は丸くなくて楕円体で近似を取っています。

そして日本の標高の基準は、測量法で平均海面と定められています。

この平均海面を仮想的に陸地へ延長した面をジオイドといいます。

![image_6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/781521/33635c0e-bb8c-09f0-6785-21190418852f.png)



つまり
![image_5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/781521/06ea3e0a-923c-96a8-e7af-4c7dc2ceb9c6.png)

EPSGコードから場所とかを参照して高さを求めてるんですね(テキトー解釈)。

なんだかとってもおもしろいよね。

:::note warn
※けっこうテキトーな事を言ってるので各自確認してください。
　今回はジオイドモデルではなく、数値標高モデルを使ってます。
:::

<a name="H_79E93D21"></a>
# 出力結果

さて結果を見てみましょう。


今回もフリーソフト [QGIS](https://qgis.org/ja/site/) で確認します。

![image_7.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/781521/3b884af5-9d7d-f45a-eb21-b13a9824b09b.png)


きちんとできていることがわかりました。

<a name="H_B39AF7A9"></a>
# おまけ

実はアプリも作ってます。


こちらは [waitbar](https://jp.mathworks.com/help/matlab/ref/waitbar.html) 等のメッセージボックスを入れてUIをあげてます。

![image_8.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/781521/64137c8f-6ab8-055d-9f75-914b10c53456.png)



こいつはライトニングトーク用かなあ。

<a name="H_57870575"></a>
# 参考
-  [国土地理院：基盤地図情報 数値標高モデル](https://fgd.gsi.go.jp/download/menu.php)
-  [ラスターの統合 - モザイク](https://fme-casestudy.blogspot.com/2015/07/blog-post_11.html) 
-  [複数の画像を統合（モザイク）](https://www.youtube.com/watch?v=la_aOPP3vuQ&ab_channel=ESRI%E3%82%B8%E3%83%A3%E3%83%91%E3%83%B3) 
-  [国土地理院：ジオイド](https://www.gsi.go.jp/buturisokuchi/grageo_geoid.html) 
-  [【GIS】測地系、投影法、座標系、EPSGコードについてまとめてみた](https://qiita.com/asahina820/items/a1707cba749c1fb2b192) 
-  [ジオードニキ：ジオードの中身は何が入ってる？](https://www.tiktok.com/@labnoz/video/7271493349197024520) 
-  [国土地理院 基盤地図情報 ～GIS データの基本～](https://speakerdeck.com/coolerking/guo-tu-di-li-yuan-ji-pan-di-tu-qing-bao-gis-detanoji-ben)