GEE代码学习 day11

使用表达式计算 BAI

Martin （1998） 开发了烧伤面积指数 （BAI） 来帮助描绘烧伤疤痕和评估烧伤严重程度。它依赖于留下灰烬和木炭的火;不会产生灰烬或木炭的火灾以及灰烬和木炭已被冲走或覆盖的旧火将无法被很好地检测到。BAI 计算每个像素到燃烧区域往往相似的光谱参考点的光谱距离。远离此参考的像素（例如，健康的植被）将具有非常小的值，而靠近此参考的像素（例如，来自火的木炭）将具有非常大的值。

这个方程中有两个常数：ρcr 是红色波段的常数，等于 0.1;ρcnir 用于 NIR 波段，等于 0.06。要检查燃烧指数，请加载 2013 年的图像，该图像显示了加利福尼亚州内华达山脉的边缘大火。我们将使用 Landsat 8 来探索这场火灾。

// Examine the true-color Landsat 8 images for the 2013 Rim Fire. 
var burnImage =  ee.ImageCollection('LANDSAT/LC08/C02/T1_TOA') 
                    .filterBounds(ee.Geometry.Point(-120.083, 37.850)) 
                    .filterDate('2013-09-15', '2013-09-27') 
                    .sort('CLOUD_COVER') 
                    .first();  
Map.centerObject(ee.Geometry.Point(-120.083, 37.850), 11);  
var rgbParams = {  bands: ['B4', 'B3', 'B2'], min: 0, max: 0.3 }; 
Map.addLayer(burnImage, rgbParams, 'True-Color Burn Image');

检查此图像的真彩色显示。你能发现火势吗？如果没有，BAI 可能会有所帮助。与 EVI 一样，使用表达式在 Earth Engine 中计算 BAI

// Calculate BAI. 
var bai = burnImage.expression(  '1.0 / ((0.1 - RED)**2 + (0.06 - NIR)**2)', { 'NIR': burnImage.select('B5'), 'RED': burnImage.select('B4'), });

显示结果。

// Display the BAI image. 
var burnPalette = ['green', 'blue', 'yellow', 'red'];  
Map.addLayer(bai, { min: 0, max: 400, palette: burnPalette }, 'BAI');

9.2 使用矩阵代数处理图像---利用矩阵代数的线性变换

缨帽改造

TC 变换是一类变换，在绘制图形时，它们看起来像一顶带有流苏的羊毛帽。最常见的实施方式是用来最大限度地分离小麦的不同生长阶段，小麦是一种具有重要经济意义的作物。随着小麦的生长，田地从裸露的土壤发展到绿色植物的发育，再到黄色的植物成熟，田间收获。将大面积的许多田地的这些阶段分开是流苏帽改造的最初目的。

Kauth 和 Thomas 通过将他们变换空间的第一轴定义为平行于下图中的土壤线来找到 R。选择第一个列以指向土壤的长轴以及从美国伊利诺伊州给定点的 Landsat 影像得出的值。第二列被选为与第一列正交，并指向他们所谓的“绿色事物”，即绿色植被。第三列与前两列正交，指向“黄色的东西”，例如，成熟的小麦和其他禾本科作物。最后一列与前三列正交，在原始派生中称为 “nonesuch” — 即类似于噪声。

其中 p0 是原始 p × 1 像素向量（该特定像素的 p 波段值堆栈作为数组），矩阵 R 是新空间的正交基，其中每列彼此正交（因此，RT 是它的转置），输出 p1 是该像素的旋转值堆栈

已为每颗 Landsat 卫星推导出 R 矩阵，包括 Landsat 5 （Crist 1985）、Landsat 7、（Huang et al. 2002）、Landsat 8 （Baig et al. 2014） 等。我们可以使用 Arrays 在 Earth Engine 中实现此转换。具体来说，让我们创建一个新脚本，并为 Landsat 5 的专题制图器 （TM） 工具制作一个 TC 系数数组：

