Skip to main content

Lab 3B - Transformations

We've gone over indices to highlight unique characteristics in our imagery by utilizing the bands outside of the visible spectrum.

Linear transforms are linear combinations of input pixel values. These can result from a variety of different strategies, but a common theme is that pixels are treated as arrays of band values, and we can use these arrays to create weighted values for specific purposes.

3B.1.1 - Tasseled cap (TC)​

Based on observations of agricultural land covers in the NIR-red spectral space, Kauth and Thomas (1976) devised a rotational transform of the form

p1=RTp0p_1 = R^T p_0

where:

  1. p0p_0 is the original pixel vector (a stack of the p band values as an Array)

  2. p1p_1 is the rotated pixel

  3. R is an orthonormal basis of the new space (inverse of RTR^T). Kauth and Thomas found R by defining the first axis of their transformed space to be parallel to the soil line in the following chart, then used the Gram-Schmidt process to find the other basis vectors.

Tassled Cap

Assuming that R is available, one way to implement this rotation in Earth Engine is with arrays. Specifically, make an array of TC coefficients. Since these coefficients are for the TM sensor, get a less cloudy Landsat 5 scene. To do the matrix multiplication, first convert the input image from a multi-band image to an array image in which each pixel position stores an array. Do the matrix multiplication, then convert back to a multi-band image.

var point = ee.Geometry.Point([-106.81, 39.19]);
var coefficients = 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]
]);  
var landsat = ee.ImageCollection("LANDSAT/LT05/C01/T1_TOA")
var image = ee.Image(landsat
  .filterBounds(point)
  .filterDate('2008-06-01', '2008-09-01')
  .sort('CLOUD_COVER')
  .first());
var bands = ['B1', 'B2', 'B3', 'B4', 'B5', 'B7'];
// Make an Array Image,  with a 1-D Array per pixel.
var arrayImage1D =  image.select(bands).toArray();
// Make an Array Image  with a 2-D Array per pixel, 6x1.
var arrayImage2D = arrayImage1D.toArray(1);  
var componentsImage = ee.Image(coefficients)
				.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.centerObject(point, 9)
Map.addLayer(componentsImage, vizParams, 'TC components');

Tassled Cap

Question 3: Upload the a tasseled cap image from a point near Blacksburg, VA and interpret the output. what are some of the values that you can extract when using Inspector? Are the results meaningful?

3B.1.2 - Principal Component Analysis (PCA)​

Like the Tasseled Cap transform, the PCA transform is a rotational transform in which the new basis is orthonormal, but the axes are determined from statistics of the input image, rather than empirical data. Specifically, the new basis is the eigenvector of the image's variance-covariance matrix. As a result, the principal components are uncorrelated. To demonstrate, use the Landsat 8 image converted to an array image. Use the reduceRegion() method to compute statistics (band covariances) for the image.

A reducer is an object that tells Earth Engine what statistic to compute. Note that the result of the reduction is an object with one property, an array, that stores the covariance matrix. The next step is to compute the eigenvectors and eigenvalues of that covariance matrix. Since the eigenvalues are appended to the eigenvectors, slice the two apart and discard the eigenvectors. Perform the matrix multiplication, as with the TC components. Finally, convert back to a multi-band image and display the first PC.

Use the layer manager to stretch the result appropriately. What do you observe? Try displaying some of the other principal components. The image parameters in the code chunk below are built specifically for Principal Component 1.

var point = ee.Geometry.Point([-80.42, 37.22]);
var landsat = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA")
var image = ee.Image(landsat
  .filterBounds(point)
  .filterDate('2013-11-17', '2014-03-27')
  .sort('CLOUD_COVER')
  .first());
var bands = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B10', 'B11'];
var arrayImage =  image.select(bands).toArray();  
var covar = arrayImage.reduceRegion({
  reducer: ee.Reducer.covariance(),
  maxPixels: 1e9
});
var covarArray = ee.Array(covar.get('array'));  
var eigens = covarArray.eigen();  
var eigenVectors = eigens.slice(1, 1);  
var principalComponents = 			ee.Image(eigenVectors).matrixMultiply(arrayImage.toArray(1));  
var pcImage = principalComponents      
				.arrayProject([0])    
