【GEE】6、在 Google 地球引擎中构建各种遥感指数

// Call in image by unique ID.
var bioClim = ee.Image("WORLDCLIM/V1/BIO");
print(bioClim, "bioClim");
1
2
3

从 print 语句中我们可以看到每个变量都存储为一个带区。ee.Image.select()我们可以使用该函数拉取我们感兴趣的特定波段。选择后，我们会将结果剪辑到我们感兴趣的区域并重命名，以便我们记住“bio12”代表什么。

// Select specific band, clip to region of interest, and rename.
var precip = bioClim
  .select("bio12")
  .clip(geometry)
  .rename("precip");
 
Map.addLayer(precip, {min:240 , max:1300}, "precip", false);
1
2
3
4
5
6
7

我们可以从生成的图像中看出，当您从东向西穿越大平原时，降水量会出现可预测的减少。我们将在该地区的西部看到更多的灌溉农业是有道理的。我们将添加 Cropland 数据层来测试这个假设。

// Call in collection base on unique ID and print.
var crops = ee.ImageCollection("USDA/NASS/CDL");
print(crops, "crops");
1
2
3

图像集合包含多个具有相同属性的图像。因此，处理图像集合在很大程度上依赖于过滤功能。过滤很有用，但应用过滤器不会更改影像的数据类型。过滤器总是返回一个图像集合对象。这很好，但是不同的 GEE 对象有不同的功能。通常您需要将图像集合转换为图像。最常见的方法之一是应用reducer。与其他遥感软件相比，减速器是使 GEE 脱颖而出的元素之一。虽然它们非常强大，但我们不会在本模块中详细介绍它们的功能。

// Filter collection by a date range and pull a specific image from the collection.
var crop2016 = crops
  .filterDate("2016-01-01","2016-12-30")
  .select(['cropland']);
print(crop2016, "crop2016");
1
2
3
4
5

正如我们所料，该crop2016对象仍然是一个图像集合对象，即使它只包含一个图像。

2016 年农田数据集的打印声明。
 

我们可以简单地通过它包含的唯一 ID 调用此图像，然后使用该Image.select()函数拉取我们感兴趣的特定波段，而不是应用 reducer。我们可以通过打印对象和读取元数据来查看图像 ID。

// Call the specific image by the unique ID, select band of interest, and clip to region of interest.
var cropType2016 =  ee.Image("USDA/NASS/CDL/2016")
  .select("cropland")
  .clip(geometry);
print(cropType2016, "cropType2016");
Map.addLayer(cropType2016, {}, "cropType2016", false);
1
2
3
4
5
6

在这里，我们调用图像，选择感兴趣的波段，并将图像剪辑到我们的输出区域。

我们可以使用这两个数据集来查看土地耕地和该地区平均降水量的大图模式。如果你真的想深入研究这些数据，你可以看看种植的作物类型。您需要通过搜索数据集从 USDA 层提取元数据。如上所述，数据集中有 255 种不同的土地覆被类型，因此它很快就会变得不堪重负。下面是我们感兴趣的地区常见作物的表格。

作物的圆形区域表示使用中心枢轴灌溉的田地。考虑到该地区的干旱程度以及流经该地区的主要河流数量，中心支点使用的水很可能来自地表以下，即奥加拉拉含水层。不圆的田地仍然可以灌溉。虽然效率可能较低，但许多地区使用漫灌来浇灌作物。如果没有对作物和该地的生长条件的广泛了解，则很难根据这两层来确定田地是否进行了洪水灌溉。我们希望农田和非农田之间光谱反射率值的差异有助于我们做出区分。

您可以打开这两个图层并调整透明度以查看降水量和农田之间的重叠。

通过增加降水图像的透明度，我们可以直观地评估大面积地理区域作物类型与年降水量之间的关系。
 

// Evapotranspiration: call in image collection by unique ID and print to understand data structure.
var et = ee.ImageCollection("MODIS/006/MOD16A2")
  .filterBounds(geometry);
print(et, "modisET");
1
2
3
4

调用新数据集并打印它是一种很好的做法，这样您就可以在尝试过滤之前查看数据结构和图像可用性。即使请求与数据集中存在的信息不匹配，GEE 也会对数据集应用过滤器，因此您需要了解您的要求是否可行。

从元数据中，我们建立了一个有意义的日期范围，并且我们知道我们对哪些波段感兴趣。现在我们将过滤集合。

// Filter by data, select band, apply reducer to get single image and clip image.
var et2 = et
  .filterDate("2016-07-01","2016-09-30")
  .select(['ET'])
  .median()
  .clip(geometry);
 
Map.addLayer(et2, {max: 357.18 , min: 29.82}, "evapotransparation", false);
1
2
3
4
5
6
7
8

