[GEEネタ]Google Earth Engineのサンプルコード写経となんちゃって解説~Image編~

はじめに

Google Earth Engineというものがあります。https://earthengine.google.com

Google Earth Engineは、数ペタバイト規模の衛星画像と地理空間データセットを惑星規模の解析機能と組み合わせて、科学者、研究者、開発者が変化を検出し、傾向をマッピングし、地球表面の違いを定量化できるようにします。 webサイトより引用

非常に興味深いので、使えるようになりたいと思っています。そこでGoogle Earth Engine（以降GEE）を使えるようになるため、Google が用意しているサンプルコードを写経していきたいと思います。ただし、ただの写経だとあんまり味気ないため、下記3つを自分に課したいと思います。
■Phase１：まずは普通に写経（コピペは厳禁）
■Phase２：コードの解説
■Phase３：写経したコードを元に少し改良する（余裕があったら）

QiitaはGoogle Earth Engineのシンタックスハイライトに対応していないので、少し見にくいです。

サンプルコードについて

Google が用意しているサンプルコードは

• Image
• Image Collection
• Feature Collection
• Charts
• Arrays
• Primitive
• Cloud Masking
• Code Editor
• Demos
• User Interface
  の10コースあります。

今回は、Imageのサンプルコードを読み解いて行きたいと思います。

Image

Imageは、一番最初のコースです。画像そのものの呼び出し方や扱い方が記載サれています。

##From Name
基本的な3つを学びます。

1. ee.Imageによるデータの呼び出し
2. Map.setCenterによる描写位置の緯度経度指定
3. Map.addLayerによるディスプレイ表示

###Phase１
コードは下記のとおり

var Image = ee.Image('CGIAR/SRTM90_V4');
Map.setCenter(-110,40,5);
Map.addLayer(Image,{min:0,max:3000},'SRTM')

###Phase２
1行目；
SRTMのデータを読み込んでいる。
SRTMについてはここを参照してください。（http://www.kashmir3d.com/srtm/）
2行目；
実行時、どの座標を中心にして描写するか決めている
引数は以下の通り；
Map.setCenter([軽度],[緯度],[拡大レベル（0-24）]);
3行目；
地図上にレイヤを追加する。
レイヤを追加するのはImage（つまりSRTMデータ）であり、表示する際のDN値を0～3000としている。

###Phase３
Map.setCenterにて指定されている場所がアメリカのため、日本にしたい。
そのため、Map.setCenter(140,35,5);とする。
また、minやmaxといった表示レベルは別に変数に格納しておき、可読性を高める。
以下が少し改良したコード

var Image = ee.Image('CGIAR/SRTM90_V4');
var lat = 140;
var lon = 35;
var z = 8;
var Params = {min:0,max:3000};
Map.setCenter(lat,lon,z);
Map.addLayer(Image,Params,'SRTM')

Normalized Difference

衛星データでは、よくRとNIRのデータを使ってNDVIを（正規化植生指標）を計算することがある。ここではNDVIの計算方法を学ぶ。
またpaletteの指定方法を学ぶ。
###Phase1
コードは下記の通り

var image = ee.Image('MODIS/006/MOD09GA/2012_03_09');
var ndvi = image.normalizedDifference(['sur_refl_b02','sur_refl_b01']);
var palette = ['FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718',
               '74A901', '66A000', '529400', '3E8601', '207401', '056201',
               '004C00', '023B01', '012E01', '011D01', '011301'];
               
Map.setCenter(-94.84497,39.01918,8);

Map.addLayer(image.select(['sur_refl_b01','sur_refl_b04','sur_refl_b03']),
            {gain:[0.1,0.1,0.1]},'MODIS bands 1/4/3')
Map.addLayer(ndvi,{min:0,max:1,palette:palette},'NDVI')

###Phase2
1行目；
MODISプロダクトのMOD09GAを呼び出している。MOD09GAについてはここを参照（https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/mod09ga_v006）
2行目；
NDVI（正規化植生指標）を計算している。
normalizedDifference([NIRを指定,REDを指定]);
とすることで、NDVIの計算が可能。

3行目；
値別の色を指定している。

4行目；
実行時、どの座標を中心にして描写するか決めている

5行目；
Selectでどのバンドのデータを使うか指定し、さらにgainにてゲインの調整をしている。
また、レイヤー名を「MODIS bands 1/4/3」としている。

6行目；
NDVIをaddLayerにて表示している。

MODIS bands 1/4/3を下記に示す。

NDVIを下記に示す。

###Phase3
少し読みやすくした

var image = ee.Image('MODIS/006/MOD09GA/2012_03_09');

var NIR = image.select('sur_refl_b02')
var RED = image.select('sur_refl_b01')
var ndvi = image.normalizedDifference(['sur_refl_b02','sur_refl_b01'])
var lat = 140
var lon = 35
var z = 9
var bandParams = {gain:[0.1,0.1,0.1]}
var ndviParams= {min:0 , max:1,palette:palette}
var palette = ['FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718',
               '74A901', '66A000', '529400', '3E8601', '207401', '056201',
               '004C00', '023B01', '012E01', '011D01', '011301'];
Map.setCenter(lat,lon,z)
Map.addLayer(image.select(['sur_refl_b01','sur_refl_b04','sur_refl_b03']),bandParams,'MODIS band 1/4/3')
Map.addLayer(ndvi,ndviParams,'NDVI')

image.normalizedDifferenceの引数をNIRとREDとしたらエラーができた。
変数に格納して、それをNDVIの計算に使うことは多分できないっぽい。

##Expression
ここでは、EVIと呼ばれる植生指標を計算しています。NDVIをは計算する際に地面が含まれていると値が著しく下がってしまうことが知られている。それを改善したのが　EVIです。
また、EVIの計算はNDVIをよりも少し複雑です。ここではexpressionを使って計算しています。
###Phase1
コードは下記の通り

var image = ee.Image('MODIS/006/MOD09GA/2012_03_09').multiply(0.0001);
var evi = image.expression(
  '2.5 * (nir -red) / (nir + 6 *red -7.5 * blue + 1)',
  {
    red:image.select('sur_refl_b01'),
    nir:image.select('sur_refl_b02'),
    blue:image.select('sur_refl_b03')
  });

Map.setCenter(-94.84497, 39.01918, 8);
Map.addLayer(image.select(['sur_refl_b01','sur_refl_b04','sur_refl_b03']),
      {min:0,max:1},'MODIS bands 1/4/3')
Map.addLayer(evi,{min:0,max:1},'EVI')

###Phase2
今までのコードとの変更点は、

1. 1行目の読み込み時に「.multiply(0.001)」がついている
2. expressionというものを使っている

