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 NIRred spectral space, Kauth and Thomas (1976) devised a rotational transform of the form
$p_1 = R^T p_0$
where:

$p_0$ is the original pixel vector (a stack of the p band values as an Array)

$p_1$ is the rotated pixel

R is an orthonormal basis of the new space (inverse of $R^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 GramSchmidt process to find the other basis vectors.
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 multiband image to an array image in which each pixel position stores an array. Do the matrix multiplication, then convert back to a multiband image.
 JavaScript
 Python
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('20080601', '20080901')
.sort('CLOUD_COVER')
.first());
var bands = ['B1', 'B2', 'B3', 'B4', 'B5', 'B7'];
// Make an Array Image, with a 1D Array per pixel.
var arrayImage1D = image.select(bands).toArray();
// Make an Array Image with a 2D 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 multiband image with TCnamed 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');
lat = 39.19; lon = 106.81;
zoom = 11
image_collection_name = "LANDSAT/LT05/C01/T1_TOA"
date_start = '20080601'
date_end = '20081201'
name = 'Snow'
point = ee.Geometry.Point([lon, lat])
image = (
ee.ImageCollection(image_collection_name)
.filterBounds(point)
.filterDate(date_start, date_end)
.sort('CLOUD_COVER')
.first()
)
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]
]);
bands = ['B1', 'B2', 'B3', 'B4', 'B5', 'B7']
# Make an Array Image, with a 1D Array per pixel.
arrayImage1D = image.select(bands).toArray();
# Make an Array Image with a 2D Array per pixel, 6x1.
arrayImage2D = arrayImage1D.toArray(1);
componentsImage = (ee.Image(coefficients).matrixMultiply(arrayImage2D)
# Get rid of the extra dimensions.
.arrayProject([0])
# Get a multiband image with TCnamed bands.
.arrayFlatten(
[['brightness', 'greenness', 'wetness', 'fourth', 'fifth', 'sixth']]
))
vizParams = {
'bands': ['brightness', 'greenness', 'wetness'],
'min': 0.1,
'max': [0.5, 0.1, 0.1]
}
map8 = build_map(lat, lon, zoom, vizParams, componentsImage, 'TC components')
map8
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 variancecovariance 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 multiband 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.
 JavaScript
 Python
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('20131117', '20140327')
.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 multiband 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');
lat = 37.22; lon = 80.42;
zoom = 11
image_collection_name = "LANDSAT/LC08/C01/T1_TOA"
date_start = '20131117'
date_end = '20140327'
name = 'PCA'
point = ee.Geometry.Point([lon, lat])
image = (
ee.ImageCollection(image_collection_name)
.filterBounds(point)
.filterDate(date_start, date_end)
.sort('CLOUD_COVER')
.first()
)
bands = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B10', 'B11']
arrayImage = image.select(bands).toArray();
covar = arrayImage.reduceRegion(
reducer = ee.Reducer.covariance(),
maxPixels = 1000000000
)
covarArray = ee.Array(covar.get('array'));
eigens = covarArray.eigen();
eigenVectors = eigens.slice(1, 1);
principalComponents = ee.Image(eigenVectors).matrixMultiply(arrayImage.toArray(1));
pcImage = (principalComponents
.arrayProject([0])
# Make the one band array image a multiband image, [] > image.
.arrayFlatten(
[['pc1', 'pc2', 'pc3', 'pc4', 'pc5', 'pc6', 'pc7', 'pc8']]
));
#Customize the visual parameters for PC1
vizParams = {
"opacity":1,
"bands": ["pc1"],
"min":420,
"max":400,
"gamma":1}
map8 = build_map(lat, lon, zoom, vizParams, pcImage.select('pc1'), 'TC components')
map8
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 = 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= 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.
 JavaScript
 Python
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('20130617', '20140327')
.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']);
lat = 48.11; lon = 123.25;
zoom = 11
image_collection_name = "LANDSAT/LC08/C01/T1_TOA"
date_start = '20130617'
date_end = '20140327'
name = 'PCA'
point = ee.Geometry.Point([lon, lat])
image = (
ee.ImageCollection(image_collection_name)
.filterBounds(point)
.filterDate(date_start, date_end)
.sort('CLOUD_COVER')
.first()
)
# Polygons of bare earth, water and vegetation
bare = (
ee.Geometry.Polygon(
[[[123.2370707334838, 48.1151452657945],
[123.2370707334838, 48.11351208612645],
[123.23410957473136, 48.11351208612645],
[123.23410957473136, 48.1151452657945]]]))
water = (
ee.Geometry.Polygon(
[[[123.2748188020549, 48.12059599002954],
[123.2748188020549, 48.118074835535865],
[123.2673086168132, 48.118074835535865],
[123.2673086168132, 48.12059599002954]]]))
vegetation = (
ee.Geometry.Polygon(
[[[123.27462568300582, 48.11533866992809],
[123.27462568300582, 48.114163936320416],
[123.27215805071212, 48.114163936320416],
[123.27215805071212, 48.11533866992809]]]))
unmixImage = image.select(['B2', 'B3', 'B4', 'B5', 'B6', 'B7'])
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:
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 1axis (or columns direction).
Turn the 6band 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 multiband 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)
 JavaScript
 Python
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');
bareMean = unmixImage.reduceRegion(
ee.Reducer.mean(), bare, 30).values();
waterMean = unmixImage.reduceRegion(
ee.Reducer.mean(), water, 30).values();
vegMean = unmixImage.reduceRegion(
ee.Reducer.mean(), vegetation, 30).values();
endmembers = ee.Array.cat([bareMean, vegMean, waterMean], 1);
arrayImage = unmixImage.toArray().toArray(1);
unmixed = ee.Image(endmembers).matrixSolve(arrayImage);
unmixedImage = unmixed.arrayProject([0]).arrayFlatten(
[['bare', 'veg', 'water']]
);
map8 = build_map(lat, lon, zoom, {}, unmixedImage, 'unmixedImage')
map8
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  HueSaturationValue Transform
The HueSaturationValue (HSV) model is a color transform of the RGB color space. Among many other things, it is useful for pansharpening. 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:
 JavaScript
 Python
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('20130617', '20140327')
.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}, 'Pansharpened');
lat = 37.22; lon = 80.42;
zoom = 11;
image_collection_name = "LANDSAT/LC08/C01/T1_TOA"
date_start = '20130617'
date_end = '20140327'
name = 'Hue Saturation'
point = ee.Geometry.Point([lon, lat])
image = (
ee.ImageCollection(image_collection_name)
.filterBounds(point)
.filterDate(date_start, date_end)
.sort('CLOUD_COVER')
.first()
)
# Convert Landsat RGB bands to HSV
hsv = image.select(['B4', 'B3', 'B2']).rgbToHsv();
# Convert back to RGB, swapping the image panchromatic band for the value.
rgb = ee.Image.cat([
hsv.select('hue'),
hsv.select('saturation'),
image.select(['B8'])]).hsvToRgb();
map8 = build_map(lat, lon, zoom, {'max': 0.4}, rgb, 'hue saturation')
map8
Question 6: Compare the pansharpened 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 lowpass filters (they let lowfrequency data pass through) and edge detection filters are called highpass 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.
 JavaScript
 Python
