GEE代码学习 day12

TC变换

TC变换要进行矩阵乘法，首先将输入图像从多波段图像（对于每个波段，每个像素存储一个值）转换为数组图像。数组图像是一种更高维度的图像，其中每个像素存储一个波段的值数组。（第 IV 部分将更详细地介绍和讨论数组图像。您将使用波段 1-5 和 7 以及 toArray 函数：

在Google Earth Engine (GEE) 中，即使您不显式地指定数组维度（即不传递参数给 toArray()），生成的数组影像在内部实际上是以二维数组（或矩阵）的形式存储的，但这里的“二维”并不是指每个像素都是一个具有多行多列的矩阵，而是指每个像素的值是一个数组，且这个数组可以被视为一个1xN（1行N列）或Nx1（N行1列）的矩阵，其中N是选定波段的数量。

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.' 
// toArray(1)方法中的参数1指定了新数组的第二个维度的大小（即列数），而第一个维度（行数）则自动由原始数组的长度（在这个例子中是6）确定。 
var arrayImage2D = arrayImage1D.toArray(1);

我们使用 matrixMultiply 函数完成流苏帽线性变换的矩阵乘法，然后使用 arrayProject 和 arrayFlatten 函数将结果转换回多波段图像：

“tasseled cap”变换（也称为K-L变换或主成分变换的特定形式）

Landsat 5的反射率影像（landsat5RT）

arrayProject方法用于减少数组影像的维度。在这个例子中，传递[0]作为参数意味着您想要保留数组的第一个维度（即，如果您有一个NxM的数组，arrayProject([0])会保留N行，但将M列“压平”为一个一维数组）。

arrayFlatten方法用于将数组影像转换为一个多波段影像，其中每个波段都有一个指定的名称。在这个例子中，您指定了六个波段名称：brightness（亮度）、greenness（绿度）、wetness（湿度）、fourth、fifth和sixth。这些名称对应于tasseled cap变换后的波段。

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] ]);
//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');

PCA变换

要演示将 PCA 应用于影像的实际应用，请导入 Landsat 8 TOA 影像，并将其命名为 imageL8。首先，我们将它转换为数组图像：

.sort('CLOUD_COVER'): 这一步将剩余的影像按照CLOUD_COVER属性进行排序。CLOUD_COVER是一个元数据属性，表示影像中云覆盖的百分比。排序默认是升序，即云覆盖最少的影像会被放在前面。

// Begin PCA example.  
// Select and map a true-color L8 image. 
var point = ee.Geometry.Point([-118.7436019417829,  47.18135755009023]); 
Map.centerObject(point, 10);
var imageL8 = ee.ImageCollection('LANDSAT/LC08/C02/T1_TOA')  
                .filterBounds(point) 
                .filterDate('2018-06-01', '2018-09-01') 
                .sort('CLOUD_COVER') 
                .first();  
// 真彩色波段
var trueColorL8 = {  bands: ['B4', 'B3', 'B2'], 
                     min: 0, 
                     max: 0.3 }; 
Map.addLayer(imageL8, trueColorL8, 'L8 true color');  
// 定义主成分分析波段. 
var PCAbands = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B10',  'B11'];  
// Convert the Landsat 8 image to a 2D array for the later matrix 
// toArray()方法将这些波段转换为一个2D数组。这个2D数组是后续进行矩阵计算（如PCA）的基础。
var arrayImage = imageL8.select(PCAbands).toArray();

在下一步中，使用 reduceRegion 方法和 ee.Reducer.covariance 函数来计算图像的统计数据（在本例中为波段的协方差）。

使用reduceRegion方法对arrayImage应用了一个协方差（covariance）归约器。reduceRegion方法会对影像的指定区域（默认情况下是整个影像）应用归约器，并返回一个包含归约结果的字典。

// Calculate the covariance using the reduceRegion method. 
var covar = arrayImage.reduceRegion({  reducer: ee.Reducer.covariance(), maxPixels: 1e9 });  // Extract the covariance matrix and store it as an array. 
var covarArray = ee.Array(covar.get('array'));

使用 ee.Array.get 函数提取协方差矩阵并将其存储为数组。现在我们有一个基于图像的协方差矩阵，我们可以执行特征分析来计算执行 PCA 所需的特征向量。为此，我们将使用 eigen 函数。

