GEE技術検証 | しきさいポータル

GEE利用事例

北海道の海氷を見てみよう：RGB合成とGEE上でのGIF作成について

海氷を衛星画像から見てみよう。

しずく（GCOM-W）が観測した、2023年の海氷移動の様子がGCOMインスタグラムに投稿されています。
（GCOMインスタURL：https://www.instagram.com/jaxajp/reel/C3pJ4YPi2l1/）

搭載センサの種類や観測メカニズムの違いがありますが、しきさい（GCOM-C）も同様に海氷を観測することが出来ます。

参考：
「地球がみえる」掲載記事：「打ち上げから5年目を迎える『しきさい』~雪氷環境編~」
https://earth.jaxa.jp/ja/earthview/2022/09/05/7209/index.html

海氷の様子を、GEEにアップロードされているGCOM-CのLevel-1B（L1B）放射輝度データからも確認してみましょう。

GEEのPlatformのコンソール画面を立ち上げたら、まずはL1Bプロダクトの前処理を行うとRGB合成をする関数を定義します。
詳しくは、チュートリアル「ImageCollectionのマップ表示(RGB①)」ページの「④Slope、Offsetの計算」・「⑤太陽天頂角の補正」をご覧ください。
ただし、今回は大気上端反射率を計算するために、プロパティは”Slope_reflectance”・”Offset_reflectance”を用います。
また、この関数では人間の目で見る様子に近いRGB合成ができるよう、VN08を赤、VN05を緑、VN03を青に割り当てる処理を行います。


// L1Bの前処理・可視光域を使用したRGB合成を行う関数を定義する
function preprocessing_RGB1(img){
  var vn03        = img.select("Lt_VN03");
  var vn05        = img.select("Lt_VN05");
  var vn08        = img.select("Lt_VN08");
  var sza         = img.select("Solar_zenith_angle");
  
  var vn03_slope  = ee.Number.parse(img.get('VNR_VN03_Slope_reflectance')).float();
  var vn05_slope  = ee.Number.parse(img.get('VNR_VN05_Slope_reflectance')).float();
  var vn08_slope  = ee.Number.parse(img.get('VNR_VN08_Slope_reflectance')).float();
  var sza_slope   = ee.Number.parse(img.get('SZA_Slope')).float();

  var vn03_offset = ee.Number.parse(img.get('VNR_VN03_Offset_reflectance')).float();
  var vn05_offset = ee.Number.parse(img.get('VNR_VN05_Offset_reflectance')).float();
  var vn08_offset = ee.Number.parse(img.get('VNR_VN08_Offset_reflectance')).float();
  var sza_offset  = ee.Number.parse(img.get('SZA_Offset')).float();
  
  var vn03_tmp    = vn03.multiply(vn03_slope).add(vn03_offset);
  var vn05_tmp    = vn05.multiply(vn05_slope).add(vn05_offset);
  var vn08_tmp    = vn08.multiply(vn08_slope).add(vn08_offset);
  var sza_cor     = sza.multiply(sza_slope).add(sza_offset);
    
  var sza_rad     = sza_cor.divide(180).multiply(Math.PI);
  
  var vn03_cor    = vn03_tmp.divide(sza_rad.cos());
  var vn05_cor    = vn05_tmp.divide(sza_rad.cos());
  var vn08_cor    = vn08_tmp.divide(sza_rad.cos());
  
  var rgb         = ee.Image.rgb(vn08_cor, vn05_cor, vn03_cor);
  
  return rgb
}

その後、GEE上のL1Bプロダクトをロードし、先ほど定義したpreprocessing_RGB1関数をImageCollection内の各Imageに適用していきます。
2023/02/07を、今回の表示の対象とし、ロード時点でデータを絞ります。


// GCOM-C L1Bをロードする。
var l1b_20230207 = ee.ImageCollection('projects/ee-gcomsgli/assets/L1BDQ')
                     .filterDate('2023-02-07T00:00:00', '2023-02-08T00:00:00');

// preprocessing_RGB1をImageに適用する
var rgb1_20230207 = l1b_20230207.map(preprocessing_RGB1);

その後視覚化パラメータを定義し、マップ表示します。
今回は北海道の北方沖が表示範囲に含まれるような設定をします。


