--- title: "How Agricultural Fields Change in AlphaEarth Foundations" url: "https://wherobots.com/blog/alphaearth-foundations-agriculture-change-detection/" description: "Detect crop rotation and field-level change in AlphaEarth Foundations satellite embeddings using PCA, cosine distance, and Zarr mosaics with RasterFlow." published: "2026-04-07" modified: "2026-10-19" --- # How Agricultural Fields Change in AlphaEarth Foundations Author: Len Strnad [AlphaEarth Foundations](https://arxiv.org/abs/2507.22291) is a [geospatial AI model](https://wherobots.com/discover/geospatial-ai/) from Google DeepMind. It compresses a year of satellite observations into 64-dimensional embedding vectors for every 10-meter pixel on Earth’s land surfaces and shallow coastal waters annually. In this notebook, we use AEF to examine how agricultural fields evolve across several regions and ask a practical question: what kinds of field-level change become visible when we inspect the embeddings directly rather than relying only on imagery? Key takeaways - AlphaEarth Foundations (Google DeepMind) compresses a year of satellite observations into 64-dimensional embedding vectors for every 10-meter pixel on Earth's land surfaces and shallow coastal waters, annually. This notebook inspects field-level change in embedding space rather than relying only on imagery. - Three agricultural AOIs are used: Iowa, California, and Japan. RasterFlow build_gti_mosaic builds a Zarr mosaic per AOI from AlphaEarth Foundations COGs on Source Cooperative. - An example Iowa mosaic is 38 GB (9 time steps, 64 bands, 2017–2025). - Three views of temporal change: the first three embedding bands as RGB, a PCA projection of all 64 bands into three components, and an experimental cosine-distance-over-time workflow. PCA is described as making field-level seasonal structure easier to see than raw RGB bands. - The goal is not a single definitive interpretation of AEF. It is a practical workflow for exploring seasonal cycles, repeated interventions, and other field-scale dynamics in embedding space, using experimental RasterFlow PCA and distance-over-time tools. We begin by defining a few agricultural areas of interest and building Zarr mosaics over each one from COGs available in [Source Cooperative](https://source.coop/tge-labs/aef). From there, we compare several views of temporal change: the first three embedding bands rendered as RGB, a Principal Component Analysis (PCA) projection that makes broader structure easier to interpret, and an experimental distance-over-time workflow that highlights where embedding values remain stable or shift over time. The goal is not to claim a single definitive interpretation of AEF. Instead, this notebook shows a practical workflow for exploring seasonal cycles, repeated interventions, and other field-scale dynamics in embedding space. ## RasterFlow Tools Used RasterFlow is Wherobots’ serverless inference engine for [Earth Observation data](https://wherobots.com/discover/earth-observation-data/). It builds inference-ready mosaics from satellite imagery and runs distributed workflows at scale. This notebook uses three RasterFlow capabilities: - [build_gti_mosaic](https://docs.wherobots.com/reference/rasterflow/client#build_gti_mosaic) for constructing Zarr mosaics from AlphaEarth Foundations data - an experimental PCA workflow for dimensionality reduction - and an experimental distance-over-time workflow for measuring embedding change between timesteps. ## Selecting Agricultural AOIs Across Iowa, California, and Japan We start with three small AOIs that contain agricultural fields in different regions. This gives us a compact but varied set of examples for exploring how AlphaEarth Foundations embeddings change over time. ``` import geopandas as gpd from shapely.geometry import box from folium import Map ``` ``` iowa1 = [-94.732361, 41.911475, -94.237289, 42.112996] california1 = [-120.11095, 36.63082, -120.018082, 36.683839] japan1 = [141.018019, 39.226668, 141.067801, 39.249803] aois = [iowa1, california1, japan1] gdf = gpd.GeoDataFrame({"geometry": [box(*aoi) for aoi in aois]}, crs="EPSG:4326") # start on the first AOI m = Map(location=gdf.iloc[0].geometry.centroid.coords[0][::-1], zoom_start=8) gdf.explore(height="120%", width="120%", m=m) ``` Make this Notebook Trusted to load map: File -> Trust Notebook ``` # store this out to parameterize which Zarr mosaics we want to build for AEF gdf.to_parquet("s3://union-sandbox-unionai-wherobots/aois/field_change_aef.parquet") ``` ### AlphaEarth Foundations Zarr Mosaics by AOI The next step is to create an AlphaEarth Foundations embeddings Zarr mosaic for each AOI. This gives us a time-indexed raster store that we can query efficiently for visualization and downstream analysis. The output here is from the [build_gti_mosaic](https://docs.wherobots.com/reference/rasterflow/client#build_gti_mosaic) RasterFlow workflow. ``` import geopandas as gpd ``` ``` mosaic_index = gpd.read_parquet( "s3://union-sandbox-unionai-wherobots/ha/wherobots/wherobots-mosaics/development/rnmkrr4sbxsdf24t288p/d8cq2y2k3b1r4p9gr0k5n1ldh/1/7q/fdjxpi5q/61a6e8f637cce7ddf177ff7a6637aceb/" ) mosaic_index ``` | | geometry | location | | --- | --- | --- | | 0 | POLYGON ((-94.23729 41.91148, -94.23729 42.113… | s3://sandbox-wherobots-mosaics-tmp/mosaics/rnm… | | 1 | POLYGON ((-120.01808 36.63082, -120.01808 36.6… | s3://sandbox-wherobots-mosaics-tmp/mosaics/rnm… | | 2 | POLYGON ((141.0678 39.22667, 141.0678 39.2498,… | s3://sandbox-wherobots-mosaics-tmp/mosaics/rnm… | The result is one Zarr mosaic per AOI, which we can load independently for inspection and comparison. ### Inspect an Example AEF Zarr Mosaic ``` import xarray as xr ``` ``` ds = xr.open_zarr(mosaic_index.iloc[0]["location"]) ds ``` ``` Size: 38GB Dimensions: (time: 9, band: 64, y: 3020, x: 5512) Coordinates: * time (time) datetime64[ns] 72B 2017-01-01 2018-01-01 ... 2025-01-01 * band (band) object 512B 'band_0' 'band_1' ... 'band_62' 'band_63' * y (y) float64 24kB 5.178e+06 5.178e+06 ... 5.148e+06 5.148e+06 * x (x) float64 44kB -1.055e+07 -1.055e+07 ... -1.049e+07 spatial_ref int64 8B ... Data variables: variables (time, band, y, x) float32 38GB dask.array ``` xarray.Dataset - Dimensions: - time: 9 - band: 64 - y: 3020 - x: 5512 - Coordinates: (5) - time (time) datetime64\[ns\] 2017-01-01 … 2025-01-01 ``` array(['2017-01-01T00:00:00.000000000', '2018-01-01T00:00:00.000000000', '2019-01-01T00:00:00.000000000', '2020-01-01T00:00:00.000000000', '2021-01-01T00:00:00.000000000', '2022-01-01T00:00:00.000000000', '2023-01-01T00:00:00.000000000', '2024-01-01T00:00:00.000000000', '2025-01-01T00:00:00.000000000'], dtype='datetime64[ns]') ``` - band (band) object ‘band_0’ ‘band_1’ … ‘band_63’ ``` array(['band_0', 'band_1', 'band_2', 'band_3', 'band_4', 'band_5', 'band_6', 'band_7', 'band_8', 'band_9', 'band_10', 'band_11', 'band_12', 'band_13', 'band_14', 'band_15', 'band_16', 'band_17', 'band_18', 'band_19', 'band_20', 'band_21', 'band_22', 'band_23', 'band_24', 'band_25', 'band_26', 'band_27', 'band_28', 'band_29', 'band_30', 'band_31', 'band_32', 'band_33', 'band_34', 'band_35', 'band_36', 'band_37', 'band_38', 'band_39', 'band_40', 'band_41', 'band_42', 'band_43', 'band_44', 'band_45', 'band_46', 'band_47', 'band_48', 'band_49', 'band_50', 'band_51', 'band_52', 'band_53', 'band_54', 'band_55', 'band_56', 'band_57', 'band_58', 'band_59', 'band_60', 'band_61', 'band_62', 'band_63'], dtype=object) ``` - y (y) float64 5.178e+06 5.178e+06 … 5.148e+06 ``` array([5177915.753919, 5177905.753919, 5177895.753919, ..., 5147745.753919, 5147735.753919, 5147725.753919], shape=(3020,)) ``` - x (x) float64 -1.055e+07 … -1.049e+07 ``` array([-10545553.188165, -10545543.188165, -10545533.188165, ..., -10490463.188165, -10490453.188165, -10490443.188165], shape=(5512,)) ``` - spatial_ref () int64 … crs_wkt : PROJCS\[“WGS 84 / Pseudo-Mercator”,GEOGCS\[“WGS 84”,DATUM\[“WGS_1984”,SPHEROID\[“WGS 84”,6378137,298.257223563,AUTHORITY\[“EPSG”,”7030″\]\],AUTHORITY\[“EPSG”,”6326″\]\],PRIMEM\[“Greenwich”,0,AUTHORITY\[“EPSG”,”8901″\]\],UNIT\[“degree”,0.0174532925199433,AUTHORITY\[“EPSG”,”9122″\]\],AUTHORITY\[“EPSG”,”4326″\]\],PROJECTION\[“Mercator_1SP”\],PARAMETER\[“central_meridian”,0\],PARAMETER\[“scale_factor”,1\],PARAMETER\[“false_easting”,0\],PARAMETER\[“false_northing”,0\],UNIT\[“metre”,1,AUTHORITY\[“EPSG”,”9001″\]\],AXIS\[“Easting”,EAST\],AXIS\[“Northing”,NORTH\],EXTENSION\[“PROJ4″,”+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs”\],AUTHORITY\[“EPSG”,”3857″\]\] spatial_ref : PROJCS\[“WGS 84 / Pseudo-Mercator”,GEOGCS\[“WGS 84”,DATUM\[“WGS_1984”,SPHEROID\[“WGS 84”,6378137,298.257223563,AUTHORITY\[“EPSG”,”7030″\]\],AUTHORITY\[“EPSG”,”6326″\]\],PRIMEM\[“Greenwich”,0,AUTHORITY\[“EPSG”,”8901″\]\],UNIT\[“degree”,0.0174532925199433,AUTHORITY\[“EPSG”,”9122″\]\],AUTHORITY\[“EPSG”,”4326″\]\],PROJECTION\[“Mercator_1SP”\],PARAMETER\[“central_meridian”,0\],PARAMETER\[“scale_factor”,1\],PARAMETER\[“false_easting”,0\],PARAMETER\[“false_northing”,0\],UNIT\[“metre”,1,AUTHORITY\[“EPSG”,”9001″\]\],AXIS\[“Easting”,EAST\],AXIS\[“Northing”,NORTH\],EXTENSION\[“PROJ4″,”+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs”\],AUTHORITY\[“EPSG”,”3857″\]\] GeoTransform : -10545558.188164568 10.0 0.0 5177920.753919057 0.0 -10.0 ``` [1 values with dtype=int64] ``` - Data variables: (1) - variables (time, band, y, x) float32 dask.array\ | Array Chunk Bytes 35.72 GiB 4.00 MiB Shape (9, 64, 3020, 5512) (1, 1, 1024, 1024) Dask graph 10368 chunks in 2 graph layers Data type float32 numpy.ndarray | | | | --- | --- | --- | | | Array | Chunk | | Bytes | 35.72 GiB | 4.00 MiB | | Shape | (9, 64, 3020, 5512) | (1, 1, 1024, 1024) | | Dask graph | 10368 chunks in 2 graph layers | | | Data type | float32 numpy.ndarray | | ``` # pull this into memory once to reuse subset = ds["variables"][:, :3, :, :].compute() ``` ``` import imageio.v3 as iio import numpy as np def to_uint8_robust(array: xr.DataArray, low: xr.DataArray, high: xr.DataArray) -> np.ndarray: arr = array.transpose("y", "x", ...).astype("float32") scale = np.where((high - low) == 0, 1.0, (high - low)) arr = (arr - low) / scale arr = np.clip(arr, 0, 1) return (arr * 255).astype(np.uint8) def robust_params(da: xr.DataArray) -> tuple[xr.DataArray, xr.DataArray]: q_low, q_high = 1, 99 lo = da.quantile(q_low / 100, dim=("time", "y", "x")) hi = da.quantile(q_high / 100, dim=("time", "y", "x")) return lo, hi ``` ``` low, high = robust_params(subset) frames = [to_uint8_robust(subset.sel(time=t), low, high) for t in subset.time] iio.imwrite("assets/rgb_animation.gif", frames, duration=1000, loop=0) ``` To establish a baseline view, we render the first three embedding bands as an RGB animation over time. The color mapping is not directly semantic, but it provides a quick way to spot recurring seasonal patterns and abrupt changes. ![Animation of the first three Alpha Earth Embedding bands over time](https://wherobots.com/wp-content/uploads/2026/03/rgb_animation.gif?ver=1794530625) ### Using PCA to Identify Seasonal Patterns in Embeddings A more interpretable view comes from projecting the 64-dimensional embeddings into their first three principal components. PCA concentrates the dominant variance into three channels, which often makes field-level change easier to see. ``` # training sample for PCA using just the first year and a spatial slice sample = ds["variables"][:, :, :1024, :1024].compute() sample = sample.transpose("band","time","y","x").data.reshape(64, -1) sample = sample[:, ~np.isnan(sample).any(axis=0)] ``` ``` from sklearn.decomposition import PCA model = PCA(n_components=3) model.fit(sample.T) ``` ``` PCA(n_components=3) ``` **In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.** PCA [?Documentation for PCA](https://scikit-learn.org/1.8/modules/generated/sklearn.decomposition.PCA.html)iFitted Parameters | | [n_components n_components: int, float or ‘mle’, default=None Number of components to keep. if n_components is not set all components are kept:: n_components == min(n_samples, n_features) If “n_components == ‘mle’“ and “svd_solver == ‘full’“, Minka’s MLE is used to guess the dimension. Use of “n_components == ‘mle’“ will interpret “svd_solver == ‘auto’“ as “svd_solver == ‘full’“. If “0 \< n_components \< 1“ and “svd_solver == ‘full’“, select the number of components such that the amount of variance that needs to be explained is greater than the percentage specified by n_components. If “svd_solver == ‘arpack’“, the number of components must be strictly less than the minimum of n_features and n_samples. Hence, the None case results in:: n_components == min(n_samples, n_features) – 1](https://scikit-learn.org/1.8/modules/generated/sklearn.decomposition.PCA.html#:~:text=n_components,-int%2C%20float%20or%20%27mle%27%2C%20default%3DNone) | 3 | | --- | --- | --- | | | [copy copy: bool, default=True If False, data passed to fit are overwritten and running fit(X).transform(X) will not yield the expected results, use fit_transform(X) instead.](https://scikit-learn.org/1.8/modules/generated/sklearn.decomposition.PCA.html#:~:text=copy,-bool%2C%20default%3DTrue) | True | | | [whiten whiten: bool, default=False When True (False by default) the \`components\_\` vectors are multiplied by the square root of n_samples and then divided by the singular values to ensure uncorrelated outputs with unit component-wise variances. Whitening will remove some information from the transformed signal (the relative variance scales of the components) but can sometime improve the predictive accuracy of the downstream estimators by making their data respect some hard-wired assumptions.](https://scikit-learn.org/1.8/modules/generated/sklearn.decomposition.PCA.html#:~:text=whiten,-bool%2C%20default%3DFalse) | False | | | [svd_solver svd_solver: {‘auto’, ‘full’, ‘covariance_eigh’, ‘arpack’, ‘randomized’}, default=’auto’ “auto” : The solver is selected by a default ‘auto’ policy is based on \`X.shape\` and \`n_components\`: if the input data has fewer than 1000 features and more than 10 times as many samples, then the “covariance_eigh” solver is used. Otherwise, if the input data is larger than 500×500 and the number of components to extract is lower than 80% of the smallest dimension of the data, then the more efficient “randomized” method is selected. Otherwise the exact “full” SVD is computed and optionally truncated afterwards. “full” : Run exact full SVD calling the standard LAPACK solver via \`scipy.linalg.svd\` and select the components by postprocessing “covariance_eigh” : Precompute the covariance matrix (on centered data), run a classical eigenvalue decomposition on the covariance matrix typically using LAPACK and select the components by postprocessing. This solver is very efficient for n_samples >> n_features and small n_features. It is, however, not tractable otherwise for large n_features (large memory footprint required to materialize the covariance matrix). Also note that compared to the “full” solver, this solver effectively doubles the condition number and is therefore less numerical stable (e.g. on input data with a large range of singular values). “arpack” : Run SVD truncated to \`n_components\` calling ARPACK solver via \`scipy.sparse.linalg.svds\`. It requires strictly \`0 \< n_components \< min(X.shape)\` “randomized” : Run randomized SVD by the method of Halko et al. .. versionadded:: 0.18.0 .. versionchanged:: 1.5 Added the ‘covariance_eigh’ solver.](https://scikit-learn.org/1.8/modules/generated/sklearn.decomposition.PCA.html#:~:text=svd_solver,-%7B%27auto%27%2C%20%27full%27%2C%20%27covariance_eigh%27%2C%20%27arpack%27%2C%20%27randomized%27%7D%2C%20%20%20%20%20%20%20%20%20%20%20%20default%3D%27auto%27) | ‘auto’ | | | [tol tol: float, default=0.0 Tolerance for singular values computed by svd_solver == ‘arpack’. Must be of range \[0.0, infinity). .. versionadded:: 0.18.0](https://scikit-learn.org/1.8/modules/generated/sklearn.decomposition.PCA.html#:~:text=tol,-float%2C%20default%3D0.0) | 0.0 | | | [iterated_power iterated_power: int or ‘auto’, default=’auto’ Number of iterations for the power method computed by svd_solver == ‘randomized’. Must be of range \[0, infinity). .. versionadded:: 0.18.0](https://scikit-learn.org/1.8/modules/generated/sklearn.decomposition.PCA.html#:~:text=iterated_power,-int%20or%20%27auto%27%2C%20default%3D%27auto%27) | ‘auto’ | | | [n_oversamples n_oversamples: int, default=10 This parameter is only relevant when \`svd_solver=”randomized”\`. It corresponds to the additional number of random vectors to sample the range of \`X\` so as to ensure proper conditioning. See :func:\`~sklearn.utils.extmath.randomized_svd\` for more details. .. versionadded:: 1.1](https://scikit-learn.org/1.8/modules/generated/sklearn.decomposition.PCA.html#:~:text=n_oversamples,-int%2C%20default%3D10) | 10 | | | [power_iteration_normalizer power_iteration_normalizer: {‘auto’, ‘QR’, ‘LU’, ‘none’}, default=’auto’ Power iteration normalizer for randomized SVD solver. Not used by ARPACK. See :func:\`~sklearn.utils.extmath.randomized_svd\` for more details. .. versionadded:: 1.1](https://scikit-learn.org/1.8/modules/generated/sklearn.decomposition.PCA.html#:~:text=power_iteration_normalizer,-%7B%27auto%27%2C%20%27QR%27%2C%20%27LU%27%2C%20%27none%27%7D%2C%20default%3D%27auto%27) | ‘auto’ | | | [random_state random_state: int, RandomState instance or None, default=None Used when the ‘arpack’ or ‘randomized’ solvers are used. Pass an int for reproducible results across multiple function calls. See :term:\`Glossary \`. .. versionadded:: 0.18.0](https://scikit-learn.org/1.8/modules/generated/sklearn.decomposition.PCA.html#:~:text=random_state,-int%2C%20RandomState%20instance%20or%20None%2C%20default%3DNone) | None | To keep the notebook lightweight, we pull all time steps and bands for a small spatial subset rather than the full mosaic. AlphaEarth Foundations embeddings contain 64 bands, so this smaller window reduces data transfer while still preserving the temporal behavior we want to inspect. The same workflow can be scaled up for larger analyses. ``` # can take 30+ seconds subset = ds["variables"][:, :, :1024, :2048].compute() ``` ``` # Transform with PCA flattened = subset.transpose("band", "time", "y", "x").data.reshape(64, -1).T transformed = model.transform(flattened) img = transformed.T.reshape(3, subset.shape[0], subset.shape[2], subset.shape[3]).transpose( 1, 0, 2, 3 ) # Use the existing xarray data structure to build a nice home for PCA data. subset_copy = subset.copy() array = subset_copy.isel(band=slice(0, 3)) # insert the transformed PCA data array.data = img array["band"] = ["pca1", "pca2", "pca3"] # we chunk only to get a better view of the dataset array.to_dataset(name='variables') ``` ``` Size: 227MB Dimensions: (time: 9, x: 2048, y: 1024, band: 3) Coordinates: * time (time) datetime64[ns] 72B 2017-01-01 2018-01-01 ... 2025-01-01 * x (x) float64 16kB -1.055e+07 -1.055e+07 ... -1.053e+07 * y (y) float64 8kB 5.178e+06 5.178e+06 ... 5.168e+06 5.168e+06 * band (band) Size: 533MB Dimensions: (time: 8, band: 1, y: 3020, x: 5512) Coordinates: * time (time) datetime64[ns] 64B 2018-01-01 2019-01-01 ... 2025-01-01 * band (band) object 8B 'cosine_distance' * y (y) float64 24kB 5.178e+06 5.178e+06 ... 5.148e+06 5.148e+06 * x (x) float64 44kB -1.055e+07 -1.055e+07 ... -1.049e+07 spatial_ref int64 8B 0 Data variables: variables (time, band, y, x) float32 533MB 0.188 0.1655 ... 0.02265 ``` xarray.Dataset - Dimensions: - time: 8 - band: 1 - y: 3020 - x: 5512 - Coordinates: (5) - time (time) datetime64\[ns\] 2018-01-01 ... 2025-01-01 ``` array(['2018-01-01T00:00:00.000000000', '2019-01-01T00:00:00.000000000', '2020-01-01T00:00:00.000000000', '2021-01-01T00:00:00.000000000', '2022-01-01T00:00:00.000000000', '2023-01-01T00:00:00.000000000', '2024-01-01T00:00:00.000000000', '2025-01-01T00:00:00.000000000'], dtype='datetime64[ns]') ``` - band (band) object 'cosine_distance' ``` array(['cosine_distance'], dtype=object) ``` - y (y) float64 5.178e+06 5.178e+06 ... 5.148e+06 ``` array([5177915.753919, 5177905.753919, 5177895.753919, ..., 5147745.753919, 5147735.753919, 5147725.753919], shape=(3020,)) ``` - x (x) float64 -1.055e+07 ... -1.049e+07 ``` array([-10545553.188165, -10545543.188165, -10545533.188165, ..., -10490463.188165, -10490453.188165, -10490443.188165], shape=(5512,)) ``` - spatial_ref () int64 0 crs_wkt : PROJCS\["WGS 84 / Pseudo-Mercator",GEOGCS\["WGS 84",DATUM\["WGS_1984",SPHEROID\["WGS 84",6378137,298.257223563,AUTHORITY\["EPSG","7030"\]\],AUTHORITY\["EPSG","6326"\]\],PRIMEM\["Greenwich",0,AUTHORITY\["EPSG","8901"\]\],UNIT\["degree",0.0174532925199433,AUTHORITY\["EPSG","9122"\]\],AUTHORITY\["EPSG","4326"\]\],PROJECTION\["Mercator_1SP"\],PARAMETER\["central_meridian",0\],PARAMETER\["scale_factor",1\],PARAMETER\["false_easting",0\],PARAMETER\["false_northing",0\],UNIT\["metre",1,AUTHORITY\["EPSG","9001"\]\],AXIS\["Easting",EAST\],AXIS\["Northing",NORTH\],EXTENSION\["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"\],AUTHORITY\["EPSG","3857"\]\] spatial_ref : PROJCS\["WGS 84 / Pseudo-Mercator",GEOGCS\["WGS 84",DATUM\["WGS_1984",SPHEROID\["WGS 84",6378137,298.257223563,AUTHORITY\["EPSG","7030"\]\],AUTHORITY\["EPSG","6326"\]\],PRIMEM\["Greenwich",0,AUTHORITY\["EPSG","8901"\]\],UNIT\["degree",0.0174532925199433,AUTHORITY\["EPSG","9122"\]\],AUTHORITY\["EPSG","4326"\]\],PROJECTION\["Mercator_1SP"\],PARAMETER\["central_meridian",0\],PARAMETER\["scale_factor",1\],PARAMETER\["false_easting",0\],PARAMETER\["false_northing",0\],UNIT\["metre",1,AUTHORITY\["EPSG","9001"\]\],AXIS\["Easting",EAST\],AXIS\["Northing",NORTH\],EXTENSION\["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"\],AUTHORITY\["EPSG","3857"\]\] GeoTransform : -10545558.188164568 10.0 0.0 5177920.753919057 0.0 -10.0 ``` array(0) ``` - Data variables: (1) - variables (time, band, y, x) float32 0.188 0.1655 ... 0.01864 0.02265 ``` array([[[[0.18804508, 0.16548318, 0.1549539 , ..., 0.6547667 , 0.6600232 , 0.6827302 ], [0.09239805, 0.08352977, 0.07957667, ..., 0.65431345, 0.6571518 , 0.6796979 ], [0.07838935, 0.07207382, 0.07041121, ..., 0.64434373, 0.643643 , 0.67302835], ..., [0.0281769 , 0.02823567, 0.02864486, ..., 0.32737243, 0.3262853 , 0.30775493], [0.04068887, 0.04150879, 0.03881335, ..., 0.3269202 , 0.3267812 , 0.31013072], [0.04347068, 0.04405069, 0.04101866, ..., 0.3367141 , 0.33538705, 0.3183689 ]]], [[[0.20149559, 0.17581868, 0.16398937, ..., 0.4470787 , 0.448758 , 0.45887315], [0.09369123, 0.08462119, 0.07994598, ..., 0.44638008, 0.44314015, 0.45811075], [0.08294475, 0.07634479, 0.07309014, ..., 0.44814014, ... [0.04521638, 0.04491764, 0.04381305, ..., 0.46263814, 0.45893395, 0.4046142 ], [0.04455835, 0.04325521, 0.04292834, ..., 0.46341467, 0.4581054 , 0.40957177]]], [[[0.49459922, 0.46871775, 0.4611882 , ..., 0.1296547 , 0.12676072, 0.12663496], [0.4966467 , 0.4909042 , 0.4853648 , ..., 0.1275351 , 0.1253922 , 0.1257289 ], [0.4927944 , 0.4910075 , 0.48507315, ..., 0.12276459, 0.11952102, 0.12123233], ..., [0.08114195, 0.08070195, 0.08397031, ..., 0.01797432, 0.01881957, 0.02318203], [0.08387047, 0.0837459 , 0.08279037, ..., 0.01825058, 0.01814717, 0.02346212], [0.08379203, 0.08498359, 0.08619195, ..., 0.01759416, 0.01863509, 0.02264881]]]], shape=(8, 1, 3020, 5512), dtype=float32) ``` ``` array = ds["variables"][:, :, :, :] low = np.nanpercentile(array[0].values, 0) high = np.nanpercentile(array[0].values, 50) frames = [to_uint8_robust(array.sel(time=t)[0], low, high)[:, :] for t in array.time] iio.imwrite("assets/distance_over_time_50p.gif", frames, duration=1000, loop=0) ``` ![Animation showing distance over time in Alpha Earth Embedding space](https://wherobots.com/wp-content/uploads/2026/03/distance_over_time_50p.gif?ver=1794530625) In the animation above, darker pixels indicate embeddings that remain more similar over time, while lighter pixels indicate larger temporal differences. Interpretation is still evolving, but the visualization suggests several useful patterns. Persistent bright pixels may mark locations with repeated or abrupt change. Alternating light-dark cycles may indicate recurring agricultural activity such as planting, harvest, or irrigation. This does not replace domain knowledge, but it provides another lens for understanding how field-scale practices appear in AlphaEarth Foundations embeddings. ## What AlphaEarth Foundations Embeddings Show About Field-Level Dynamics AlphaEarth Foundations embeddings make it possible to examine agricultural change as a time series in embedding space rather than only as raw imagery. In this notebook, we built Zarr mosaics for several AOIs, visualized temporal variation in the first three embedding bands, projected the embeddings with PCA, and then compared that view with a distance-over-time workflow. Together, these helped draw insight into a few of the many ways AlphaEarth Foundations embeddings (and others) can provide insight into field-level dynamics. Check out [Change Detection Using AlphaEarth Foundations (Part 2)](https://wherobots.com/blog/alphaearth-foundations-change-detection-part-2/). Get Started with RasterFlow [Start Building](https://wherobots.com/rasterflow-preview/) ## Frequently Asked Questions ### What is AlphaEarth Foundations? A geospatial AI model from Google DeepMind. It compresses a year of satellite observations into 64-dimensional embedding vectors for every 10-meter pixel on Earth’s land surfaces and shallow coastal waters, produced annually. ### How do you visualize agricultural change in AEF embeddings? This notebook compares three views: the first three embedding bands rendered as RGB over 2017–2025, a PCA projection of the 64-d vectors into three components, and an experimental cosine-distance-over-time raster. PCA is presented as the more interpretable view for field-level seasonal patterns. ### Which regions are used in the AEF agriculture notebook? Three small AOIs with agricultural fields: one in Iowa, one in California, and one in Japan. Each gets its own Zarr mosaic from RasterFlow build_gti_mosaic so they can be inspected independently. ### Does this notebook claim to detect crop type from AlphaEarth? No. The author states the goal is not a definitive interpretation of AEF. The notebook shows a workflow for exploring seasonal cycles and field-scale dynamics in embedding space. Complementary posts cover zonal statistics on field polygons and LOO-medoid change scoring.