1について、これは掛け算を意味している呼び出したMODISのデータに0.0001をかけていることを意味している。（なぜ0.0001をかけるのがはわからない。EVIを計算する上での係数？）

2について、複雑な計算をする場合は「expression」を使ったほうがよいとしている。expressionの第一引数は計算式を示し、第二引数は第一引数で使われる変数を定義する。
演算子等はここを参考（https://developers.google.com/earth-engine/image_math）

下記にEVIを示す。

###Phase3

var image = ee.Image('MODIS/006/MOD09GA/2012_03_09').multiply(0.0001);
var lon = -94.84497
var lat = 39.01918
var z = 8
var Params = {min:0,max:1}
var evi = image.expression(
  '2.5 * (nir -red) / (nir + 6 *red -7.5 * blue + 1)',
  {
    red:image.select('sur_refl_b01'),
    nir:image.select('sur_refl_b02'),
    blue:image.select('sur_refl_b03')
  });
  
var ndvi = image.expression(
  'nir - red / nir + red',
  {
    nir:image.select('sur_refl_b02'),
    red:image.select('sur_refl_b01')
  })

Map.setCenter(lon,lat,z);
Map.addLayer(image.select(['sur_refl_b01','sur_refl_b04','sur_refl_b03']),
      Params,'MODIS bands 1/4/3')
Map.addLayer(evi,Params,'EVI')
Map.addLayer(ndvi,Params,'NDVI')

ぶっちゃけあんまり変わってない。
expressionを使ってNDVIを計算してみた。
個人的にexpressionを使ったほうがわかりやすいかもしれない。
※expression＝何かの計算をしている　ということがすぐに把握できるから。

Hillshade

ここでは陰影起伏を計算します。
また、GEE上での自作関数やfor文についても触れます。

###Phase1
コードは下記の通り

function radians(image){
  return image.toFloat().multiply(Math.PI).divide(180);
  }
  
function hillshade(az,ze,slope,aspect){
  var azimuth = radians(ee.Image(az));
  var zenith = radians(ee.Image(ze));
  return azimuth.subtract(aspect).cos()
    .multiply(slope.sin())
    .multiply(zenith.sin())
    .add(
      zenith.cos().multiply(slope.cos()));
}

var terrain = ee.Algorithms.Terrain(ee.Image('CGIAR/SRTM90_V4'));
var slope = radians(terrain.select('slope'));
var aspect = radians(terrain.select('aspect'));

Map.setCenter(-127.767,46.852,11);
for (var i = 0; i < 360; i+=60){
  Map.addLayer(hillshade(i,60,slope,aspect),{},i+'deg');
}

そもそもHillshadeを知らないので、調べました。
http://www.gsi.go.jp/bousaichiri/hillshademap.html
https://www.esrij.com/gis-guide/other-dataformat/elevation-data/
何に使われるというと、「尾根線、谷線の判別や断層の判読」にできるとのことです。
また、「地形の騎乗を立体的に表現したい」場合に適しているとのことです。

###Phase2
今までとは違い、関数が使われています。

function radians(img) {
  return img.toFloat().multiply(Math.PI).divide(180);
}

これは度からラジアンに変換する関数ですね。

function hillshade(az,ze,slope,aspect){
  var azimuth = radians(ee.Image(az));
  var zenith = radians(ee.Image(ze));
  return azimuth.subtract(aspect).cos()
    .multiply(slope.sin())
    .multiply(zenith.sin())
    .add(
      zenith.cos().multiply(slope.cos()));
}

これは、Hillshadeを計算するための関数のようです。
Hillshadeの計算式は下記の通りです
Hillshade = cos(Azimuth - Aspect) * sin(Slope) * sin(Zenith) + cos(Zenith) * cos(Slope)
これを計算するための関数のようです。

var terrain = ee.Algorithms.Terrain(ee.Image('CGIAR/SRTM90_V4'));
var slope = radians(terrain.select('slope'));
var aspect = radians(terrain.select('aspect'));

ee.Algorithms.Terrainとはなんでしょうか。
APIリファレンスを確認すると
「地形DEMから勾配、アスペクト、単純なHillShadeを計算する」とあります。
勾配は、傾斜角度、アスペクトは傾斜方向を意味します。
傾斜角度は、山岳域等は大きくなる傾向になり、都市域のような平坦な場所では小さくなる傾向になると想像は付きます。
アスペクト（傾斜方向）は0度（真北）～90度（真東）～180度（真南）～270度（真西）～360度（真北）と表すみたいです。

var slopeとvar aspectはslopeとaspectを選択（select）し、それらを度数からラジアンに変換して取得しているみたいですね。

for (var i = 0; i < 360; i += 60) {
  Map.addLayer(hillshade(i, 60, slope, aspect), {}, i + ' deg');
}

このコードは60度ごとにMap.addLayerをしています。ただ一つわからないのが中括弧がある部分です。これが何を示しているのかがよくわかりません。
もし分かる人がいたらコメント欄で教えていただけるとありがたいです。。。

300度のときの陰影起伏結果を示します。

Phase3

汚いコードになりましたが、下記です。

function radians(image){
  return image.toFloat().multiply(Math.PI).divide(180);
  }
  
function hillshade(az,ze,slope,aspect){
  var azimuth = radians(ee.Image(az));
  var zenith = radians(ee.Image(ze));
  return azimuth.subtract(aspect).cos()
    .multiply(slope.sin())
    .multiply(zenith.sin())
    .add(
      zenith.cos().multiply(slope.cos()));
}

function shadow(hillshadeimager){
}

var terrain = ee.Algorithms.Terrain(ee.Image('CGIAR/SRTM90_V4'));
var slope = radians(terrain.select('slope'));
var aspect = radians(terrain.select('aspect'));
var lat = 140;
var lon = 35;
var z = 8;

Map.setCenter(lat,lon,z);
for (var i = 0; i < 360; i+=60){
  var hillshadeimage = hillshade(i,60,slope,aspect)
  var image_mask = hillshadeimage.lt(0.1)
  var masks = image_mask.eq(1);
  var maskedComposite = masks.updateMask(masks);
  Map.addLayer(maskedComposite,{palette:['000000','FF0000']},i+'deg')
}

var hillshadeimage_0 = hillshade(0,60,slope,aspect)
var hillshadeimage_60 = hillshade(60,60,slope,aspect)
var hillshadeimage_120 = hillshade(120,60,slope,aspect)
var hillshadeimage_180 = hillshade(180,60,slope,aspect)
var hillshadeimage_240 = hillshade(240,60,slope,aspect)
var hillshadeimage_300 = hillshade(300,60,slope,aspect)
var hillshadeimage_360 = hillshade(360,60,slope,aspect)