// 視覚化パラメータを定義する
var imageVisParam_rgb1 = {
  min: [  0,   0,   0],
  max: [0.8, 0.8, 0.8]
}

// 表示位置の中心を設定して、マップ表示する。
Map.setCenter(146.0, 45.5, 7);
Map.addLayer(rgb1_20230207, imageVisParam_rgb1, 'RGB1_20230207');

北海道北方沖のRGBを表示することができました。
しかし、この画像をみると雲と海氷の両方とも白色で表示されています。

雲と海氷の区別をつきやすくするために、赤・緑・青に割り当てるバンドを変えていきます。オホーツク海氷モニタで公開されているGCOM-C画像の１つに、VN05を赤、VN11を緑、SW03を青に割り当てているものがあるため、今回はこのバンドの組み合わせを用いたRGB合成を行います。
また、バンドの組み合わせを変えたため新たな視覚化パラメータを定義し、表示させます。


// L1Bの前処理・雪氷を目立つようにしたRGB合成を行う関数を定義する
function preprocessing_ohk_RGB(img){
  var vn05        = img.select("Lt_VN05");
  var vn11        = img.select("Lt_VN11");
  var sw03        = img.select("Lt_SW03");
  var sza         = img.select("Solar_zenith_angle");
  
  var vn05_slope  = ee.Number.parse(img.get('VNR_VN03_Slope_reflectance')).float();
  var vn11_slope  = ee.Number.parse(img.get('VNR_VN11_Slope_reflectance')).float();
  var sw03_slope  = ee.Number.parse(img.get('IRS_SW03_Slope_reflectance')).float();
  var sza_slope   = ee.Number.parse(img.get('SZA_Slope')).float();

  var vn05_offset = ee.Number.parse(img.get('VNR_VN03_Offset_reflectance')).float();
  var vn11_offset = ee.Number.parse(img.get('VNR_VN11_Offset_reflectance')).float();
  var sw03_offset = ee.Number.parse(img.get('IRS_SW03_Offset_reflectance')).float();
  var sza_offset  = ee.Number.parse(img.get('SZA_Offset')).float();
  
  var vn05_tmp    = vn05.multiply(vn05_slope).add(vn05_offset);
  var vn11_tmp    = vn11.multiply(vn11_slope).add(vn11_offset);
  var sw03_tmp    = sw03.multiply(sw03_slope).add(sw03_offset);
  var sza_cor     = sza.multiply(sza_slope).add(sza_offset);
    
  var sza_rad     = sza_cor.divide(180).multiply(Math.PI);
  
  var vn05_cor    = vn05_tmp.divide(sza_rad.cos());
  var vn11_cor    = vn11_tmp.divide(sza_rad.cos());
  var sw03_cor    = sw03_tmp.divide(sza_rad.cos());
  
  var rgb         = ee.Image.rgb(vn05_cor, vn11_cor, sw03_cor);
  
  return rgb
}

// 海氷RGBの視覚化パラメータを定義する
var imageVisParam_rgb_ohk = {
  min: [ 0,   0,  0],
  max: [1.2, 0.8, 0.5]
}

// preprocessing_ohk_RGBをImageに適用する
var rgb_ohk_20230207 = l1b_20230207.map(preprocessing_ohk_RGB);

// 表示位置の中心を設定して、マップ表示する。
Map.setCenter(146.0, 45.5, 7);
Map.addLayer(rgb_ohk_20230207, imageVisParam_rgb_ohk, 'RGB_ohk_20230207');

今度のRGB合成では、海氷が濃い黄色、雲が明るい黄色から水色、と区別できる色付けになりました。
このように多波長を観測する衛星データを扱う際は、対象物の反射特性に合わせたバンドを組み合わせてRGB合成を行うことで、対象物をより目立たせるような画像を作ることができます。

GCOM-Cに搭載されているSGLIセンサは、今回用いたバンド以外にもさまざまな波長帯を観測するバンドがあります。
皆さんも、いろんなバンドの組み合わせによるRGB合成を試してみてください。

「SGLIの観測波長と主なプロダクト」
https://suzaku.eorc.jaxa.jp/GCOM_C/instruments/product_j.html

