GeoTIFFを読み込めない

こんにちは。

GeoTIFFを読み込み、地図を表示させたいのですが、エラーになってしまいます。

どうしたらよいでしょうか？

 [A,R] = readgeoraster("Dynamic Surface.tif");
 mapshow(A,R)
 

Answer 1

[Kojiro Saito]

GeoTIFF画像の座標系が地理座標系か投影座標系で使う表示用の関数が変わります。

readgeorasterの第2出力(R)がMapCellsReferenceとかMapPostingsReferenceという投影座標系ではmapshowを使います。

例えばドキュメントの例にあるboston.tifは投影座標系の画像です。

[A,R] = readgeoraster('boston.tif');
whos R

  Name      Size            Bytes  Class                              Attributes

  R         1x1             16888  map.rasterref.MapCellsReference              

mapshow(A,R)

一方でRが GeographicCellsReferenceやGeographicPostingsReferenceの場合はgeoshowで可視化できます。

[A2,R2] = readgeoraster('n39_w106_3arc_v2.dt1','OutputType','double');
whos R2

  Name      Size            Bytes  Class                                        Attributes

  R2        1x1              1848  map.rasterref.GeographicPostingsReference              

latlim = R2.LatitudeLimits;
lonlim = R2.LongitudeLimits;
figure
usamap(latlim,lonlim)
geoshow(A2,R2,'DisplayType','surface')
demcmap(A2)

mapshowで出ているエラーから推測すると、お使いのデータの座標系(R)がMap**の投影座標系になっていないと思われます。

Geographic**の座標系になっていたら、geoshowをお使いになってはいかがでしょうか。

[Yu]

コメントくださりありがとうございます。

実はgeoshowも試したのですが、同じくエラーが出てしまい困っておりました。

出力元のソフトは他のデータも投影座標系でデータを出力しているので、このtiffデータもそうだと思うのですが…。

[Kojiro Saito]

座標系は投影座標系になっているのでmapshowで良さそうです。

エラーメッセージからすると、RGB画像なのにM×N×3のサイズになっていないのが原因のようです。

whos A

で画像部分の配列サイズを確認していただけますか？

ファイル名から推測するとこのGeoTIFF画像は地表面モデル(DSM)でしょうか。グレースケールの画像だと色情報が無いのでM×N×3の「×3」の部分が欠けているのかもしれません。

[Yu]

添付のように、カラーの海底地形モデルのgeotiffを表示させたいです。

Sizeの最後が4となっているのが原因ということでしょうか。

tiffを出力する際には何も変更ができません。

tiff以外には、xyz、KMLなどのフォーマットで出力も可能ですが、他のフォーマットでも対応可能でしょうか？

どうぞ宜しくお願い致します。

[Kojiro Saito]

3次元目のサイズが4ですが、色情報(R、G、Bの輝度値)が3列で、4列目は海底面高さでしょうか。

手元にデータが無いので確認はできませんが、表示するなら

mapshow(A(:,:,1:3), R)

で色情報が表示できると思います。

もし4列目が高さでしたら、

mapshow(A(:,:,4), R)

や

mapshow(A(:,:,4), R, 'DisplayType','surface')

で標高のテキスチャが表示できるかもしれません。

[Yu]

mapshow(A(:,:,1:3), R )

で前回添付したものと同じカラー地図が表示できました！

mapshow(A(:,:,4), R)

　→これは地形部分が真っ白に、

mapshow(A(:,:,4), R, 'DisplayType','surface')

　→こちらはエラーとなりました。

表示したところ、XY座標であったため、緯度・経度（°）に変換したいと思っています。

geotiffからXYZ座標を抽出するにはどうすればよいでしょうか？

このスレッドのままお聞きしてよいかわかりませんが、もしご存知でしたらご教示いただけますと幸いです。

[Kojiro Saito]

すると、3次元目のデータの1列から3列がRGBの情報のようですね。

mapshow(A(:,:,4), R)

で白くなってしまったのは、4列目のデータの値が高い値に偏っているためと思われます。

 histogram(A(:,:,4))

で散らばりが確認できますが、横軸の右側に偏っていませんでしょうか？

データを正規化して色の散らばりを変えてみると見やすくなるかもしれません。

 mapshow(normalize(double(A(:,:,4)), 'range')*256, R)

また、mapshowはXY座標なので、lat、lonで表示するにはprojinvで座標系を変換してgeoshowで表す必要があります。

参考：ドキュメント

[A,R] = readgeoraster("Dynamic Surface.tif");
p = R.ProjectedCRS;
[x,y] = worldGrid(R);
[lat,lon] = projinv(p,x,y);
geoshow(lat,lon, A)

[Yu]

返信が遅くなり申し訳ありません。

【１】

histogram(A(:,:,4))
で散らばりが確認できますが、横軸の右側に偏っていませんでしょうか？
データを正規化して色の散らばりを変えてみると見やすくなるかもしれません。

→確認してみると、左右に偏っていました。正規化したものも、地形部分が以前と変わらず真っ白でした。

【２】lat lonへの変換はできたのですが、geoshowでエラーが発生してしまいました。geoshowのドキュメントも読んだのですが、解決方法が分かりませんでした…。

何度も申し訳ございませんが、ご教示いただけますと幸いです。

[Kojiro Saito]

@Yu　さん

使用しているGeoTIFFデータをシェアいただくのは可能でしょうか？

ドキュメントの例にあるGeoTIFFのサンプルと異なっているようなので、実際にデータを見させてもらったほうがスムーズそうです。

[Yu]

ありがとうございます。添付にサンプルデータを添付を致しましたので、お手数ですがご確認のほど宜しくお願いいたします。

[Kojiro Saito]

@Yu

添付ありがとうございます。

A(:, :, 4)は0か255の2値になっていて、ピクセルの輝度値が画像の領域内なら255、領域外なら0になっているだけでしたね。

A(:, :, 1:3)の色情報を緯度、経度にマッピングするには不要の情報でした。

[A, R] = readgeoraster("Sample.tif");
p = R.ProjectedCRS;
[x,y] = worldGrid(R);
[lat,lon] = projinv(p,x,y);
geoshow(lat,lon, A(:,:,1:3))

geoshowの中に緯度経度座標に投影したAの色情報(1:3)だけ入れればできました。

[Yu]

なるほど、輝度値というのがあるのですね。勉強になりました。

こちらでも、元のデータで無事緯度経度に変換できました。

何度も丁寧にお答えいただき大変助かりました！ありがとうございました。