// Make the one band  array image a multi-band image, [] -> image.    
				.arrayFlatten(
          [['pc1', 'pc2', 'pc3', 'pc4', 'pc5', 'pc6', 'pc7', 'pc8']]
        );      
// Customize the visual parameters for PC1
var imageVisParam = {
  "opacity":1,
  "bands":["pc1"],
  "min":-420,"max":-400,
  "gamma":1};
Map.centerObject(point, 10);    
Map.addLayer(pcImage.select('pc1'), imageVisParam, 'PC');

Tassled Cap

Question 4: How much did you need to stretch the results to display outputs for principal component 2?

Display and upload images of each the other principal components, stretching each band as needed for visual interpretation and indicating how you selected the min and max values.

How do you interpret each PC band? On what basis do you make that interpretation?

3B.1.3 - Spectral Unmixing​

The linear spectral mixing model is based on the assumption that each pixel is a mixture of "pure" spectra. The pure spectra, called endmembers, are from land cover classes such as water, bare land, vegetation. The goal is to solve the following equation for f, the Px1 vector of endmember fractions in the pixel:

Sf=pSf = p

where S is a BxP matrix in which the columns are P pure endmember spectra (known) and p is the Bx1 pixel vector when there are B bands (known). In this example, B=6B= 6:

The first step is to get the endmember spectra, which we can do by computing the mean spectra in polygons around regions of pure land cover. In this example, we will use a location in northern Washington State and use the geometry tools to select homogeneous areas of bare land, vegetation and water.

Using the geometry drawing tools, make three layers clicking + new layer. In the first layer, digitize a polygon around pure bare land, in the second layer make a polygon of pure vegetation, and in the third layer, make a water polygon. Name the imports bare, water and, and vegetation, respectively.

Note: For a starting point, we included some basic polygons but feel free to replace in a region of your choice.

var point = ee.Geometry.Point([-123.25, 48.11]);
var landsat = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA")
var image = ee.Image(landsat
  .filterBounds(point)
  .filterDate('2013-06-17', '2014-03-27')
  .sort('CLOUD_COVER')
  .first());
// Polygons of bare earth, water and vegetation
var bare = 
    ee.Geometry.Polygon(
        [[[-123.2370707334838, 48.1151452657945],
          [-123.2370707334838, 48.11351208612645],
          [-123.23410957473136, 48.11351208612645],
          [-123.23410957473136, 48.1151452657945]]], null, false);
var water = 
    ee.Geometry.Polygon(
        [[[-123.2748188020549, 48.12059599002954],
          [-123.2748188020549, 48.118074835535865],
          [-123.2673086168132, 48.118074835535865],
          [-123.2673086168132, 48.12059599002954]]], null, false);
var vegetation = 
    ee.Geometry.Polygon(
        [[[-123.27462568300582, 48.11533866992809],
          [-123.27462568300582, 48.114163936320416],
          [-123.27215805071212, 48.114163936320416],
          [-123.27215805071212, 48.11533866992809]]], null, false);
var unmixImage = image.select(['B2', 'B3', 'B4', 'B5', 'B6', 'B7']);

Check the polygons you made by charting mean spectra in them using Chart.image.regions():

Your chart should look something like:

Tassled Cap

var unmixImage = image.select(['B2', 'B3', 'B4', 'B5', 'B6', 'B7']);  
print(Chart.image.regions(unmixImage, ee.FeatureCollection([
    ee.Feature(bare, {label: 'bare'}), 
    ee.Feature(water, {label: 'water'}),
    ee.Feature(vegetation, {label: 'vegetation'})]), 
  ee.Reducer.mean(), 30, 'label', [0.48, 0.56, 0.65, 0.86, 1.61, 2.2]));
"""
var = (
print(Chart.image.regions(unmixImage, ee.FeatureCollection([
    ee.Feature(bare, {label: 'bare'}), 
    ee.Feature(water, {label: 'water'}),
    ee.Feature(vegetation, {label: 'vegetation'})]), 
  ee.Reducer.mean(), 30, 'label', [0.48, 0.56, 0.65, 0.86, 1.61, 2.2])))
  """