///// 
// Manipulating images with matrices 
/////  
// Begin Tasseled Cap example. 
var landsat5RT = ee.Array([  [0.3037, 0.2793, 0.4743, 0.5585, 0.5082, 0.1863], 
                             [-0.2848, -0.2435, -0.5436, 0.7243, 0.0840, -0.1800], 
                             [0.1509, 0.1973, 0.3279, 0.3406, -0.7112, -0.4572], 
                             [-0.8242, 0.0849, 0.4392, -0.0580, 0.2012, -0.2768], 
                             [-0.3280, 0.0549, 0.1075, 0.1855, -0.4357, 0.8085], 
                             [0.1084, -0.9022, 0.4120, 0.0573, -0.0251, 0.0238] ]);  
print('RT for Landsat 5', landsat5RT);

请注意，我们刚刚创建的结构是一个包含六个列表的列表，然后将其转换为 Earth Engine ee.Array 对象。6 x 6 的值数组对应于 TM 仪器的六个非热波段值的线性组合：波段 1-5 和 7。

由于这些系数适用于卫星反射率（大气顶部）处的 TM 传感器，因此我们将访问云量较少的 Landsat 5 场景。我们将访问 Landsat 5 影像的集合，对其进行过滤，然后按增加的云量进行排序，并获取第一个影像

You can search for “Odessa, WA, USA” in the search bar.

// Define a point of interest in Odessa, Washington, USA. 
var point = ee.Geometry.Point([-118.7436019417829,  47.18135755009023]); Map.centerObject(point, 10);  
// Filter to get a cloud free image to use for the TC. 
var imageL5 = ee.ImageCollection('LANDSAT/LT05/C02/T1_TOA')  .filterBounds(point) .filterDate('2008-06-01', '2008-09-01') .sort('CLOUD_COVER') .first();  
//Display the true-color image. 
var trueColor = {  bands: ['B3', 'B2', 'B1'], min: 0, max: 0.3 }; 
Map.addLayer(imageL5, trueColor, 'L5 true color');

要进行矩阵乘法，首先将输入图像从多波段图像（对于每个波段，每个像素存储一个值）转换为数组图像。数组图像是一种更高维度的图像，其中每个像素存储一个波段的值数组。

var bands = ['B1', 'B2', 'B3', 'B4', 'B5', 'B7'];  
// Make an Array Image, with a one dimensional array per pixel. 
// This is essentially a list of values of length 6, 
// one from each band in variable 'bands.' 
var arrayImage1D = imageL5.select(bands).toArray(); 
 // Make an Array Image with a two dimensional array per pixel, 
// of dimensions 6x1. This is essentially a one column matrix with 
// six rows, with one value from each band in 'bands.' 
// This step is needed for matrix multiplication (p0). 
var arrayImage2D = arrayImage1D.toArray(1);

1 指的是列（Earth Engine 中的“第一个”轴），为 p0 创建一个 6 行 x 1 列的数组，接下来，我们使用 matrixMultiply 函数完成流苏帽线性变换的矩阵乘法，然后使用 arrayProject 和 arrayFlatten 函数将结果转换回多波段图像

.arrayFlatten() 方法用于将数组图像的波段展平成一个多波段图像。您提供了一个包含波段名称的二维数组（尽管在这个上下文中它实际上是一个列表的列表，其中每个内部列表只包含一个元素），并且 GEE 会根据这个数组中的名称来命名新的多波段图像中的波段。这些波段对应于 tasseled cap 变换的不同组件。

定义了一个名为 vizParams 的对象，用于指定图像的可视化参数。这里，您选择了 brightness、greenness 和 wetness 这三个波段进行可视化，并为它们分别设置了最小值和最大值。然后，您使用 Map.addLayer() 方法将变换后的图像添加到地图上，并为其指定了一个名称 'TC components'。

//Multiply RT by p0. 
var tasselCapImage = ee.Image(landsat5RT)  // Multiply the tasseled cap coefficients by the array// made from the 6 bands for each pixel. 
                       .matrixMultiply(arrayImage2D) // Get rid of the extra dimensions. 
                       .arrayProject([0]) // Get a multi-band image with TC-named bands. 
                       .arrayFlatten( [  ['brightness', 'greenness', 'wetness', 'fourth', 'fifth', 'sixth' ] ]);
var vizParams = {  bands: ['brightness', 'greenness', 'wetness'], 
                   min: -0.1, max: [0.5, 0.1, 0.1] }; 
Map.addLayer(tasselCapImage, vizParams, 'TC components');