var hillshadeimage_0_mask = hillshadeimage_0.lt(0.1)
var masks_0 = hillshadeimage_0_mask.eq(1);
var maskedComposite_0 = masks_0.updateMask(masks_0);

var hillshadeimage_60_mask = hillshadeimage_60.lt(0.1)
var masks_60 = hillshadeimage_60_mask.eq(1);
var maskedComposite_60 = masks_60.updateMask(masks_60);

var hillshadeimage_120_mask = hillshadeimage_120.lt(0.1)
var masks_120 = hillshadeimage_120_mask.eq(1);
var maskedComposite_120 = masks_120.updateMask(masks_120);

var hillshadeimage_180_mask = hillshadeimage_180.lt(0.1)
var masks_180 = hillshadeimage_180_mask.eq(1);
var maskedComposite_180 = masks_180.updateMask(masks_180);

var hillshadeimage_240_mask = hillshadeimage_240.lt(0.1)
var masks_240 = hillshadeimage_240_mask.eq(1);
var maskedComposite_240 = masks_240.updateMask(masks_240);

var hillshadeimage_300_mask = hillshadeimage_300.lt(0.1)
var masks_300 = hillshadeimage_300_mask.eq(1);
var maskedComposite_300 = masks_300.updateMask(masks_300);

var result = maskedComposite_0.eq(maskedComposite_60)
                              .eq(maskedComposite_120)
                              .eq(maskedComposite_180)
                              .eq(maskedComposite_240)
                              .eq(maskedComposite_300)
Map.addLayer(result,{palette:['000000','00BFFF']},'result')

影となっている箇所のみ、マスクをかけるということをしました。
影か、影じゃないかの判別は適当にしきい値をきめました。
Hillshadeの値をみてみたところ、だいたい1~0の値をもち、0.1以下だと影になってそうだと判明したため、0.1以下をマスクとして処理しています。
resultには、0度、60度、120度、180度、240度、300度の全方位すべて影になっているところのみマスクとする処理をしています。
個人的には、山岳域がたくさんあるような箇所はマスク処理されるだろうと予想していましたが、ほとんどありませんでした。

##Landcover Cleanup
ここでは膨張圧縮処理や再投影などを学びます。

###Phase1
コードは下記の通り

var SCALE = 500;

var image1 = ee.Image('MODIS/051/MCD12Q1/2001_01_01');
var image2 = image1.select(['Land_Cover_Type_1']);
var image3 = image2.reproject('EPSG:4326',null,SCALE);
var image4 = image3.focal_mode()
var image5 = image4.focal_max(3).focal_min(5).focal_max(3)
var image6 = image5.reproject('EPSG:4326',null,SCALE);

var PALETTE = [
    'aec3d4', // water
    '152106', '225129', '369b47', '30eb5b', '387242', // forest
    '6a2325', 'c3aa69', 'b76031', 'd9903d', '91af40', // shrub, grass, savannah
    '111149', // wetlands
    'cdb33b', // croplands
    'cc0013', // urban
    '33280d', // crop mosaic
    'd7cdcc', // snow and ice
    'f7e084', // barren
    '6f6f6f'  // tundra
].join(','); 

var vis_params = {min: 0, max: 17, palette: PALETTE};

Map.setCenter(-113.41842, 40.055489, 6);
Map.addLayer(image2,vis_params,'IGBP classification');

Map.addLayer(image3,vis_params,'Reprojected')
Map.addLayer(image4,vis_params,'Mode')
Map.addLayer(image5,vis_params,'Smooth')
Map.addLayer(image6,vis_params,'Smooth')

###Phase2

var SCALE = 500;

とありますが、これは500m分解能にするための変数のようです。後段の処理で、再投影などの処理をするのですが、そのときにSCALEを指定しないと変わってしまいます。そのため、再投影時にMODISの分解能に合わせるため、ここでSCALEとして500を保持しています。

var image1 = ee.Image('MODIS/051/MCD12Q1/2001_01_01');
var image2 = image1.select(['Land_Cover_Type_1']);
var image3 = image2.reproject('EPSG:4326',null,SCALE);
var image4 = image3.focal_mode()
var image5 = image4.focal_max(3).focal_min(5).focal_max(3)
var image6 = image5.reproject('EPSG:4326',null,SCALE);

image1では、MCD12Q1を読み込んでいます。（MCD12Q1：https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/mcd12q1）

image2では、Land_Cover_Type_1を選択していますね。ちなみに選択したデータは「国際地球圏 - 生物圏プログラム（IGBP）の全球植生分類計画」というデータとなっています。詳細はこちらにアクセスしてください（https://lpdaac.usgs.gov/about/news_archive/modisterra_land_cover_types_yearly_l3_global_005deg_cmg_mod12c1）

image3ではEPSGをWGS84に再投影します。
なぜ、再投影するのかというと、今回読み込んでいるMCD12Q1はサンソン図法というものを採用しています。GoogleEarthで使われている投影はWGS84といわれるもの**（たぶん。だれか指摘してくださると幸いです。）**です。この投影が変わると歪んでしまうため、投影法をWGS84に再投影します。

image4では、スムージングフィルタをかけています。

imge5では、膨張・収縮フィルタをかけています。

image6では、いままでかけてきたフィルタを施してきたデータにして、再投影をしています。（image3ですでに再投影しているので、する必要はないと思うのですが、再投影しています。）

var PALETTE = [
    'aec3d4', // water
    '152106', '225129', '369b47', '30eb5b', '387242', // forest
    '6a2325', 'c3aa69', 'b76031', 'd9903d', '91af40', // shrub, grass, savannah
    '111149', // wetlands
    'cdb33b', // croplands
    'cc0013', // urban
    '33280d', // crop mosaic
    'd7cdcc', // snow and ice
    'f7e084', // barren
    '6f6f6f'  // tundra
].join(',');

色を16進数で指定しています。今回使用している「Land_Cover_Type_1」は16分類のデータです。
公式サイトを確認すると、下記のように定義されています。

Land Cover	Class	Color
Fill Value	255	#000000
Water	0	#000080
Evergreen Needleleaf Forest	1	#008000
Evergreen Broadleaf Forest	2	#00FF00
Deciduous Needleleaf Forest	3	#99CC00
Deciduous Broadleaf Forest	4	#99FF99
Mixed Forest	5	#339966
Closed Shrubland	6	#993366
Open Shrubland	7	#FFCC99
Woody Savannas	8	#CCFFCC
Savannas	9	#FFCC00
Grasslands	10	#FF9900
Permanent Wetlands	11	#006699
Croplands	12	#FFFF00
Urban and Built-Up	13	#FF0000
Cropland/Natural Vegetation Mosaic	14	#999966
Snow and Ice	15	#FFFFFF
Barren or Sparsely Vegetated	16	#808080
Unclassified	254	#000000

