RGB画像の作成方法

RGB画像

RGB画像を作る時は三つのバンドを利用する。
それぞれのバンドを青、緑、赤に着色して合成カラー画像を作成する。
ここではバイナリデータ形式の大気上端反射率を使用している。
Sentinel2Aにおいて私が使ったバンドは下記の通りである。

VN02を青の波長帯
VN03を緑の波長帯
VN04を赤の波長帯
VN08を近赤外の波長帯
VN11,SW03なども利用した。

#ture color 本来の色が見える
rgb_3d = np.zeros((nxh,nyh,3), dtype=np.uint8)
rgb_3d[:, :, 0]=brefv04[:, :] #Red
rgb_3d[:, :, 1]=brefv03[:, :] #Green
rgb_3d[:, :, 2]=brefv02[:, :] #Blue
img = Image.fromarray(rgb_3d, 'RGB')
# natural color 植物の特徴的に示すことができる
rgb_3d = np.zeros((nxh,nyh,3), dtype=np.uint8)
rgb_3d[:, :, 0]=brefv04[:, :] #Red
rgb_3d[:, :, 1]=brefv08[:, :] #Green
rgb_3d[:, :, 2]=brefv03[:, :] #Blue
img = Image.fromarray(rgb_3d, 'RGB')
#false color 
rgb_3d = np.zeros((nxh,nyh,3), dtype=np.uint8)
rgb_3d[:, :, 0]=brefv08[:, :] #Red
rgb_3d[:, :, 1]=brefv04[:, :] #Green
rgb_3d[:, :, 2]=brefv03[:, :] #Blue
img = Image.fromarray(rgb_3d, 'RGB')
#よく利用したパターン(色が鮮やかに出る)
rgb_3d = np.zeros((nxh,nyh,3), dtype=np.uint8)
rgb_3d[:, :, 0]=brefsw03[:, :] #Red
rgb_3d[:, :, 1]=brefv11[:, :] #Green
rgb_3d[:, :, 2]=brefv08[:, :] #Blue
img = Image.fromarray(rgb_3d, 'RGB')

plt.imshow(img, origin='upper')
plt.show()

今日はここまで!

衛星データの内容を知るためには

データの参照方法

衛星データがどのような形式、大きさであるか把握していることは
画像加工する上で重要である。
特にエラーが起きた時はこれらを確認しておくと、エラーの原因がわかったりする。

img = file_name_v02
px = img[100,100]
print(px) #指定した画素中の値
print(type(img)) #画像タイプ
print(img.dtype) #データタイプ
print(img.shape) #データの大きさ:例えば(100,100)など

||

衛星画像の切り取り方法

衛星画像のデータが大きいと画像処理に何かと負荷がかかったり、時間がかかる。

そこで、観測地域が限られている時は画像の切り取りを行う。



衛星画像の切り取り

方法は主に二つある。
①緯度経度を使用する場合

ul_lon=135 #画像左上経度 
ul_lat= 40 #画像左上緯度
lr_lon=145 #画像右下経度
lr_lat= 30 #画像右下緯度

img_extent = [ul_lon, lr_lon, lr_lat, ul_lat]
plt.imshow(img, extent=img_extent)  

②緯度経度を使用しない場合(画素の数で切り取り)
下の画像で、青部分は元画像、橙部分は切り取り後の画像
切り取りたい画像部分の左上端を起点に幅と高さを決める
f:id:pmaryretro:20210203185834p:plain

left = 8235 #左上端のx座標
top = 8235 #左上端のy座標
width = 2745 #横幅
height = 2745 #高さ
img02 = img01[top:top+height,left:left+width]
plt.imshow(img02)  

今日はここまで!

SNAPの使い方(Sentinel加工)

「SNAP」とは

ESA (European Space Agency, 欧州宇宙機構)が所有している衛星の他、さまざまな衛星データを扱えるオープンソフトの衛星画像処理ソフトである。

 

衛星データの大まかな内容を把握するのに楽!

より詳細な加工を行うのはpythonの方がやりやすい。

Sentinelのデータだと大気上端反射率(0-1)で表示・情報をみることができる。

 

Sentinel1、Sentinel2、Sentinel3や日本のALOSでも利用できる。

 

ここには以前使い方をまとめたスライドを載せておく

一枚目は基本事項の説明。

f:id:pmaryretro:20210203180517p:plain

 

二枚目はリサンプリングの説明である。

 

f:id:pmaryretro:20210203180524p:plain

 

三枚目はさまざまな指数の説明である。

 

f:id:pmaryretro:20210203180531p:plain

 

四枚目はバンド計算の説明である。

 

f:id:pmaryretro:20210203180519p:plain

 

今日はここまで!

 

Sentinel2Aのデータ構造

Sentinel2Aの処理レベル

Sentinel2Aは「1Cレベル」という処理レベルになっている

簡単にまとめると「ジオメトリック、ラジオメトリック、幾何補正などを行なったのちに

大気上端反射率を計算している」状態。

そのため、ほとんど補正を自分の手で行わずに次の段階に進むことができる。

読み込むときにはバイナリデータ(0-255)になっている。

 

Sentinel2Aのデータダウンロード

私は普段earthexplorerで衛星データのダウンロードを行なっている。

ダウンロード方法は多くのネット記事があるので省略。

アドレスは下記の通りである。

https://earthexplorer.usgs.gov/

ダウンロードしたファイルはこんな感じになっている。

f:id:pmaryretro:20210203173259p:plain

データファイルの詳細

「GRANULE」を開いたあと、さらに「IMG_DATA」を開くと画像データが得られる。

f:id:pmaryretro:20210203173631p:plain

GRANULEの詳細

 

きょうはここまで!

 

バイナリデータ(0-255)から大気上端反射率(0-1)に直すために係数「slop」,「offset」

を利用したいのだが、これがどこにかかれているのかがわからないので分かり次第再記

Sentinel2の衛星画像(.jp2)の読み込むためには

 

 

Sentinel2Aの画像はjpeg2000の形式になっている。

そのため、それに対応した読込を行わなければいけないのだが、、

これが面倒くさい

ネットにはpng,jpeg形式の読込方法はよく転がってるのだが、

jpeg2000になるとその数がぐんと減少する

 

今回はOpenCVを用いる方法を残していく

import cv2

file_name = cv2.imread('filename.jp2',cv2.IMREAD_GRAYSCALE)

結構簡単で適応力があるので、OpenCVはおすすめ

因みに

import cv2

file_name = cv2.imread('filename.jp2',cv2.IMREAD_COLOR)

file_name = cv2.imread('filename.jp2',cv2.IMREAD_GRAYSCALE)

file_name = cv2.imread('filename.jp2',cv2.IMREAD_UNCHANGED)

file_name = cv2.imread('filename.jp2',1)

file_name = cv2.imread('filename.jp2',0)

file_name = cv2.imread('filename.jp2',-1)

IMREAD_COLOR (1)             カラーで読み込む

IMREAD_GRAYSCALE (0)    グレースケールで読み込む

IMREAD_UNCHANGED  (-1) アルファスケールで読み込む(透過率の配列を含む)