目次
はじめに
今回やること
フォルダを読み込む
2次メッシュに当てはめるための格子を作成する
3次メッシュ番号を抜き出す
3次メッシュデータを決まった場所に格納する
データのない箇所を探す
データを読み込む
データの欠損部分をNaNで埋める
画像データを上下逆にする
モザイク結合
出力
閑話休題：EPSG
出力結果
おまけ
参考
はじめに
本件は日本の標高データをMATLABで活用 - GeoTiff変換 -の続きです。
先回はDEMデータをGeoTiffファイルにしました。
しかしきれいに繋がらなかったり、そもそも3次メッシュが一部欠損してたり…
これらを(0, 0) ~ (9,9)までの区画に当てはめて2次メッシュをきれいに結合します。
モザイク結合についてはラスターの統合 - モザイクを見てください。
今回やること
3次メッシュを100枚集めて結合し2次メッシュにします。
イメージとしては以下のように3次メッシュ5841-46-43を(4, 3)に当てはめていく感じです。
フォルダを読み込む
先回作った5841-46-XXのメッシュでXXは(0, 0) ~ (9, 9)にそれぞれ当てはまる合計100このtifファイルがあるはずです。
まずは先回作ったtmpファイル内の中を見ていきましょう。
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枚のファイルを並べてモザイク結合を…
fprintf("tifファイルは%0iこ",fsize)
ん…？
tifファイルは82こ
そうです。3次メッシュが100こない場合もあります。
今回のファイルだと00の次が10になっていることがわかります。
これは後で全てNaNのダミーデータを作って無理やり当てはめます。
2次メッシュに当てはめるための格子を作成する
先程の画像のような格子をcell配列で作り、各cellに3次メッシュ情報をぶち込めばOKです。
この時にGeoTiffファイルがない部分には同じサイズのNaN行列を入れ込みます。
3次メッシュ番号を抜き出す
まずは番地に当てはめるために番号を抽出しましょう。
[~,name,~] = fileparts(file);
name = cell2mat(name);
rr = str2double(string(name(:,end-1:end)));
disp(rr(1:4))
0
10
11
14
3次メッシュデータを決まった場所に格納する
次に全部Falseにした10×10の配列を作成します。
これは後でどこに3次メッシュデータがあるかを探索するように使います。
IDX1 = false(10);
IDX2 = IDX1;
次に1桁目と2桁目を分離させ、行番号と列番号を抽出します。
10a + b ⇒ (a, b) にする感じです。
またMATLABなので(1, 1)からスタートする必要があるため1を各桁に加算します。
10の位と1の位を分割し、1からスタートさせるよう加算します。
digits = [floor(rr./10), rem(rr,10)] + 1;
disp(digits(1:4,:))
1 1
2 1
2 2
2 5
GeoTiffファイルを入れるため10×10のcell配列を作ります。
Rasters = cell(10);
最初のcell部分はTrueにしておきます。
今回は(0, 0)スタートなので以下の部分がTrueになる感じです。
IDX1(digits(:,1), digits(:,2)) = true;
データのない箇所を探す
データのあるファイルの出力先をIDX2であぶりだします。
idx1 = sub2ind(size(IDX1),digits(:,1),digits(:,2));
IDX2(idx1) = true;
わかりやすいように図示します。
データがない箇所を黒、データがある場所が白です。
figure
image(IDX2)
axis xy
colormap([0, 0, 0; 1, 1, 1])
colorbar('Ticks',[0,1],'TickLabels',{'データ無','データ有'})
右下がデータないんですね。多分海だからでしょうね。
データを読み込む
各GeoTIFFファイルを読み込み、緯度経度高度を各配列に格納します。
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
データの欠損部分をNaNで埋める
排他的論理和でFalseとTrueを反転させます。
IDX3 = xor(IDX1,IDX2);
Raster内で埋める必要があるところを750×1125のNaN配列で埋めます。
idx2 = find(IDX3);
Rasters(idx2) = {nan(RSzs(1,:))};
残りは画像として排除していいのでNAN埋めします。
idx3 = find(~IDX1);
Rasters(idx3) = {nan(RSzs(1,:))};
画像データを上下逆にする
上下を反対にして一枚のdoubleデータにします。
MATLABでは左上が(0, 0)ですが、2次メッシュの当てはめ方は左下が(0, 0)です。
なので flip関数 で順序を入れ替える必要があります。
check = flip(Rasters);
Rasters = cell2mat(flip(Rasters));
Rasters = Rasters(:,~all(isnan(Rasters),1));
Rasters = Rasters(~all(isnan(Rasters),2),:);
モザイク結合
緯度・経度・高度情報(750×1125)をcell配列の決められた番地に格納していきます。
結合は georefcells で出来ます。
LatsLim = [min(Lats(:,1)), max(Lats(:,2))];
LonsLim = [min(Lons(:,1)), max(Lons(:,2))];
R = georefcells(LatsLim, LonsLim,size(Rasters),"ColumnsStartFrom","north");
MATLAB上で問題ないか見てみましょう。
geoshow(uint8(Rasters),R)
できてそうですね。
1m間隔のメッシュをgeoshowで確認しようとするとめっちゃメモリ使います。
特にライブスクリプトはクラッシュもあります。
気を付けましょう。
出力
さて出力します。
geotiffwrite('GeoTiff_Merge.tif',single(Rasters),R,"CoordRefSysCode",'EPSG:6668')
disp('嗚呼…モザイク結合完了…')
嗚呼…モザイク結合完了…
嗚呼…モザイク結合完了…(画像略)
閑話休題：EPSG
EPSGとはGISで使用される様々な要素に必要なパラメータを1つにまとめたものです。
それらの集合体同士を区別するためにEPSGコードで割り振ってるんですね。
実は地球は丸くなくて楕円体で近似を取っています。
そして日本の標高の基準は、測量法で平均海面と定められています。
この平均海面を仮想的に陸地へ延長した面をジオイドといいます。
つまり
EPSGコードから場所とかを参照して高さを求めてるんですね(テキトー解釈)。
なんだかとってもおもしろいよね。
※けっこうテキトーな事を言ってるので各自確認してください。
今回はジオイドモデルではなく、数値標高モデルを使ってます。
出力結果
さて結果を見てみましょう。
今回もフリーソフト QGIS で確認します。
きちんとできていることがわかりました。
おまけ
実はアプリも作ってます。
こちらは waitbar 等のメッセージボックスを入れてUIをあげてます。
こいつはライトニングトーク用かなあ。