//Compute and extract the eigenvectors 
var eigens = covarArray.eigen();

特征函数同时输出特征向量和特征值。由于我们需要 PCA 的特征向量，因此我们可以使用数组的 slice 函数来提取它们。特征向量存储在 1 轴的第 0 个位置。

var eigenVectors = eigens.slice(1, 1);

我们使用这些 eigenVector 和我们之前创建的 arrayImage 执行矩阵乘法。这与我们对流苏帽组件使用的流程相同。每次乘法都会产生一个主成分。

// Perform matrix multiplication 
var principalComponents = ee.Image(eigenVectors)  
                            .matrixMultiply(arrayImage.toArray(1));

最后，转换回多波段图像并显示第一个主成分 （pc1）：

var pcImage = principalComponents  // Throw out an unneeded dimension, [[]] -> []. 
             .arrayProject([0]) // Make the one band array image a multi-band image, [] -> image. 
             .arrayFlatten([ ['pc1', 'pc2', 'pc3', 'pc4', 'pc5', 'pc6', 'pc7', 'pc8'] ]);  // Stretch this to the appropriate scale. 
Map.addLayer(pcImage.select('pc1'), {}, 'pc1');

arrayProject() 方法用于从多维数组中选择特定的维度（或轴）

arrayFlatten() 方法用于将数组图像的波段“展平”成一个多波段图像。

确定要绘制的波段后，输入每个波段的最小值和最大值，确保它们的顺序正确。

//The min and max values will need to change if you map different bands or locations. 
var visParamsPCA = {  bands: ['pc1', 'pc3', 'pc4'], 
                      min: [-455.09, -2.206, -4.53], 
                      max: [-417.59, -1.3, -4.18] };  
Map.addLayer(pcImage, visParamsPCA, 'PC_multi');

端元混合分析是一种技术，用于从多光谱或高光谱遥感影像中分离出不同的地表覆盖类型（如植被、水体、裸土等）的贡献比例。

如果我们考虑数据集中的单个像素（例如，对应于 Landsat 像素的 30 × 30 m 空间），则它可能表示地面上的多个物理对象。因此，像素的光谱特征是该空间中存在的每个物体的“纯”光谱的混合物。例如，考虑森林的 Landsat 像素。像素的光谱特征是树木、林下、树木投射的阴影以及通过树冠可见的土壤块的混合。线性谱解混模型基于这个假设（Schultz et al. 2016， Souza 2005）。纯光谱（称为端元）来自土地覆被类，例如水、裸露土地和植被。这些端元表示来自地面特征（例如仅裸露地面）的纯光谱的光谱特征。目标是求解 ƒ 的以下方程，即像素中端元分数的 P × 1 向量：

S 是 B × P 矩阵，其中 B 是波段数，列是 P 纯端元光谱，p 是当有 B 波段时× 1 像素向量的 B（图 9.7）。我们知道 p，我们可以定义端元谱来得到 S，这样我们就可以求解 ƒ

// Specify which bands to use for the unmixing. 
var unmixImage = imageL8.select(['B2', 'B3', 'B4', 'B5',  'B6', 'B7']);

第一步是定义端成员，以便我们可以定义 S。我们将通过计算围绕净土覆盖区域划定的多边形中的平均光谱来实现这一点。将地图缩放到裸露土地、植被和水体均质区域的位置（机场可用作合适的位置）。将 Landsat 8 图像可视化为伪彩色合成：

// Use a false color composite to help define polygons of 'pure' land cover. 
Map.addLayer(imageL8, 
              { bands: ['B5', 'B4', 'B3'], 
               min: 0.0, max: 0.4 }, 
               'false color');

为了更快地渲染，您可能需要注释掉之前添加到地图中的图层。通常，执行此操作的方法是在纯土地覆被区域周围绘制多边形，以定义这些土地覆被的光谱特征。如果您想自己执行此操作，请按以下步骤操作。使用几何图形绘制工具，通过选择多边形工具，然后单击 + 新建图层，创建三个新图层（因此，P = 3）。在第一个图层中，数字化纯裸土地周围的多边形;在第二层中，制作一个纯植被的多边形;在第三个图层中，创建一个 Water 多边形。将导入分别命名为 bare、water 和 veg。您需要使用设置（齿轮图标）来重命名几何图形。

您还可以使用此代码指定裸露、水和植被的预定义区域。这仅适用于此示例。