どうやら色指定の16進数は違うみたいですね。

PALETTEの最後にjoin(',')をしていますが、これはなくても動作します。
なぜjoinしているのかがよくわかりませんでした。

###Phase3
このコードは特に思いつくことがなかったため、なし。

##Reduce_Region
ここでは画像表示ではなく、ee.Geometry.Rectangle()とee.Reducer.max()の使い方を学び、Console上に出力させる方法を学びます。
###Phase1
コードは下記の通り

var Image = ee.Image('CGIAR/SRTM90_V4')

var poly = ee.Geometry.Rectangle([-109.05, 41, -102.05, 37])
var max = Image.reduceRegion({
  reducer:ee.Reducer.max(),
  geometry:poly,
  scale:200,
  })
print(max)

###Phase2

var poly = ee.Geometry.Rectangle([-109.05, 41, -102.05, 37])

GISでいうところのシェイプを設定している。
sampleで指定している緯度経度は米国のコロラド州、デンバーあたりを示している。

var max = Image.reduceRegion({
  reducer:ee.Reducer.max(),
  geometry:poly,
  scale:200,
  })

はじめての関数「Image.reduceRegion」を使っている。
これは指定した地域（もしくは複数の地域）の統計値を計算する関数のようです。
画像領域のピクセル値の統計を取得するときに使います。
この中はmax()としているので最大値を取得しますが、他にもmean(),
median()などなど、たくさんあるみたいです。
指定する領域はgeometryで決定します。

scaleはおそらく解像度のことを指しているのだと思います。

###Phase3
Import機能を使いました。
下記コードになります。

var Image = ee.Image('CGIAR/SRTM90_V4')
var lat = 35.2
var lon = 138.45
var z = 10
Map.setCenter(lon,lat,z)

var max = Image.reduceRegion({
  reducer:ee.Reducer.max(),
  geometry:Fuji_poly,
  scale:200,
  })

var min = Image.reduceRegion({
  reducer:ee.Reducer.min(),
  geometry:Fuji_poly,
  scale:200,
  })

var mean = Image.reduceRegion({
  reducer:ee.Reducer.mean(),
  geometry:Fuji_poly,
  scale:200,
  })

print(max)
print(min)
print(mean)

コード上にはかいていないのですが　、polyの範囲をGeometry Importsでマウスより指定しました。
場所は富士山周辺です。
富士山周辺の標高が「最も高い標高値」、「最も低い標高値」、「標高値の平均」の３つを求めました。

Canny Edge Detectoror

ここではCanny法によるエッジ抽出を学びます。
###Phase1
下記コードとなります。

var Image = ee.Image('LANDSAT/LT05/C01/T1_TOA/LT05_031034_20110619')
var ndvi = Image.normalizedDifference(['B4','B3'])
var canny = ee.Algorithms.CannyEdgeDetector(ndvi,0.7)
canny = canny.updateMask(canny)
Map.setCenter(-101.05259, 37.93418, 13);
Map.addLayer(ndvi, {min:0,max:1},'Landsat NDVI')
Map.addLayer(canny,{min:0,max:1,palette:'FF0000'},'Canny Edges')

###Phase2

var canny = ee.Algorithms.CannyEdgeDetector(ndvi,0.7)

「CannyEdgeDetectory」という新しい関数がでてきました。これはCanny法によるエッジ検出をしているのだと思われます。（参考：https://ja.wikipedia.org/wiki/%E3%82%A8%E3%83%83%E3%82%B8%E6%A4%9C%E5%87%BA#%E3%82%AD%E3%83%A3%E3%83%8B%E3%83%BC%E6%B3%95_(Canny_edge_detector)）

canny = canny.updateMask(canny)

updateMaskを使っていますね。
これはピクセルがゼロ以外のもの（つまり1）をマスク画像とし、他の値のものは透明に設定するという関数です。
cannyにて計算されたエッジ検出画像は0か1かの二値画像のため、canny.updateMaskをすることでマスク画像を作成することが簡単みたいですね。
衛星リモートセンシングではマスク画像は多用すると思うので、updateMaskは重要だと感じます。

###Phase3
今回使用しているデータは、雲領域がありました。光学データだと、付き合い続ける嫌な領域だと思います。その雲領域を視覚化したいと思います。

また、サンプルで使用しているLandsat5のデータの諸元は下記の通りです。

Name	Resolution	Wavelength	Description
B1	30 meters	0.45 - 0.52 µm	Blue
B2	30 meters	0.52 - 0.60 µm	Green
B3	30 meters	0.63 - 0.69 µm	Red
B4	30 meters	0.76 - 0.90 µm	Near infrared
B5	30 meters	1.55 - 1.75 µm	Shortwave infrared 1
B6	30 meters	10.40 - 12.50 µm	Thermal Infrared 1. Resampled from 60m to 30m.
B7	30 meters	2.08 - 2.35 µm	Shortwave infrared 2
BQA			Landsat Collection 1 QA Bitmask (See Landsat QA page)

このうち、BQAは以下のようになっています。

• Bit 0: Designated Fill
  • 0: No
  • 1: Yes
• Bit 1: Designated Pixel
  • 0: No
  • 1: Yes
• Bits 2-3: Radiometric Saturation
  • 0: No bands contain saturation
  • 1: 1-2 bands contain saturation
  • 2: 3-4 bands contain saturation
  • 3: 5 or more bands contain saturation
• Bit 4: Cloud
  • 0: No
  • 1: Yes
• Bits 5-6: Cloud Confidence
  • 0: Not Determined / Condition does not exist.
  • 1: Low, (0-33 percent confidence)
  • 2: Medium, (34-66 percent confidence)
  • 3: High, (67-100 percent confidence)
• Bits 7-8: Cloud Shadow Confidence
  • 0: Not Determined / Condition does not exist.
  • 1: Low, (0-33 percent confidence)
  • 2: Medium, (34-66 percent confidence)
  • 3: High, (67-100 percent confidence)
• Bits 9-10: Snow / Ice Confidence
  • 0: Not Determined / Condition does not exist.
  • 1: Low, (0-33 percent confidence)
  • 2: Medium, (34-66 percent confidence)
  • 3: High, (67-100 percent confidence)

以上より、4ビット目を見れば雲があるかないかを判定することが可能です。この4ビット目を見てマスク画像を作っていきましょう。

下記コードです。

var LANDSAT = ee.Image('LANDSAT/LT05/C01/T1_TOA/LT05_031034_20110619')


