diff --git a/samples/python/guides/arrays01.py b/samples/python/guides/arrays01.py new file mode 100644 index 000000000..ba011f90f --- /dev/null +++ b/samples/python/guides/arrays01.py @@ -0,0 +1,71 @@ +# Copyright 2024 The Google Earth Engine Community Authors +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# https://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + + +# Earth Engine Developer's Guide examples +# from 'Arrays - Arrays and Array Images' section + +# [START earthengine__arrays01__tc_array] +# Create an Array of Tasseled Cap coefficients. +coefficients = ee.Array([ + [0.3029, 0.2786, 0.4733, 0.5599, 0.508, 0.1872], + [-0.2941, -0.243, -0.5424, 0.7276, 0.0713, -0.1608], + [0.1511, 0.1973, 0.3283, 0.3407, -0.7117, -0.4559], + [-0.8239, 0.0849, 0.4396, -0.058, 0.2013, -0.2773], + [-0.3294, 0.0557, 0.1056, 0.1855, -0.4349, 0.8085], + [0.1079, -0.9023, 0.4119, 0.0575, -0.0259, 0.0252], +]) +# [END earthengine__arrays01__tc_array] + +# [START earthengine__arrays01__dimensions] +# Print the dimensions. +display(coefficients.length()) # [6,6] +# [END earthengine__arrays01__dimensions] + +# [START earthengine__arrays01__slice] +# Get the 1x6 greenness slice, display it. +greenness = coefficients.slice(axis=0, start=1, end=2, step=1) +display(greenness) +# [END earthengine__arrays01__slice] + +# [START earthengine__arrays01__array_image] +# Load a Landsat 8 image, select the bands of interest. +image = ee.Image('LANDSAT/LC08/C02/T1_TOA/LC08_044034_20140318').select( + ['B2', 'B3', 'B4', 'B5', 'B6', 'B7'] +) + +# Make an Array Image, with a 1-D Array per pixel. +array_image_1d = image.toArray() + +# Make an Array Image with a 2-D Array per pixel, 6x1. +array_image_2d = array_image_1d.toArray(1) +# [END earthengine__arrays01__array_image] + +# [START earthengine__arrays01__multiplication] +# Do a matrix multiplication: 1x6 times 6x1. +# Cast the greenness Array to an Image prior to multiplication. +greenness_array_image = ee.Image(greenness).matrixMultiply(array_image_2d) +# [END earthengine__arrays01__multiplication] + +# [START earthengine__arrays01__array_get] +# Get the result from the 1x1 array in each pixel of the 2-D array image. +greenness_image = greenness_array_image.arrayGet([0, 0]) + +# Display the input imagery with the greenness result. +m = geemap.Map() +m.set_center(-122.3, 37.562, 10) +m.add_layer(image, {'bands': ['B5', 'B4', 'B3'], 'min': 0, 'max': 0.5}, 'image') +m.add_layer(greenness_image, {'min': -0.1, 'max': 0.13}, 'greenness') +m +# [END earthengine__arrays01__array_get] diff --git a/samples/python/guides/arrays02.py b/samples/python/guides/arrays02.py new file mode 100644 index 000000000..4b00e909c --- /dev/null +++ b/samples/python/guides/arrays02.py @@ -0,0 +1,27 @@ +# Copyright 2024 The Google Earth Engine Community Authors +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# https://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + + +# Earth Engine Developer's Guide examples +# from 'Arrays - Arrays and Array Images' section + +# [START earthengine__arrays02__toy_arrays] +array_1d = ee.Array([1, 2, 3]) # [1,2,3] +array_2d = ee.Array.cat([array_1d], 1) # [[1],[2],[3]] +# [END earthengine__arrays02__toy_arrays] + +display( + 'Arrays:', + {'1-D array': array_1d.getInfo(), '2-D array': array_2d.getInfo()}, +) diff --git a/samples/python/guides/arrays03.py b/samples/python/guides/arrays03.py new file mode 100644 index 000000000..083fd5e90 --- /dev/null +++ b/samples/python/guides/arrays03.py @@ -0,0 +1,63 @@ +# Copyright 2024 The Google Earth Engine Community Authors +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# https://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + + +# Earth Engine Developer's Guide examples +# from 'Arrays - Arrays and Array Images' section + +# [START earthengine__arrays03__tc_transform] +# Define an Array of Tasseled Cap coefficients. +coefficients = ee.Array([ + [0.3029, 0.2786, 0.4733, 0.5599, 0.508, 0.1872], + [-0.2941, -0.243, -0.5424, 0.7276, 0.0713, -0.1608], + [0.1511, 0.1973, 0.3283, 0.3407, -0.7117, -0.4559], + [-0.8239, 0.0849, 0.4396, -0.058, 0.2013, -0.2773], + [-0.3294, 0.0557, 0.1056, 0.1855, -0.4349, 0.8085], + [0.1079, -0.9023, 0.4119, 0.0575, -0.0259, 0.0252], +]) + +# Load a Landsat 8 image, select the bands of interest. +image = ee.Image('LANDSAT/LC08/C02/T1_TOA/LC08_044034_20140318').select( + ['B2', 'B3', 'B4', 'B5', 'B6', 'B7'] +) + +# Make an Array Image, with a 1-D Array per pixel. +array_image_1d = image.toArray() + +# Make an Array Image with a 2-D Array per pixel, 6x1. +array_image_2d = array_image_1d.toArray(1) + +# Do a matrix multiplication: 6x6 times 6x1. +components_image = ( + ee.Image(coefficients) + .matrixMultiply(array_image_2d) + # Get rid of the extra dimensions. + .arrayProject([0]) + .arrayFlatten( + [['brightness', 'greenness', 'wetness', 'fourth', 'fifth', 'sixth']] + ) +) + +# Display the first three bands of the result and the input imagery. +viz_params = { + 'bands': ['brightness', 'greenness', 'wetness'], + 'min': -0.1, + 'max': [0.5, 0.1, 0.1], +} +m = geemap.Map() +m.set_center(-122.3, 37.562, 10) +m.add_layer(image, {'bands': ['B5', 'B4', 'B3'], 'min': 0, 'max': 0.5}, 'image') +m.add_layer(components_image, viz_params, 'components') +m +# [END earthengine__arrays03__tc_transform] diff --git a/samples/python/guides/arrays04.py b/samples/python/guides/arrays04.py new file mode 100644 index 000000000..ff2655525 --- /dev/null +++ b/samples/python/guides/arrays04.py @@ -0,0 +1,187 @@ +# Copyright 2024 The Google Earth Engine Community Authors +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# https://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +# Arrays +# Earth Engine Developer's Guide examples +# Array transformation section + +# [START earthengine__arrays04__harmonic_model] +import math + + +# Scales and masks Landsat 8 surface reflectance images. +def prep_sr_l8(image): + # Develop masks for unwanted pixels (fill, cloud, cloud shadow). + qa_mask = image.select('QA_PIXEL').bitwiseAnd(int('11111', 2)).eq(0) + saturation_mask = image.select('QA_RADSAT').eq(0) + + # Apply the scaling factors to the appropriate bands. + optical_bands = image.select('SR_B.').multiply(0.0000275).add(-0.2) + thermal_bands = image.select('ST_B.*').multiply(0.00341802).add(149.0) + + # Replace the original bands with the scaled ones and apply the masks. + return ( + image.addBands(optical_bands, None, True) + .addBands(thermal_bands, None, True) + .updateMask(qa_mask) + .updateMask(saturation_mask) + ) + + +# Load a Landsat 8 surface reflectance image collection. +collection = ( + ee.ImageCollection('LANDSAT/LC08/C02/T1_L2') + # Filter to get only two years of data. + .filterDate('2019-04-01', '2021-04-01') + # Filter to get only imagery at a point of interest. + .filterBounds(ee.Geometry.Point(-122.08709, 36.9732)) + # Prepare images by mapping the prep_sr_l8 function over the collection. + .map(prep_sr_l8) + # Select NIR and red bands only. + .select(['SR_B5', 'SR_B4']) + # Sort the collection in chronological order. + .sort('system:time_start', True) +) + + +# This function computes the predictors and the response from the input. +def make_variables(image): + # Compute time of the image in fractional years relative to the Epoch. + year = ee.Image(image.date().difference(ee.Date('1970-01-01'), 'year')) + # Compute the season in radians, one cycle per year. + season = year.multiply(2 * math.pi) + # Return an image of the predictors followed by the response. + return ( + image.select() + .addBands(ee.Image(1)) # 0. constant + .addBands(year.rename('t')) # 1. linear trend + .addBands(season.sin().rename('sin')) # 2. seasonal + .addBands(season.cos().rename('cos')) # 3. seasonal + .addBands(image.normalizedDifference().rename('NDVI')) # 4. response + .toFloat() + ) + + +# Define the axes of variation in the collection array. +image_axis = 0 +band_axis = 1 + +# Convert the collection to an array. +array = collection.map(make_variables).toArray() + +# Check the length of the image axis (number of images). +array_length = array.arrayLength(image_axis) +# Update the mask to ensure that the number of images is greater than or +# equal to the number of predictors (the linear model is solvable). +array = array.updateMask(array_length.gt(4)) + +# Get slices of the array according to positions along the band axis. +predictors = array.arraySlice(band_axis, 0, 4) +response = array.arraySlice(band_axis, 4) +# [END earthengine__arrays04__harmonic_model] + +# Solve for linear regression coefficients in three different ways. +# All three methods produce equivalent results, but some are easier. +# [START earthengine__arrays04__hard_way] +# Compute coefficients the hard way. +coefficients_1 = ( + predictors.arrayTranspose() + .matrixMultiply(predictors) + .matrixInverse() + .matrixMultiply(predictors.arrayTranspose()) + .matrixMultiply(response) +) +# [END earthengine__arrays04__hard_way] + +# [START earthengine__arrays04__easy_way] +# Compute coefficients the easy way. +coefficients_2 = predictors.matrixPseudoInverse().matrixMultiply(response) +# [END earthengine__arrays04__easy_way] + +# [START earthengine__arrays04__easiest_way] +# Compute coefficients the easiest way. +coefficients_3 = predictors.matrixSolve(response) +# [END earthengine__arrays04__easiest_way] + +# [START earthengine__arrays04__image_flatten] +# Turn the results into a multi-band image. +coefficients_image = ( + coefficients_3 + # Get rid of the extra dimensions. + .arrayProject([0]).arrayFlatten([['constant', 'trend', 'sin', 'cos']]) +) +# [END earthengine__arrays04__image_flatten] + +import pandas as pd +import matplotlib.pyplot as plt + +# Use this mask for cartographic purposes, to get rid of water areas. +mask = ee.Image('CGIAR/SRTM90_V4').mask() + +# Display the result. +m = geemap.Map() +m.set_center(-122.08709, 36.9732, 13) +m.add_layer( + coefficients_image.updateMask(mask), + { + 'bands': ['sin', 'trend', 'cos'], + 'min': [-0.05, -0.1, -0.05], + 'max': [0.05, 0.1, 0.05], + }, +) +display(m) + + +# Map a function over the collection to compute the fitted values at each time. +def add_fitted_trend(image): + coeffs = coefficients_image.select(['constant', 'trend', 'sin', 'cos']) + predicted = ( + image.select(['constant', 't', 'sin', 'cos']) + .multiply(coeffs) + .reduce('sum') + .rename('fitted') + ) + + return image.select('NDVI').addBands(predicted) + + +fitted = collection.map(make_variables).map(add_fitted_trend) + +# Define a point that expresses a seasonal NDVI pattern. +roi = ee.Geometry.Point(-122.08709, 36.9732).buffer(120) + + +def extract_chart_data(image): + date = image.date().format('YYYY-MM-dd') + values = image.reduceRegion(reducer=ee.Reducer.mean(), geometry=roi, scale=30) + + return ee.Feature(roi, values).set('date', date) + + +fitted = fitted.map(extract_chart_data) + +df = ee.data.computeFeatures( + {'expression': fitted, 'fileFormat': 'PANDAS_DATAFRAME'} +) + +plt.figure(figsize=(10, 6)) +plt.plot(df['date'], df['NDVI'], label='NDVI') +plt.plot(df['date'], df['fitted'], label='Fitted') +plt.xticks(rotation=90) +plt.xlabel('Date') +plt.ylabel('Value') +plt.legend() +plt.title('NDVI and Fitted Trend') +plt.grid(True) +plt.show() diff --git a/samples/python/guides/arrays05.py b/samples/python/guides/arrays05.py new file mode 100644 index 000000000..7c7efb352 --- /dev/null +++ b/samples/python/guides/arrays05.py @@ -0,0 +1,119 @@ +# Copyright 2024 The Google Earth Engine Community Authors +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# https://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + + +# Earth Engine Developer's Guide examples +# from 'Arrays - Eigen analysis' section + +# Use these bands. +band_names = ee.List(['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B10', 'B11']) + +# Load a landsat 8 image and select the bands of interest. +image = ee.Image('LANDSAT/LC08/C02/T1/LC08_044034_20140318').select(band_names) + +# Display the input imagery and the region in which to do the PCA. +region = image.geometry() +m = geemap.Map() +m.center_object(region, 10) +m.add_layer(ee.Image().paint(region, 0, 2), {}, 'Region') +m.add_layer( + image, + {'bands': ['B5', 'B4', 'B2'], 'min': 0, 'max': 20000}, + 'Original Image', +) +display(m) + +# Set an appropriate scale for Landsat data. +scale = 30 + +# Mean center the data to enable a faster covariance reducer +# and an SD stretch of the principal components. +mean_dict = image.reduceRegion( + reducer=ee.Reducer.mean(), geometry=region, scale=scale, maxPixels=1e9 +) +means = mean_dict.toImage(band_names) +centered = image.subtract(means) + + +# This helper function returns a list of new band names. +def get_new_band_names(prefix): + seq = ee.List.sequence(1, band_names.length()) + + def add_prefix_and_number(b): + return ee.String(prefix).cat(ee.Number(b).int()) + + return seq.map(add_prefix_and_number) + + +# This function accepts mean centered imagery, a scale and +# a region in which to perform the analysis. It returns the +# Principal Components (PC) in the region as a new image. +# [START earthengine__arrays05__principal_components] +def get_principal_components(centered, scale, region): + # Collapse bands into 1D array + arrays = centered.toArray() + + # Compute the covariance of the bands within the region. + covar = arrays.reduceRegion( + reducer=ee.Reducer.centeredCovariance(), + geometry=region, + scale=scale, + maxPixels=1e9, + ) + + # Get the 'array' covariance result and cast to an array. + # This represents the band-to-band covariance within the region. + covar_array = ee.Array(covar.get('array')) + + # Perform an eigen analysis and slice apart the values and vectors. + eigens = covar_array.eigen() + + # This is a P-length vector of Eigenvalues. + eigen_values = eigens.slice(1, 0, 1) + # This is a PxP matrix with eigenvectors in rows. + eigen_vectors = eigens.slice(1, 1) + + # Convert the array image to 2D arrays for matrix computations. + array_image = arrays.toArray(1) + + # Left multiply the image array by the matrix of eigenvectors. + principal_components = ee.Image(eigen_vectors).matrixMultiply(array_image) + + # Turn the square roots of the Eigenvalues into a P-band image. + sd_image = ( + ee.Image(eigen_values.sqrt()) + .arrayProject([0]) + .arrayFlatten([get_new_band_names('sd')]) + ) + + # Turn the PCs into a P-band image, normalized by SD. + return ( + # Throw out an an unneeded dimension, [[]] -> []. + principal_components.arrayProject([0]) + # Make the one band array image a multi-band image, [] -> image. + .arrayFlatten([get_new_band_names('pc')]) + # Normalize the PCs by their SDs. + .divide(sd_image) + ) + + +# [END earthengine__arrays05__principal_components] + +# Get the PCs at the specified scale and in the specified region +pc_image = get_principal_components(centered, scale, region) + +# Plot each PC as a new layer +for i in range(len(band_names.getInfo())): + band = pc_image.bandNames().get(i).getInfo() + m.add_layer(pc_image.select([band]), {'min': -2, 'max': 2}, band) diff --git a/samples/python/guides/arrays06.py b/samples/python/guides/arrays06.py new file mode 100644 index 000000000..210fb7ebf --- /dev/null +++ b/samples/python/guides/arrays06.py @@ -0,0 +1,92 @@ +# Copyright 2024 The Google Earth Engine Community Authors +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# https://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + + +# Earth Engine Developer's Guide examples +# from 'Arrays - Sorting and Reducing' section + + +# [START earthengine__arrays06__sort_reduce] +# Define a function that scales and masks Landsat 8 surface reflectance images +# and adds an NDVI band. +def prep_sr_l8(image): + # Develop masks for unwanted pixels (fill, cloud, cloud shadow). + qa_mask = image.select('QA_PIXEL').bitwiseAnd(int('11111', 2)).eq(0) + saturation_mask = image.select('QA_RADSAT').eq(0) + + # Apply the scaling factors to the appropriate bands. + optical_bands = image.select('SR_B.').multiply(0.0000275).add(-0.2) + thermal_bands = image.select('ST_B.*').multiply(0.00341802).add(149.0) + + # Calculate NDVI. + ndvi = optical_bands.normalizedDifference(['SR_B5', 'SR_B4']).rename('NDVI') + + # Replace the original bands with the scaled ones and apply the masks. + return ( + image.addBands(optical_bands, None, True) + .addBands(thermal_bands, None, True) + .addBands(ndvi) + .updateMask(qa_mask) + .updateMask(saturation_mask) + ) + + +# Define an arbitrary region of interest as a point. +roi = ee.Geometry.Point(-122.26032, 37.87187) + +# Load a Landsat 8 surface reflectance collection. +collection = ( + ee.ImageCollection('LANDSAT/LC08/C02/T1_L2') + # Filter to get only imagery at a point of interest. + .filterBounds(roi) + # Filter to get only six months of data. + .filterDate('2021-01-01', '2021-07-01') + # Prepare images by mapping the prep_sr_l8 function over the collection. + .map(prep_sr_l8) + # Select the bands of interest to avoid taking up unneeded memory. + .select('SR_B.|NDVI') +) + +# Convert the collection to an array. +array = collection.toArray() + +# Label of the axes. +image_axis = 0 +band_axis = 1 + +# Get the NDVI slice and the bands of interest. +band_names = collection.first().bandNames() +bands = array.arraySlice(band_axis, 0, band_names.length()) +ndvi = array.arraySlice(band_axis, -1) + +# Sort by descending NDVI. +sorted = bands.arraySort(ndvi.multiply(-1)) + +# Get the highest 20% NDVI observations per pixel. +num_images = sorted.arrayLength(image_axis).multiply(0.2).int() +highest_ndvi = sorted.arraySlice(image_axis, 0, num_images) + +# Get the mean of the highest 20% NDVI observations by reducing +# along the image axis. +mean = highest_ndvi.arrayReduce(reducer=ee.Reducer.mean(), axes=[image_axis]) + +# Turn the reduced array image into a multi-band image for display. +mean_image = mean.arrayProject([band_axis]).arrayFlatten([band_names]) +m = geemap.Map() +m.center_object(roi, 12) +m.add_layer( + mean_image, {'bands': ['SR_B6', 'SR_B5', 'SR_B4'], 'min': 0, 'max': 0.4} +) +m +# [END earthengine__arrays06__sort_reduce]