// Define polygons of bare, water, and vegetation. 
var bare = /* color: #d63000 */ ee.Geometry.Polygon(  [ [  [-119.29158963591193, 47.204453926034134], [-119.29192222982978, 47.20372502078616], [-119.29054893881415, 47.20345532330602], [-119.29017342955207, 47.20414049800489] ] ]), water = /* color: #98ff00 */ ee.Geometry.Polygon( [ [  [-119.42904610218152, 47.22253398528318], [-119.42973274768933, 47.22020224831784], [-119.43299431385144, 47.21390604625894], [-119.42904610218152, 47.21326472446865], [-119.4271149116908, 47.21868656429651], [-119.42608494342907, 47.2217470355224] ] ]), veg = /* color: #0b4a8b */ ee.Geometry.Polygon( [ [  [-119.13546041722502, 47.04929418944858], [-119.13752035374846, 47.04929418944858], [-119.13966612096037, 47.04765665820436], [-119.13777784581389, 47.04408900535686] ] ]);

//Print a chart.
 var lcfeatures = ee.FeatureCollection([  ee.Feature(bare, {label: 'bare'}), ee.Feature(water, {label: 'water'}), ee.Feature(veg, {label: 'vegetation'}) ]);  
print( ui.Chart.image.regions({ image: unmixImage, regions: lcfeatures, reducer: ee.Reducer.mean(), scale: 30, seriesProperty: 'label', xLabels: [0.48, 0.56, 0.65, 0.86, 1.61, 2.2] }) .setChartType('LineChart') .setOptions({ title: 'Image band values in 3 regions', hAxis: { title: 'Wavelength' }, vAxis: { title: 'Mean Reflectance' } }));

xLabels 代码行在六个波段中每个波段的光谱中点处获取每个多边形（要素）的平均值。数字 （[0.48， 0.56， 0.65， 0.86， 1.61， 2.2]） 表示这些光谱中点。

使用 reduceRegion 方法计算您创建的多边形内每个波段的平均值。请注意，reduceRegion 的返回值是一个数字 Dictionary，用于汇总多边形中的值，输出按带级名称编制索引。通过在计算平均值后调用 values 函数来获取 List 形式的平均值。请注意，值按字母数字顺序返回结果，这些顺序按键排序。这之所以有效，是因为 B2 − B7 已经按字母数字排序，但在它们尚未排序的情况下将不起作用。在这些情况下，请指定乐队名称列表，以便您首先按已知顺序获取它们。

// Get the means for each region.
var bareMean = unmixImage 
              .reduceRegion(ee.Reducer.mean(), bare, 30).values(); 
var waterMean = unmixImage  
              .reduceRegion(ee.Reducer.mean(), water, 30).values(); 
var vegMean = unmixImage  
              .reduceRegion(ee.Reducer.mean(), veg, 30).values();

这三个列表中的每一个都代表一个平均频谱向量，它是上面定义的 S 矩阵的列之一。通过沿 1 轴（列）连接向量，将它们堆叠成 6 × 3 个末端成员数组：

// Stack these mean vectors to create an Array. 
var endmembers = ee.Array.cat([bareMean, vegMean,  waterMean], 1); 
print(endmembers);

现在将 6 波段输入图像转换为每个像素都是 1D 向量 （toArray） 的图像，然后转换为每个像素都是 6 × 1 矩阵的图像 （toArray（1））。这将创建 p，以便我们可以为每个像素求解上述方程。

// Convert the 6-band input image to an image array. 
var arrayImage = unmixImage.toArray().toArray(1);

现在我们已经准备好了一切，对于每个像素，我们求解 ƒ 的方程

// Solve for f. 
var unmixed = ee.Image(endmembers).matrixSolve(arrayImage);

最后，将结果从二维数组图像转换为一维数组图像 （arrayProject），然后转换为更熟悉的零维多波段图像 （arrayFlatten）。这与我们在前面几节中使用的方法相同。这三个波段对应于 ƒ 中裸露、植被和水分数的估计值：

显示结果，其中裸露为红色，植被为绿色，水为蓝色

// Convert the result back to a multi-band image. 
var unmixedImage = unmixed  
                   .arrayProject([0]) 
                   .arrayFlatten([ ['bare', 'veg', 'water'] ]);
Map.addLayer(unmixedImage, {}, 'Unmixed');