var getQABits = function(Image, start, end, newName) {
    var pattern = 0;
    for (var i = start; i <= end; i++) {
       pattern += Math.pow(2, i);
    }
    return Image.select([0], [newName])
                  .bitwiseAnd(pattern)
                  .rightShift(start);
};

var cloud_check = function(Image) {
  var QA = Image.select(['BQA']);
  return getQABits(QA, 4,4, 'Cloud').eq(1);
};

var maskClouds = function(Image) {
  var c = cloud_check(Image);
  return c.updateMask(c);
};

var mask_cloud = maskClouds(LANDSAT);

var ndvi = LANDSAT.normalizedDifference(['B4','B3']);
var canny = ee.Algorithms.CannyEdgeDetector(ndvi,0.7);
canny = canny.updateMask(canny);
Map.setCenter(-101.05259, 37.93418, 13);
Map.addLayer(ndvi, {min:0,max:1},'Landsat NDVI');
Map.addLayer(canny,{min:0,max:1,palette:'FF0000'},'Canny Edges');
print(mask_cloud)
Map.addLayer(mask_cloud,{palette:'FFF000'},'Cloud_mask');

var getQABits = function(Image, start, end, newName) {
    var pattern = 0;
    for (var i = start; i <= end; i++) {
       pattern += Math.pow(2, i);
    }
    return Image.select([0], [newName])
                  .bitwiseAnd(pattern)
                  .rightShift(start);
};

getQABitsが指定したビットの値を求める関数です。start とendの間のビットを調べます。

var maskClouds = function(Image) {
  var c = cloud_check(Image);
  return c.updateMask(c);
};

var mask_cloud = maskClouds(LANDSAT);

ここでマスク画像を作っています。

Map.addLayer(mask_cloud,{palette:'FFF000'},'Cloud_mask');

最後にaddLayerです。黄色のパレットに指定して、表示させています。

bit値の取得は衛星データでよくすることなので、多用していくと思いました。

##Center Pivot Irrigation Detector
ここではセンターピボットの灌漑抽出を行います。ee.Kernel.circle ee.Kernel.square convolveを学びます。
##Phase1
下記コードです。

Map.setCenter(-106.06,37.71,12)

var palette = [
  'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718',
  '74A901', '66A000', '529400', '3E8601', '207401', '056201',
  '004C00', '023B01', '012E01', '011D01', '011301'];
  
var image = ee.Image('LANDSAT/LC08/C01/T1_TOA/LC08_034034_20170608');
var ndvi = image.normalizedDifference(['B5','B4']);
var ndviParams = {min:0.0,max:1.0,palette:palette}

Map.addLayer(ndvi,ndviParams,'Landsat NDVI')

var farmSize = 500;
var circleKernel = ee.Kernel.circle(farmSize,'meters');
var squareKernel = ee.Kernel.square(farmSize,'meters');
var circles = ndvi.convolve(circleKernel);
var squares = ndvi.convolve(squareKernel);
var diff = circles.subtract(squares);

var diff = diff.abs().multiply(100).toByte();
var max = diff.focal_max({radius:farmSize * 1.8,units:'meters'});
var local = diff.where(diff.neq(max),0);
var thresh = local.gt(2);

var peaks = thresh.focal_max({kernel:circleKernel});
Map.addLayer(peaks.updateMask(peaks),{palette:'FF3737'},'Kernel Peaks')

var canny = ee.Algorithms.CannyEdgeDetector(ndvi,0);
canny = canny.gt(0.3);

var inner = ee.Kernel.circle(farmSize -20, 'meters',false,-1)
var outer = ee.Kernel.circle(farmSize + 20, 'meters',false,-1)
var ring = outer.add(inner,true);

var centers = canny.convolve(ring).gt(0.5).focal_max({kernel:circleKernel});
Map.addLayer(centers.updateMask(centers),{palette:'4285FF'},'Ring Centers');

##Phase2
円形の灌漑を検出するプログラムになります。

Map.setCenter(-106.06, 37.71, 12);

// A nice NDVI palette.
var palette = [
  'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718',
  '74A901', '66A000', '529400', '3E8601', '207401', '056201',
  '004C00', '023B01', '012E01', '011D01', '011301'];

// Just display the image with the palette.
var image = ee.Image('LANDSAT/LC08/C01/T1_TOA/LC08_034034_20170608');
var ndvi = image.normalizedDifference(['B5','B4']);

Map.addLayer(ndvi, {min: 0, max: 1, palette: palette}, 'Landsat NDVI');

Map.setCenterで指定している緯度経度を確認すると、円形の耕作地が確認できます。
センターピボットというのは主に米国のグレートプレーンズやサウジアラビア、エジプトなどの乾燥地域にて行われている農法のようです。詳しくはここを参照してください。（https://ja.wikipedia.org/wiki/%E3%82%BB%E3%83%B3%E3%82%BF%E3%83%BC%E3%83%94%E3%83%9C%E3%83%83%E3%83%88）

paletteではNDVIのカラーパレットを指定しています。

衛星データはLandsat-8の2017年06月08二値のデータを使用しているみたいですね。

その後、NDVIを算出してMap.addLayerをしています。

var farmSize = 500;
var circleKernel = ee.Kernel.circle(farmSize, 'meters');
var squareKernel = ee.Kernel.square(farmSize, 'meters');
var circles = ndvi.convolve(circleKernel);
var squares = ndvi.convolve(squareKernel);
var diff = circles.subtract(squares);

耕作地のサイズを半径500と定義しています。
（なぜ500かというと、予め今回ターゲットsとしている猟奇は半径500mの円形農場で占めていることが知られているみたいです。）
今回はじめてKernelというものがでてきました。Kernelとは、円形や四角形などの領域のことだと理解しています。（geometryのregionのような理解）
今回は、KernelをCircle。つまり円形のKernelと定義し、大きさは500を指定しています。同様に、四角形のKernelも定義しています。

その後、convolveをしています。convolveは畳み込みの関数で、引数はKernelのみです。
ndvi画像とcircleの畳み込み、ndvi画像とsquareの畳み込みを作成しています。

最後にcircleで畳み込んだ結果からsquareで畳み込んだ結果を引きます。
わかりにくいので、画像を示します。

まずはNDVI画像

次にNDVIをCircleで畳み込んだ画像

次にNDVIをSquareで畳み込んだ画像

最後にCircleで畳み込んだ結果からsquareで畳んだ結果を引いた画像

var diff = diff.abs().multiply(100).toByte();
var max = diff.focal_max({radius:farmSize * 1.8,units:'meters'});
var local = diff.where(diff.neq(max),0);
var thresh = local.gt(2);

ここでは、円の中心を求めていると思われます。