またチュートリアル「ImageCollectionのマップ表示(RGB②)」では、今回とは違うバンドの組み合わせでRGＢ合成をしていますので、こちらもぜひご一読ください。

チュートリアル「ImageCollectionのマップ表示(RGB②)」
https://shikisai.jaxa.jp/GEE/Tutorial/tutorial_003.html

海氷の移動の様子をとらえよう

前項では、１日分のデータについてRGB合成を表示しました。
本項目では、海氷の移動や時間経過による変化の様子を確認するために、複数日分のデータについてRGB合成を行い、連続して表示できるようGIF画像に変換していきます。

まず、今回GIF画像を作りたい領域を定義し、その範囲を含むようなGCOM-C L1Bデータをロードします。
今回は雲の少ないデータのみでGIF画像を作りたいため、あらかじめ目視で2023年1月から3月前半で雲の少ない日を選び、それぞれ個別のImageCollectionとしてロードします。


// 今回海氷の移動を見たい領域（GIF画像の範囲; 北海道の北方沖）を定義する
var geometry = ee.Geometry.Rectangle([142, 47.72, 146, 43.45]);

// 目視で雲の少ない日を選び、それぞれの日について
// GCOM-C L1Bデータをロードする。
// 2023/01
var l1b_20230127 = ee.ImageCollection('projects/ee-gcomsgli/assets/L1BDQ')
                     .filterBounds(geometry)
                     .filterDate('2023-01-27T00:00:00', '2023-01-28T00:00:00');
var l1b_20230128 = ee.ImageCollection('projects/ee-gcomsgli/assets/L1BDQ')
                     .filterBounds(geometry)
                     .filterDate('2023-01-28T00:00:00', '2023-01-29T00:00:00');
var l1b_20230131 = ee.ImageCollection('projects/ee-gcomsgli/assets/L1BDQ')
                     .filterBounds(geometry)
                     .filterDate('2023-01-31T00:00:00', '2023-02-01T00:00:00');

// 2023/02       
var l1b_20230201 = ee.ImageCollection('projects/ee-gcomsgli/assets/L1BDQ')
                     .filterBounds(geometry)
                     .filterDate('2023-02-01T00:00:00', '2023-02-02T00:00:00');
var l1b_20230203 = ee.ImageCollection('projects/ee-gcomsgli/assets/L1BDQ')
                     .filterBounds(geometry)
                     .filterDate('2023-02-03T00:00:00', '2023-02-04T00:00:00');
var l1b_20230207 = ee.ImageCollection('projects/ee-gcomsgli/assets/L1BDQ')
                     .filterBounds(geometry)
                     .filterDate('2023-02-07T00:00:00', '2023-02-08T00:00:00');
var l1b_20230208 = ee.ImageCollection('projects/ee-gcomsgli/assets/L1BDQ')
                     .filterBounds(geometry)
                     .filterDate('2023-02-08T00:00:00', '2023-02-09T00:00:00');
var l1b_20230211 = ee.ImageCollection('projects/ee-gcomsgli/assets/L1BDQ')
                     .filterBounds(geometry)
                     .filterDate('2023-02-11T00:00:00', '2023-02-12T00:00:00');
var l1b_20230222 = ee.ImageCollection('projects/ee-gcomsgli/assets/L1BDQ')
                     .filterBounds(geometry)
                     .filterDate('2023-02-22T00:00:00', '2023-02-23T00:00:00');

                   
// 2023/03
var l1b_20230303 = ee.ImageCollection('projects/ee-gcomsgli/assets/L1BDQ')
                     .filterBounds(geometry)
                     .filterDate('2023-03-03T00:00:00', '2023-03-04T00:00:00');
var l1b_20230306 = ee.ImageCollection('projects/ee-gcomsgli/assets/L1BDQ')
                     .filterBounds(geometry)
                     .filterDate('2023-03-06T00:00:00', '2023-03-07T00:00:00');

その後、前処理・RGB合成の関数を定義します。内容は前項目ほぼ同様ですが、後続の処理のため、この関数の時点で視覚化パラメータを適用します。この関数を各日のImageCollectionに適用すると同時に、mosaic関数も適用し、1日分のRGBデータのImageに変換します。このImageの変換も後続の処理のために行います。