蒸发蒸腾层作为图像集提供。我们在这里应用一个中值缩减器，因为我们想要代表该地区生长季节的广义值。中值缩减器将获取一堆图像，计算每个单元格的中值，并将所有这些单元格编译回单个图像。应用中值缩减器后，您的图像不再代表单个时间点，而是一段时间内的信息集合。

我们也将对“初级生产力总产值”层重复此过程。

// Gross Primary Productivity.
var GPP = ee.ImageCollection("UMT/NTSG/v2/MODIS/GPP");
print(GPP, "modisGPP");
 
var GPP2 = GPP
  .filterDate("2016-07-01","2016-09-30")
  .select(['GPP'])
  .median()
  .clip(geometry);
 
Map.addLayer(GPP2, {max: 1038.5 , min: 174.5}, "gross Primary Productivity", false);
1
2
3
4
5
6
7
8
9
10
11

之前的两个数据集来自 MODIS 图像。我们使用的 NDVI 层来自不同的传感器 Landsat。虽然两种传感器之间存在许多差异，但主要因素是空间分辨率（像素区域）和时间分辨率（图像重新捕获周期）。Landsat 的空间分辨率要好得多（30 m，而 MODIS 为 500 m），但时间分辨率要低得多（返回时间为 16 天，而 MODIS 为 1 天）。始终建议比较多个传感器，因为每个传感器都有优点和缺陷。您需要了解哪种工具最适合您的应用程序。

// Pulling in NDVI data.
var LS8NDVI = ee.ImageCollection("LANDSAT/LC8_L1T_32DAY_NDVI")
print(LS8NDVI, "ls8NDVI");
// Filter, reduce, then clip.
var LS8NDVI2 = LS8NDVI
  .filterDate("2016-07-01","2016-09-30")
  .median()
  .clip(geometry);
 
Map.addLayer(LS8NDVI2, {max: 0.66 , min: -0.23}, "NDVI", false);
1
2
3
4
5
6
7
8
9
10

加载地图上的三幅图像需要一些时间来探索在整个景观中识别的视觉模式。观察图像之间反射率差异的好地方是在科罗拉多州的乔斯附近。您可以像在 Google 地图中一样在搜索栏中搜索“Joes, Colorado”来找到它。在乔斯方圆几英里内有一些明显的中心枢纽灌溉区、代表自然​​景观的未开发土地和河岸带。

min我们使用定义的可视化参数和 将所有图像添加到地图中max。这些值是通过应用拉伸和记录脚本中的最小值和最大值来创建的。这是一种使数据中的差异更容易看到的方法，并且不会改变数据中的基础值。

科罗拉多州乔斯周围地区的 NDVI 图像。明亮区域代表具有高 NDVI 值的位置，并与灌溉土地和河岸地区相关。
 

// Call in Landsat 8 surface reflectance Image Collection, filter by region and cloud cover, and reduce to single image.
var LS8_SR2 = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
.filterDate("2017-07-01","2017-09-30")
  .filterBounds(geometry)
  .filterMetadata('CLOUD_COVER', 'less_than', 20)
  .median();
 
print(LS8_SR2,'LS8_SR2');
1
2
3
4
5
6
7
8

按元数据过滤是一个新命令。您可以通过检查图像的文档和函数本身来了解它是如何工作的。

有多种方法可以定义索引。我们将在本节中介绍三个。其中两种方法需要定义单独的波段，一种需要从图像中选择波段。我们的第一步是从波段中定义单个图像。

// Pull bands from the LS8 image and save them as images.
var red = LS8_SR2.select('B4').rename("red");
var nir = LS8_SR2.select('B5').rename("nir");
print(red);
1
2
3
4

归一化差异植被指数 (NDVI) 是一种非常常见的植被测量方法。使用以下等式计算：

NDVI = (nir - red)/(nir + red)

生成指数的第一种方法是使用数学函数来重现 NDVI 方程。您不能在 GEE 中单独使用符号 (+,-,/,*) 作为数学表达式。相反，您需要写出操作。

// Generate NDVI based on predefined bands.
var ndvi1 = nir.subtract(red).divide(nir.add(red)).rename('ndvi');
 
Map.addLayer(ndvi1, {max: 0.66 , min: -0.23},'NDVI1',false);
1
2
3
4

直接调用数学函数很快就会变得很麻烦。另一种选择是依赖内置函数。在我们的例子中，有一个特征normalizedDifference()函数ee.Image。我们在 Landsat 8 图像集合中应用了缩减器，因此我们得到了可以使用的图像。因此，调用此函数就像将其指向正确的波段一样简单。

// Define LS8 NDVI using built in normalized difference function.
var ndvi2 = LS8_SR2.normalizedDifference(['B5', 'B4']);
Map.addLayer(ndvi2, {max: 0.66 , min: -0.23},'NDVI2',false);
1
2
3

可用于生成索引的最后一种方法是构建表达式。鉴于 GEE 中内置了归一化差分函数，因此将这种方法用于 NDVI 是没有意义的。我们在这里以它为例。对于更复杂的指数，例如流苏帽转换或增强植被指数 (EVI)，该过程是必要的。