diffでは、まず絶対値（abs()）にし、100を乗算し（multiply(100)）、データ型をByte型（.toByte()）にしています。

maxでは、focal_maxを使って膨張させています。膨張させた結果は以下の通りです。

localではwhereを使った条件付き演算を行っています。もしdiffの値とmaxの値が一緒じゃない場合、0にします。

threrhでは、2よりも大きい値を1にし、それ以外を0にしています。いわゆる二値化画像を作成しています。

var peaks = thresh.focal_max({kernel:circleKernel});
Map.addLayer(peaks.updateMask(peaks),{palette:'FF3737'},'Kernel Peaks')

var canny = ee.Algorithms.CannyEdgeDetector(ndvi,0);
canny = canny.gt(0.3);

peaksでは、focal_maxを使って、circleKernelで指定した大きさ（今回は500）まで膨張させています。さらにそれをupdateMaskをしてマスク画像としています。マスクのカラーパレットは赤色を指定しています。

面白い処理だなと思いました。何かに応用させてみたいです。
しかし、注意が必要なのはthreshにて求めた値を、circleKernelの値まで膨張させているので、もしセンターピボットの半径が500ではない場合はその都度farmSizeを変更させる必要があります。

cannyではndviデータのエッジ抽出をし、gtにて二値画像を作成しています。

var inner = ee.Kernel.circle(farmSize -20, 'meters',false,-1)
var outer = ee.Kernel.circle(farmSize + 20, 'meters',false,1)
var ring = outer.add(inner,true);

再び、Kernel.circleがでてきました。今回は第四引数まで書いてあります。
一度ここで整理してみましょう。

• 第一引数：半径
  • 生成するカーネルの半径
• 第二引数：単位
  • 「ピクセル」または「メートル」を指定。カーネルがメートル単位で指定されている場合は、ズームレベルが変更されたときにサイズが変更されます。
• 第三引数：正規化（デフォルト：true）
  • カーネル値を1になるように正規化する
• 第四引数：絶対値（デフォルト：1）
  • この値で各値を調整する

innerには500-20としたfarmSizeにし、単位をmeterにし、正規化はしていません。
outerには500+20としたfarmSizeにし、単位をmeterにし、正規化はしていません。ただし絶対値としています。

その後、addをしています。
outerにinnerを追加しカーネルを正規化しています。
(このあたりはだんだんわからなくなってきました。。。)

var centers = canny.convolve(ring).gt(0.5).focal_max({kernel: circleKernel});
Map.addLayer(centers.updateMask(centers), {palette: '4285FF'}, 'Ring centers');

ここは、半径が480から520に収まっている円形の部分のみを青色で抽出しています。

※innerとouterの処理の部分は理解がおいつかず。。。

###Phase3
このコードは特に思いつくことがなかったため、なし。

##Clamp
ここではClampを学びます。
###Phase1
下記コードです。

var Image = ee.Image('CGIAR/SRTM90_V4');
var clamped = Image.clamp(1000,2000);
var vis_Params = {min:0,max:4300};

Map.setCenter(-121.753, 46.855, 9);
Map.addLayer(Image,vis_Params,'Full Stretch');
Map.addLayer(clamped,vis_Params,'Clamped');

##Phase2

var Image = ee.Image('CGIAR/SRTM90_V4');
var clamped = Image.clamp(1000,2000);
var vis_Params = {min:0,max:4300};

SRTMを読み込み、clamp関数を使っています。

clamp関数は、指定された範囲内に収まるように画像のすべてのバンドの値を固定してくれます。
引数は下記の通りです。

• 第一引数
  • 範囲内の最小値（Float）
• 第二引数
  • 範囲内の最大値（Float）

Map.setCenter(-121.753, 46.855, 9);
Map.addLayer(Image,vis_Params,'Full Stretch');
Map.addLayer(clamped,vis_Params,'Clamped');

指定した緯度経度は下記になります。
レーニア山国立公園になります。

実行結果は以下の通りです。

Full Stretch（Clampなし）

Clamped（Clampedあり）

Clampなしの場合、レーニア山の山頂付近がサチレーションを起こしているように見えます。
Clampありの場合はサチレーションはしていませんが、表現度が落ちているような気がします。また、西の方も一面が黒色になっているため、標高の高低が表現できませんでした。

###Phase3
このコードは特に思いつくことがなかったため、なし。

##Connected Pixel Count
ここではconnectedPixelCountを学びます。

###Phase1
下記コードです。

var Image = ee.Image('MODIS/006/MOD09GA/2012_03_09')
      .select('sur_refl_b01')
      .multiply(0.0001)
      
var bright = Image.gt(0.3);

var conn = bright.connectedPixelCount({
  maxSize:30,
  eightConnected:true
});

var smallClusters = conn.lt(30)

Map.setCenter(-107.24304, 35.78663, 8);
Map.addLayer(Image,{min:0,max:1},'original')
Map.addLayer(smallClusters.updateMask(smallClusters),
          {min:0,max:1,palette:'FF0000'},'cc')

###Phase2

var Image = ee.Image('MODIS/006/MOD09GA/2012_03_09')
      .select('sur_refl_b01')
      .multiply(0.0001)
      
var bright = Image.gt(0.3);

MODISの画像を読み込み、さらに'sur_refl_b01'を選択しています。
これは波長が「620-670nm」のデータを指しており、”赤”を指しています。
さらに、multiplyで0.0001をしています。
前もmulitplyで0.0001を乗算しているデータがありました。なぜこれをかけるのだろうかわかりませんでしたが、datasetsを確認したらScaleに0.0001と書いてあったので、これのことだと思われます。

その後、gtで0.3以上のものを1それ以外を0とする二値画像を作成します。

var conn = bright.connectedPixelCount({
  maxSize:30,
  eightConnected:true
});

var smallClusters = conn.lt(30)

新しい関数が出てきました。connectedPixelCountです。
この関数について調べてみると、下記のようです。

各ピクセルに4または8近傍が含まれる画像を生成します。

• 第一引数:input_image
  • 入力画像
• 第二引数:maxSize（デフォルト：100）
  • 近傍の最大サイズ
• 第三引数:eightConnected（デフォルト：8（True））
  • 4近傍を使うのか8近傍を使うのか

maxSize（整数、デフォルト：100）：
近傍の最大サイズ（ピクセル単位）。

eightConnected（ブール、デフォルト：true）：
8近傍ではなく4近傍のルールを使用するかどうか。
※近傍についてはここを参照してください（https://www58.atwiki.jp/dooooornob/pages/53.html）

maxサイズを30としているため、30ピクセル以下の塊を強調表示することを意味しています。