var point = ee.Geometry.Point([80.42, 37.22]);
var naip = ee.ImageCollection("USDA/NAIP/DOQQ")
var image = ee.Image(naip
.filterBounds(point)
.filterDate('20130617', '20170327')
.first());
// Print a uniform kernel to see its weights.
print('A uniform kernel:', ee.Kernel.square(2));
lat = 37.22; lon = 80.42;
zoom = 14
image_collection_name = "USDA/NAIP/DOQQ"
date_start = '20130617'
date_end = '20170327'
name = 'NAIP'
point = ee.Geometry.Point([lon, lat])
image = (
ee.ImageCollection(image_collection_name)
.filterBounds(point)
.filterDate(date_start, date_end)
.sort('CLOUD_COVER')
.first()
)
# Print a uniform kernel to see its weights.
# print('A uniform kernel:', ee.Kernel.square(2));
Define a kernel with 2meter 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:
 JavaScript
 Python
// 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');
# Define a square, uniform kernel.
uniformKernel = ee.Kernel.square(
radius = 2,
units = 'meters',
)
# Filter the image by convolving with the smoothing filter.
smoothed = image.convolve(uniformKernel)
vizParams = {
'min': 0.0,
'max': 255,
}
map9 = build_map(lat, lon, zoom, vizParams, smoothed, 'Smoothed')
map9
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.
 JavaScript
 Python
// 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');
# Print a Gaussian kernel to see its weights.
#print('A Gaussian kernel:', ee.Kernel.gaussian(2))
# Define a square Gaussian kernel:
gaussianKernel = ee.Kernel.gaussian(
radius = 2,
units = 'meters',
)
# Filter the image by convolving with the smoothing filter.
gaussiansmooth = image.convolve(gaussianKernel)
map9 = build_map(lat, lon, zoom, vizParams, gaussiansmooth, 'Smoothed')
map9
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 edgedetection 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.
 JavaScript
 Python
var point = ee.Geometry.Point([80.42, 37.22]);
var naip = ee.ImageCollection("USDA/NAIP/DOQQ")
var image = ee.Image(naip
.filterBounds(point)
.filterDate('20130617', '20170327')
.first());
// Define a Laplacian, or edgedetection kernel.
var laplace = ee.Kernel.laplacian8({ normalize: false });
// Apply the edgedetection 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');
# Define a Laplacian, or edgedetection kernel.
laplace = ee.Kernel.laplacian8(normalize= False)
edges = image.convolve(laplace)
vizParams = {
'min': 0,
'max': 255,
'format': 'png'
}
zoom = 18
map10 = build_map(lat, lon, zoom, vizParams, edges, 'Edge Detection')
map10
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 23 publications that have used NDVI and twothree 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 12 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.