Skip to content

Commit

Permalink
Convert API code snippets from JS to Py
Browse files Browse the repository at this point in the history
PiperOrigin-RevId: 628993393
  • Loading branch information
jdbcode authored and copybara-github committed Apr 29, 2024
1 parent 450879e commit 7baab22
Show file tree
Hide file tree
Showing 6 changed files with 559 additions and 0 deletions.
71 changes: 71 additions & 0 deletions samples/python/guides/arrays01.py
Original file line number Diff line number Diff line change
@@ -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]
27 changes: 27 additions & 0 deletions samples/python/guides/arrays02.py
Original file line number Diff line number Diff line change
@@ -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()},
)
63 changes: 63 additions & 0 deletions samples/python/guides/arrays03.py
Original file line number Diff line number Diff line change
@@ -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]
187 changes: 187 additions & 0 deletions samples/python/guides/arrays04.py
Original file line number Diff line number Diff line change
@@ -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()
Loading

0 comments on commit 7baab22

Please sign in to comment.