// Use the expression function to generate NDVI.
var ndvi3 = LS8_SR2.expression("(nir - red) / (nir+red)", {
  "nir" : nir,
  "red" : red }
)
Map.addLayer(ndvi3, {max: 0.66 , min: -0.23},'NDVI3',false);
1
2
3
4
5
6

表达式要复杂得多。我们在字符串中定义数学表达式。然后我们使用字典来引用表达式中的变量所指的内容。我们已经将近似推断（‘B5’）和红带（‘B4’）保存为变量，因此我们可以直接调用它们。

生成索引的三种不同方法的相同 NDVI 值 该图像以内布拉斯加州奥加拉拉镇为中心。中心的黑暗特征是麦康纳湖。

如果你比较三个 NDVI 图像，你会发现它们都是一样的；到同一终点的不同路线。您使用什么方法取决于您拥有的数据集以及您尝试创建的索引。那里有数百种选择。索引数据库是一个很好的参考，可以了解每个特定传感器索引数据库的可能性。

// Call in image by unique ID.
var bioClim = ee.Image("WORLDCLIM/V1/BIO");
print(bioClim, "bioClim");
 
// Select specific band, clip to region of interest, and rename.
var precip = bioClim
  .select("bio12")
  .clip(geometry)
  .rename("precip");
 
Map.addLayer(precip, {min:240 , max:1300}, "precip", false);
 
// Call in collection base on unique ID and print.
var crops = ee.ImageCollection("USDA/NASS/CDL");
print(crops, "crops");
 
// Filter collection by a date range and pull a specific image from the collection.
var crop2016 = crops
  .filterDate("2016-01-01","2016-12-30")
  .select(['cropland']);
print(crop2016, "crop2016");
 
// Call the specific image by the unique ID, select band of interest, and clip to region of interest.
var cropType2016 =  ee.Image("USDA/NASS/CDL/2016")
  .select("cropland")
  .clip(geometry);
print(cropType2016, "cropType2016");
Map.addLayer(cropType2016, {}, "cropType2016", false);
 
// Evapotranspiration: call in image collection by unique ID and print to understand data structure.
var et = ee.ImageCollection("MODIS/006/MOD16A2")
  .filterBounds(geometry);
print(et, "modisET");
 
// Filter by data, select band, apply reducer to get single image and clip image.
var et2 = et
  .filterDate("2016-07-01","2016-09-30")
  .select(['ET'])
  .median()
  .clip(geometry);
 
Map.addLayer(et2, {max: 357.18 , min: 29.82}, "evapotransparation", false);
 
// Gross Primary Productivity.
var GPP = ee.ImageCollection("UMT/NTSG/v2/MODIS/GPP");
print(GPP, "modisGPP");
 
var GPP2 = GPP
  .filterDate("2016-07-01","2016-09-30")
  .select(['GPP'])
  .median()
  .clip(geometry);
 
Map.addLayer(GPP2, {max: 1038.5 , min: 174.5}, "gross Primary Productivity", false);
 
// Pulling in NDVI data.
var LS8NDVI = ee.ImageCollection("LANDSAT/LC8_L1T_32DAY_NDVI");
print(LS8NDVI, "ls8NDVI");
// Filter, reduce, then clip.
var LS8NDVI2 = LS8NDVI
  .filterDate("2016-07-01","2016-09-30")
  .median()
  .clip(geometry);
 
Map.addLayer(LS8NDVI2, {max: 0.66 , min: -0.23}, "NDVI", false);
 
// Call in Landsat 8 surface reflectance Image Collection, filter by region and cloud cover, and reduce to single image.
var LS8_SR2 = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
.filterDate("2017-07-01","2017-09-30")
  .filterBounds(geometry)
  .filterMetadata('CLOUD_COVER', 'less_than', 20)
  .median();
 
print(LS8_SR2,'LS8_SR2');
 
// Pull bands from the LS8 image and save them as images.
var red = LS8_SR2.select('B4').rename("red");
var nir = LS8_SR2.select('B5').rename("nir");
print(red);
 
// Generate NDVI based on predefined bands.
var ndvi1 = nir.subtract(red).divide(nir.add(red)).rename('ndvi');
 
Map.addLayer(ndvi1, {max: 0.66 , min: -0.23},'NDVI1',false);
 
// Define LS8 NDVI using built in normalized difference function.
var ndvi2 = LS8_SR2.normalizedDifference(['B5', 'B4']);
Map.addLayer(ndvi2, {max: 0.66 , min: -0.23},'NDVI2',false);
 
// Use the expression function to generate NDVI.
var ndvi3 = LS8_SR2.expression("(nir - red) / (nir+red)", {
  "nir" : nir,
  "red" : red }
)
Map.addLayer(ndvi3, {max: 0.66 , min: -0.23},'NDVI3',false);
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95