その後、var smallClusters = conn.lt(30)にて30以下のピクセル値は1、それ以外は0という二値画像を作成しています。

Map.setCenter(-107.24304, 35.78663, 8);
Map.addLayer(![image.png](https://qiita-image-store.s3.amazonaws.com/0/284911/4d7d2374-e2bf-70f7-0d38-96959bac7c5d.png)
Image,{min:0,max:1},'original')
Map.addLayer(smallClusters.updateMask(smallClusters),
          {min:0,max:1,palette:'FF0000'},'cc')

指定された緯度経度では、ナバホ・ネイション
を指しています。
その後、addLayerにて画面表示しています。

###Phase3
このコードは特に思いつくことがなかったため、なし。

##HSV Pan Sharpening
ここではHSVを使ったパンシャープン法を学びます。
##Phase1
下記コードです。

var Image = ee.Image(ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
  .filterDate('2017-01-01','2017-12-31')
  .filterBounds(ee.Geometry.Point(-122.0808, 37.3947))
  .sort('CLOUD_COVER')
  .first());

var rgb = Image.select('B4','B3','B2')
var pan = Image.select('B8')

var huesat = rgb.rgbToHsv().select('hue','saturation')
var upres = ee.Image.cat(huesat,pan).hsvToRgb()

Map.setCenter(-122.0808, 37.3947, 14);

Map.addLayer(rgb,{max:0.3},'Original');
Map.addLayer(upres,{max:0.3},'Panshapened')

###Phase2

var Image = ee.Image(ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
  .filterDate('2017-01-01','2017-12-31')
  .filterBounds(ee.Geometry.Point(-122.0808, 37.3947))
  .sort('CLOUD_COVER')
  .first());

新しくImageCollectionがでてきました。
これは、画像の時系列データを得ることができます。今まで使っていたImageと異なり、たくさんのデータを扱うことができます。

ImageCollectionの中に、
filterDate
filterBounds
sort
first()
があります。

filterDate；
これは期間を示しています。
LANDSAT/LC08/C01/T1_TOAを使用するとしていますが、
その中でも、「2017-01-01」～「2017-12-31」のデータを指定しています。

filterBounds；
これは領域を示しています。
これを指定しないと、LANDSAT/LC08/C01/T1_TOAでかつ「2017-01-01」～「2017-12-31」のデータを表示してしまう事になり、全球レベルでのデータを表示してしまいます。
そのため、filterBoundsにて「-122.0808, 37.3947」の領域を表示するように設定していまs。
sort；
指定したプロパティ（今回だと’CLOUD_COVER’）で並び替えをします。デフォルトだと昇順となっています。
first()；
Collectionは複数の画像を表示します。その中でも最も最初に表示される　画像のみを表示します。

var rgb = Image.select('B4','B3','B2')
var pan = Image.select('B8')

var huesat = rgb.rgbToHsv().select('hue','saturation')

ここではパンシャープン処理をしています。
Landsat8のバンド8は、パンクロマチック画像で分解能が高いです。これを使って高解像度パンクロマチック画像を作成します。

var huesat = rgb.rgbToHsv().select('hue','saturation')で使用しているrgbToHsv()はRGB色空間からHSV色空間に変換しています。
また、その後.select('hue','saturation')としていますが、HSV色空間に変換すると、変換後のデータはhue,saturation,valueの値を返します。それらのデータをselectしています。

var upres = ee.Image.cat(huesat,pan).hsvToRgb()

ee.Image.catは画像を連結する関数です。huesatとpanを連結させています。しかし、連結させただけだとHSV色空間のままなので、hsvToRgb()をしてRGB色空間に変換しています。

Map.setCenter(-122.0808, 37.3947, 14);

Map.addLayer(rgb,{max:0.3},'Original');
Map.addLayer(upres,{max:0.3},'Panshapened')

画像を表示させています。

以下はパンシャープン前の画像です。

以下はパンシャープン後の画像です。

明らかに分解能が高くなったのが　確認できます。パンシャープンをするとNIR等が使えないため、NDVIなどの計算ができません。

###Phase3
このコードは特に思いつくことがなかったため、なし。

##Hough
ここではHough変換（直線検出）を学びます。
###Phase1
下記コードです。

var image = ee.Image('LANDSAT/LC08/C01/T1_TOA/LC08_033032_20170719')
var ndvi = image.normalizedDifference(['B5','B4'])
var img_Params = {bands:['B4','B3','B2'],max:0.3}
var ca_Params = {min:0,max:1,palette:'blue'}
var h_Params = {min:0,max:1,palette:'red'}
var canny = ee.Algorithms.CannyEdgeDetector({
  image: ndvi,
  threshold:0.4
}).multiply(255)

var h = ee.Algorithms.HoughTransform({
  image:canny,
  gridSize:256,
  inputThreshold:50,
  lineThreshold:100
})

Map.setCenter(-103.80140, 40.21729, 13);
Map.addLayer(image,img_Params,'source image')
Map.addLayer(canny.updateMask(canny),ca_Params,'canny')
Map.addLayer(h.updateMask(h),h_Params,'hough')

###Phase2

var image = ee.Image('LANDSAT/LC08/C01/T1_TOA/LC08_033032_20170719')
var ndvi = image.normalizedDifference(['B5','B4'])

画像を読み込み、NDVIを計算しています。

var img_Params = {bands:['B4','B3','B2'],max:0.3}
var ca_Params = {min:0,max:1,palette:'blue'}
var h_Params = {min:0,max:1,palette:'red'}

Map.addLayerを使うときのパラメータを定義しました。

var canny = ee.Algorithms.CannyEdgeDetector({
  image: ndvi,
  threshold:0.4
}).multiply(255)

canny法によるエッジ抽出を指定まs。
その後、multiplyで255を乗算しています。
なぜ255を乗算するのかというと、後続のhough変換の際に、255をかけないと使えないみたいです。

var h = ee.Algorithms.HoughTransform({
  image:canny,
  gridSize:256,
  inputThreshold:50,
  lineThreshold:100
})

今回のキモとなる関数です。入力画像にhough変換を適用します。hough変換はこちらを参照してください（https://ja.wikipedia.org/wiki/%E3%83%8F%E3%83%95%E5%A4%89%E6%8F%9B）

Map.setCenter(-103.80140, 40.21729, 13);
Map.addLayer(image,img_Params,'source image')
Map.addLayer(canny.updateMask(canny),ca_Params,'canny')
Map.addLayer(h.updateMask(h),h_Params,'hough')

画像を表示します。

canny法を適用した画像を表示します。

Hough変換を適用した画像を表示します。

###Phase3
このコードは特に思いつくことがなかったため、なし。

##MODIS
このコードは、**Canny Edge Detectoror
**のときにかなり似たコードを書いているので、Phase1だけ。。。
やっていることは、MODISのデータないにあるQAビットを取得するということです。雲領域を透過処理させています。

###Phase1

var getQABits = function(image,start,end,newName){
    var pattern = 0
    for(var i = start;i <= end; i++){
      pattern += Math.pow(2,i)
    }
    return image.select([0],[newName])
                .bitwiseAnd(pattern)
                .rightShift(start);
}

var image = ee.Image('MODIS/006/MOD09GA/2012_10_11');
var QA = image.select('state_1km')

var cloud = getQABits(QA,0,1,'cloud_state')
              .expression("b(0) == 1 || b(0) == 2");
              
var landWaterFlag = getQABits(QA,3,5,'land_water_flag');

var mask = landWaterFlag.neq(7).and(cloud.not())

Map.addLayer(image.updateMask(mask),
{
  bands:['sur_refl_b01','sur_refl_b04','sur_refl_b03'],
  min:-100,
  max:2000
},'MOD09GA 143'
);

似たようなことを前述して説明しているのであまり書くことはないです。。。

##Pixel
ここではPixelArea()を学びます。単位がピクセルのデータをメートルなどに変換するときに使います。
##Phase1

var image = ee.Image.pixelArea()
Map.setCenter(0,0,3)
var vis_Params = {min: 2e8, max: 4e8, opacity: 0.85}
Map.addLayer(image,vis_Params,'Pixel area')

###Phase2

Image.pixelArea()という新しい関数がでました。
これは面積の計算を行うために使います。通常衛星画像では、１ピクセルの単位はあくまでもピクセルのままで計算していることが多いです。ピクセル単位ではなく、メートル単位などの単位に変換するときに使用します。この面積の画像に損失の画像を乗算してから結果を合計すると、面積の測定値が得られます。例えば、1ピクセル100mのものを、1ピクセル1mにしたい場合は100を乗算する必要があります。そういうときは、ee.Image.pixelArea().multiply(100)`のようにします。

特定の場所における特定の土地被覆の増減を計算するときなどに使えるみたいです。
なお、returnとしてはImageが返ってきます。

この記事がわかりやしと思います（https://ourcodingclub.github.io/2018/11/26/earth-engine.html）

地球観測衛星画像は、衛星軌道により高緯度では歪みが生じます。歪みが生じることで1ピクセルの面積が減少していきます。（分解能が変わってくるということです）

（この説明文はあっているか不明。だれか指摘していただけるとありがたいです。）

###Phase3
このコードは特に思いつくことがなかったため、なし。

##Pixel Lon Lat
ここでは緯度と経度という名前を持つ画像を作成するpixelLonLat()を学びます。
###Phase1
下記コードです。

var image = ee.Image.pixelLonLat().multiply(60.0)
image = image.subtract(image.floor()).lt(0.05)
var grid = image.select('latitude').or(image.select('longitude'))
Map.setCenter(-122.09228, 37.42330, 12);
Map.addLayer(grid.updateMask(grid),{palette:'000000'},'Graticule')

###Phase2

ee.Image.pixelLonLat().multiply(60.0)

ee.Image.pixelLonLat()という新しい関数が出てきました。
これは、軽度と緯度という２つのバンドを持つ画像を作成する関数です。
今回は60度グリッドを作成するため、multiply(60)にて60を乗算しています。

image = image.subtract(image.floor()).lt(0.05)

ここで「image　ー　imageの0.05以下の値」をしています。pixelLonLatにて作成した緯度経度ラインと、それ以外とでは値に差があり、上述している数式を適用することで緯度経度のライン以外の者は0になるのかなと理解しております。

image.select('latitude').or(image.select('longitude'))
Map.setCenter(-122.09228, 37.42330, 12);
Map.addLayer(grid.updateMask(grid),{palette:'000000'},'Graticule')

ここでorを使って'latitude'と'longitude'を選択しています。

selectで'latitude'だけ選択したものと、'longitude'だけ選択したものと、上述しているor演算子にて両方選択したものを見比べてみます。

下記はlatitudeだけ選択した画像

下記はlongitudeだけ選択した画像

下記はo演算子にて両方選択した画像

orを使うことで、グリッド線がかけるとわかりました。

###Phase3
このコードは特に思いつくことがなかったため、なし。

##Polynomial1
ここでは衛星画像に多項式を適用させる.polynomialを学びます。
##Phase1
下記コードです。

var image = ee.Image('MODIS/006/MOD09GA/2012_03_09')
          .select(['sur_refl_b01','sur_refl_b04','sur_refl_b03'])
          .multiply(0.0001)

var adj = image .polynomial([-0.2,2,4,-1.2])
var vis_Params = {min:0,max:1}

Map.setCenter(-107.24304, 35.78663, 8);
Map.addLayer(image,vis_Params,'original')
Map.addLayer(adj,vis_Params,'adjusted')

##Phase2

var adj = image .polynomial([-0.2,2,4,-1.2])
var vis_Params = {min:0,max:1}

polynomialという新しい関数が出てきました。
これは画像に非線形コントラストを適用する関数のようです。
今回の場合、-0.2 + 2.4x - 1.2x^2を適用しています。

適用前の画像は以下の通りです。

適用後の画像は以下の通りです。

かなり見やすくなりました。

###Phase3
このコードは特に思いつくことがなかったため、なし。

##Zero
ここではゼロと交差する場所を検出するzeroCrossing()を学びます。
###Phase1
下記コードです。

var elev = ee.Image('CGIAR/SRTM90_V4')
var image = elev.subtract(1000).zeroCrossing()
var vis_Params = {min:0,max:1,opacity:0.5}
Map.setCenter(-121.68148, 37.50877, 13)
Map.addLayer(image,vis_Params,'Crossing 1000m')
var exact = elev.eq(1000)
Map.addLayer(exact.updateMask(exact),{palette:'red'},'Exactly 1000m')

###Phase2

var image = elev.subtract(1000).zeroCrossing()
新しい関数が出てきました。.zeroCrossing()です。これは、値がゼロになった箇所を検出する関数です。
今回の場合は、elevから1000を引き、zeroCrossingをしているので標高が1000の部分のみ検出しています。等高線みたいな感じですね。

Map.addLayer(exact.updateMask(exact),{palette:'red'},'Exactly 1000m')

このコードではちょうど1000mのところのみを検出しています

###Phase3
このコードは特に思いつくことがなかったため、なし。

#おわりに
写経してきましたが、それでもまだ自由自在に使えるようには全くなれていませんし、不十分なところがあります。自分で問題を仮設してコーディングするなどして応用力を身に着けていきたいです。