Use the reduceRegion() method to compute mean spectra in the polygons you made. Note that the return value of reduceRegion() is a Dictionary, with reducer output keyed by band name. Get the means as a list by calling values(). Each of these three lists represents a mean spectrum vector. Stack the vectors into a 6x3 Array of endmembers by concatenating them along the 1-axis (or columns direction).

Turn the 6-band input image into an image in which each pixel is a 1D vector (toArray()), then into an image in which each pixel is a 6x1 matrix (toArray(1)). Now that the dimensions match, in each pixel, solve the equation for f. Finally, convert the result from a 2D array image into a 1D array image (arrayProject()), then to a multi-band image (arrayFlatten()). The three bands correspond to the estimates of bare, vegetation and water fractions in f. Display the result where bare is red, vegetation is green, and water is blue (the addLayer() call expects bands in order, RGB)

var unmixImage = image.select(['B2', 'B3', 'B4', 'B5', 'B6', 'B7']);  
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(), vegetation, 30).values();  
var endmembers = ee.Array.cat([bareMean,  vegMean, waterMean], 1);  
var arrayImage = unmixImage.toArray().toArray(1);
var unmixed =  ee.Image(endmembers).matrixSolve(arrayImage);
var unmixedImage = unmixed.arrayProject([0])
				.arrayFlatten(
          [['bare', 'veg', 'water']]
        );
Map.centerObject(point, 11)
Map.addLayer(unmixedImage, {}, 'Unmixed');

Tassled Cap

Question 5: Repeat this process for Blacksburg, VA. Upload the mean spectra chart you generated for bare, water, and land and the map. Interpret the output of the image by selecting different pixels with Inspector

3B.1.4 - Hue-Saturation-Value Transform​

The Hue-Saturation-Value (HSV) model is a color transform of the RGB color space. Among many other things, it is useful for pan-sharpening. This involves converting an RGB to HSV, swapping the panchromatic band for the value (V), then converting back to RGB. For example, using the Landsat 8 scene:

var point = ee.Geometry.Point([-80.42, 37.22]);
var landsat = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA")
var image = ee.Image(landsat
  .filterBounds(point)
  .filterDate('2013-06-17', '2014-03-27')
  .sort('CLOUD_COVER')
  .first());
// Convert Landsat RGB bands to HSV   
var hsv = image.select(['B4', 'B3', 'B2']).rgbToHsv();
// Convert back to RGB,  swapping the image panchromatic band for the value.
var rgb = ee.Image.cat([
  hsv.select('hue'),
  hsv.select('saturation'),
  image.select(['B8'])]).hsvToRgb();
Map.centerObject(point, 12);
Map.addLayer(rgb, {max: 0.4}, 'Pan-sharpened');

Tassled Cap

Question 6: Compare the pan-sharpened image with the original image. What do you notice that's different? The same?

3B.2 - Spectral Transformation​

3B.2.1 - Linear Filtering​

In the present context, linear filtering (or convolution) refers to a linear combination of pixel values in a 'neighborhood', or kernel, where the weights of the kernel determine the coefficients in the linear combination (for this lab, the terms kernel and filter are interchangeable.) Filtering an image can be useful for extracting image information at different spatial frequencies by reducing noise. For this reason, smoothing filters are called low-pass filters (they let low-frequency data pass through) and edge detection filters are called high-pass filters. To implement filtering in Earth Engine use image.convolve() with an ee.Kernel for the argument.

3B.2.2 - Smoothing​

Smoothing means to convolve an image with a smoothing kernel.

A simple smoothing filter is a square kernel with uniform weights that sum to one. Convolving with this kernel sets each pixel to the mean of its neighborhood. Print a square kernel with uniform weights (this is sometimes called a "pillbox" or "boxcar" filter):