// 視覚化パラメータを定義する
var imageVisParam_rgb_ohk = {
  min: [ 0,    0,  0],
  max: [1.2, 0.8, 0.5]
}

// L1Bの前処理・雪氷を目立つようにしたRGB合成を行う関数を定義する
function preprocessing_RGB_ohk(img){
  var vn05        = img.select("Lt_VN05");
  var vn11        = img.select("Lt_VN11");
  var sw03        = img.select("Lt_SW03");
  var sza         = img.select("Solar_zenith_angle");
  
  var vn05_slope  = ee.Number.parse(img.get('VNR_VN03_Slope_reflectance')).float();
  var vn11_slope  = ee.Number.parse(img.get('VNR_VN11_Slope_reflectance')).float();
  var sw03_slope  = ee.Number.parse(img.get('IRS_SW03_Slope_reflectance')).float();
  var sza_slope   = ee.Number.parse(img.get('SZA_Slope')).float();

  var vn05_offset = ee.Number.parse(img.get('VNR_VN03_Offset_reflectance')).float();
  var vn11_offset = ee.Number.parse(img.get('VNR_VN11_Offset_reflectance')).float();
  var sw03_offset = ee.Number.parse(img.get('IRS_SW03_Offset_reflectance')).float();
  var sza_offset  = ee.Number.parse(img.get('SZA_Offset')).float();
  
  var vn05_tmp    = vn05.multiply(vn05_slope).add(vn05_offset);
  var vn11_tmp    = vn11.multiply(vn11_slope).add(vn11_offset);
  var sw03_tmp    = sw03.multiply(sw03_slope).add(sw03_offset);
  var sza_cor     = sza.multiply(sza_slope).add(sza_offset);
    
  var sza_rad     = sza_cor.divide(180).multiply(Math.PI);
  
  var vn05_cor    = vn05_tmp.divide(sza_rad.cos());
  var vn11_cor    = vn11_tmp.divide(sza_rad.cos());
  var sw03_cor    = sw03_tmp.divide(sza_rad.cos());
  
  //前項目との変更点
  var rgb         = ee.Image.rgb(vn05_cor, vn11_cor, sw03_cor).visualize(imageVisParam_rgb_ohk);
  
  return rgb
}

var rgb_ohk_20230127 = l1b_20230127.map(preprocessing_RGB_ohk)
                       .mosaic();

※上記のコードの最後に、2023/01/27１日分のみpreprocessing_RGB_ohk関数とモザイク処理のコード表記をしていますが、実際には今回ロードしたL1BのImageCollectionすべてに同じ処理を行います。

今回ロードしたL1BのImageCollectionすべてをImageに変換出来たら、各日のRGBのImageを1つのImageCollectionにまとめます。


// 各日のImageをまとめて1つのImageCollectionにする.
var rgb_ohk_202301_202303 = ee.ImageCollection.fromImages(
                       [rgb_ohk_20230127,
                        rgb_ohk_20230128,
                        rgb_ohk_20230131,
                        rgb_ohk_20230201,
                        rgb_ohk_20230203,
                        rgb_ohk_20230207,
                        rgb_ohk_20230208,
                        rgb_ohk_20230211,
                        rgb_ohk_20230222,
                        rgb_ohk_20230303,
                        rgb_ohk_20230306]);

そしてGIF化パラメータを定義し、RGBのImageをまとめたImageCollectionにgetVideoThumbURL関数を適用します。


// GIF化パラメータを定義する
var gifParams = {
  'region': geometry,
  'dimensions': 1200,
  'crs': 'EPSG:4326',
  'framesPerSecond': 1.8,
  'format': 'gif'
};

// GIF化する。ConsoleにGIF画像を表示・ダウンロードできるURLが表示される。
print(rgb_ohk_202301_202303.getVideoThumbURL(gifParams))

今回のコードを実行し、Consoleに表示されるURLから表示できるGIF画像は以下です。

日の間隔がまちまちのため、滑らかな海氷移動の様子は見られませんでしたが、海氷が南下し、その後東へ移動する様子は確認することができました。