Expand the kernel object in the console to see the weights. This kernel is defined by how many pixels it covers (i.e. radius is in units of 'pixels'). A kernel with radius defined in 'meters' adjusts its size in pixels, so you can't visualize its weights, but it's more flexible in terms of adapting to inputs of different scale. In the following, use kernels with radius defined in meters except to visualize the weights.

var point = ee.Geometry.Point([-80.42, 37.22]);
var naip = ee.ImageCollection("USDA/NAIP/DOQQ")
var image = ee.Image(naip
  .filterBounds(point)
  .filterDate('2013-06-17', '2017-03-27')
  .first());
// Print a uniform kernel to see its weights.
print('A uniform kernel:', ee.Kernel.square(2));

Tassled Cap

Define a kernel with 2-meter radius (which corresponds to how many pixels in the NAIP image? Hint: try projection.nominalScale()), convolve the image with the kernel and compare the input image with the smoothed image:

// Define a square, uniform kernel.
var uniformKernel = ee.Kernel.square({
 radius: 2,
 units: 'meters',
});
// Filter the image by convolving with the smoothing filter.
var smoothed = image.convolve(uniformKernel);
var trueColorVis = {
  min: 0.0,
  max: 255.0,
};
Map.centerObject(point, 12);
Map.addLayer(smoothed, trueColorVis, 'smoothed image');

Tassled Cap

To make the image even smoother, try increasing the size of the neighborhood by increasing the pixel radius. to the human eye, the image is blurrier, but in many Machine Learning and Computer Vision algorithms, this process improves our output by reducing noise.

A Gaussian kernel can also be used for smoothing. Think of filtering with a Gaussian kernel as computing the weighted average in each pixel's neighborhood.

// Print a Gaussian kernel to see its weights.
print('A Gaussian kernel:', ee.Kernel.gaussian(2));
// Define a square Gaussian kernel:
var gaussianKernel = ee.Kernel.gaussian({
 radius: 2,
 units: 'meters',
});
// Filter the image by convolving with the Gaussian filter.
var gaussian = image.convolve(gaussianKernel);
var trueColorVis = {
  min: 0.0,
  max: 255.0,
};
Map.centerObject(point, 12);
Map.addLayer(gaussian, trueColorVis, 'smoothed image');

Question 7: What happens as you increase the pixel radius for each smoothing? What differences can you discern between the weights and the visualizations of the two smoothing kernels?

3B.2.3 - Edge Detection​

Convolving with an edge-detection kernel is used to find rapid changes in values that usually signify the edges of objects in the image data.

A classic edge detection kernel is the Laplacian kernel. Investigate the kernel weights and the image that results from convolving with the Laplacian. Other edge detection kernels include the Sobel, Prewitt and Roberts kernels. Learn more about additional edge detection methods in Earth Engine.

var point = ee.Geometry.Point([-80.42, 37.22]);
var naip = ee.ImageCollection("USDA/NAIP/DOQQ")
var image = ee.Image(naip
  .filterBounds(point)
  .filterDate('2013-06-17', '2017-03-27')
  .first());
// Define a Laplacian, or edge-detection kernel.
var laplace = ee.Kernel.laplacian8({ normalize: false });
// Apply the edge-detection kernel.
var edges = image.convolve(laplace);
var trueColorVis = {
  min: 0.0,
  max: 255.0,
  format: 'png'
};
Map.centerObject(point, 12);
Map.addLayer(edges, trueColorVis,'edges');

Tassled Cap

Question 8: Choose another edge detection method and upload the image - what characteristics do you see, and how does it compare to Laplacian Edge detection?

Additional Exercises​

Question 9: Look in google scholar to identify 2-3 publications that have used NDVI and two-three that used EVI. For what purposes were these indices used and what was the justification provided for that index?

Question 10: Discuss a spectral index that we did not cover in this lab relates to your area of research/interest. What is the the name of the spectral index, the formula used to calculate it, and what is it used to detect? Provide a citation of an academic article that has fruitfully used that index.

Question 11: Find 1-2 articles that use any of the linear transformation methods we practiced in this lab in the service of addressing an important social issue (e.g., one related to agriculture, environment, or development). Provide the citations and discussed how the transformation is used and how it's justified in the article.