ventuari.dev

Tutoriales / Vegetación

Serie temporal de NDVI con Sentinel-2 en Google Earth Engine: fenología de cultivos

Cómo calcular el NDVI mes a mes con Sentinel-2 en Google Earth Engine, convertirlo en una serie temporal por píxel con array images y obtener la amplitud anual, el mes de máximo verdor y la duración de la temporada de cada parcela. El mismo flujo en el Code Editor (JavaScript) y en Colab (Python).

Por VentuariPublicado el 11 min de lectura

NDVI de la zona agrícola al sur de Neiva, Huila: parcelas en rojo con NDVI bajo entre vegetación en tonos de verde.

Qué es el NDVI y cómo se calcula

La vegetación sana absorbe casi toda la luz roja para la fotosíntesis y refleja con fuerza el infrarrojo cercano. El NDVI (índice de vegetación de diferencia normalizada) aprovecha ese contraste:

NDVI=NIR−RojoNIR+Rojo\mathrm{NDVI} = \dfrac{\mathrm{NIR} - \text{Rojo}}{\mathrm{NIR} + \text{Rojo}}

En Sentinel-2 el rojo es la banda B4 (664,5 nm) y el infrarrojo cercano es B8 (835,1 nm), ambas a 10 m de resolución, según el catálogo de Earth Engine para Sentinel-2A. El resultado va de −1 a 1: cuanto más alto, más vegetación verde y activa. Para leer esos valores sin umbrales engañosos, mira qué es el NDVI y cómo interpretar sus valores.

Un solo NDVI es una foto: dice cuánto verdor hay en una fecha, pero no cómo cambia. Un cultivo de ciclo corto y un pastizal pueden tener el mismo NDVI en mayo y comportarse de forma muy distinta el resto del año. Con un NDVI por mes, cada píxel tiene su curva anual, y de esa curva salen la fenología y la duración de la temporada.

Un compuesto sin nubes por mes con Cloud Score+

Las nubes y sus sombras hunden el NDVI y crean falsos mínimos en la curva, así que hay que enmascararlas antes de nada. Cloud Score+ es un producto de Google disponible en Earth Engine que asigna a cada píxel una puntuación de 0 a 1 según qué tan despejado está (1 es completamente despejado); filtrar con un umbral de 0.6 suele dejar pasar menos bordes de nubes y sombras que las bandas QA60 o SCL.

El script arma la mediana de las imágenes despejadas de cada mes calendario. Si un mes no tiene ninguna imagen, devuelve una imagen vacía y enmascarada en lugar de fallar, y la capa Valid months muestra cuántos meses con dato tiene cada píxel.

Array images: de 12 bandas a un vector por píxel

Los doce NDVI se apilan en una imagen de doce bandas (ndvi_01 a ndvi_12). Con toArray(), cada píxel pasa a guardar un solo vector con sus doce valores:

[ NDVI01, NDVI02, …, NDVI12 ]\left[\, \mathrm{NDVI}_{01},\ \mathrm{NDVI}_{02},\ \ldots,\ \mathrm{NDVI}_{12} \,\right]

Sobre ese vector se trabaja con operaciones de array, que recorren los doce meses de cada píxel a la vez:

  • arrayReduce(ee.Reducer.max(), [0]) y min() dan el NDVI máximo y mínimo del año.
  • arrayArgmax() da la posición del máximo; sumando 1, el mes de máximo verdor.
  • arrayGet([0]) convierte el resultado, que sigue siendo un array de un elemento, en una imagen normal que se puede restar y visualizar.

Hay un detalle que decide si el resultado es correcto: toArray() enmascara el píxel entero si falta un solo mes. Por eso el script rellena los meses vacíos con −2 para calcular el máximo y con 2 para el mínimo. Como el NDVI va de −1 a 1, esos valores nunca pueden ganar, y un mes nublado no altera el resultado.

Amplitud y umbral fenológico

La amplitud es la diferencia entre el máximo y el mínimo anual. Un bosque o una zona urbana cambian poco a lo largo del año; un cultivo pasa de suelo casi desnudo a vegetación densa y vuelve a bajar. El script solo analiza los píxeles con al menos 8 meses válidos y una amplitud de 0.15 o más, y deja fuera el resto.

Para cada píxel calcula un umbral relativo a su propia curva:

umbral=NDVImin⁡+0.25×amplitud\text{umbral} = \mathrm{NDVI}_{\min} + 0.25 \times \text{amplitud}

Con un mínimo de 0.20 y un máximo de 0.80, el umbral queda en 0.35. Los meses por encima del umbral cuentan como meses de vegetación activa: el primero es el inicio, el último es el fin y la duración es fin − inicio + 1. Un umbral relativo funciona igual en una parcela de arroz y en un pastizal seco, porque cada píxel se compara con su propio rango y no con un valor fijo.

Inicio y fin describen la fenología espectral: cuándo el verdor sube y baja. No equivalen a las fechas de siembra y cosecha: la siembra suele ocurrir semanas antes de que el verdor suba, y la cosecha, hacia el final de la bajada o después.

El script

Pega esto en el Code Editor y cambia el año si quieres analizar otro. El área de análisis es un polígono dibujado en el mapa: usa la herramienta de dibujo para trazar tu zona agrícola; el polígono aparece en Imports con el nombre geometry, que es el que usa el script.

ndvi_fenologia.js
/**
 * Google Earth Engine - Array Image example
 *
 * Annual crop phenology using monthly Sentinel-2 NDVI.
 *
 * Each pixel is represented as a temporal vector:
 *
 *     [NDVI_Jan, NDVI_Feb, ..., NDVI_Dec]
 *
 * The Array Image is used to derive:
 *
 *     - Minimum annual NDVI
 *     - Maximum annual NDVI
 *     - Month of maximum NDVI
 *     - Annual NDVI amplitude
 *
 * A relative phenological threshold is then used to estimate:
 *
 *     - Start month of vegetation activity
 *     - End month of vegetation activity
 *     - Approximate duration of the active season
 *
 * Important:
 *     Start and end months represent spectral phenology.
 *     They should not be interpreted directly as exact sowing
 *     or harvest dates.
 */


/**
 * Configuration
 */
var year = 2025;

var clear_threshold = 0.60;

var phenology_fraction = 0.25;

var min_amplitude = 0.15;

var min_valid_months = 8;


/**
 * Area of interest
 */
var aoi = geometry;


/**
 * Collections
 */
var s2 = ee.ImageCollection(
    'COPERNICUS/S2_SR_HARMONIZED'
);

var cs_plus = ee.ImageCollection(
    'GOOGLE/CLOUD_SCORE_PLUS/V1/S2_HARMONIZED'
);


/**
 * Months
 */
var months = [
    1, 2, 3, 4, 5, 6,
    7, 8, 9, 10, 11, 12,
];

var month_names = [
    'January',
    'February',
    'March',
    'April',
    'May',
    'June',
    'July',
    'August',
    'September',
    'October',
    'November',
    'December',
];


/**
 * Annual temporal window
 */
var year_start = ee.Date.fromYMD(
    year,
    1,
    1
);

var year_end = year_start.advance(1, 'year');

print(
    'Temporal window',
    year_start.format('YYYY-MM-dd'),
    year_end.advance(-1, 'day').format('YYYY-MM-dd')
);


/**
 * Returns a cloud-filtered Sentinel-2 collection.
 *
 * @param {ee.Date} start_date Start date.
 * @param {ee.Date} end_date End date.
 * @return {ee.ImageCollection} Cloud-filtered collection.
 */
var get_clear_collection = function (start_date, end_date) {
    return s2
        .filterBounds(aoi)
        .filterDate(start_date, end_date)
        .linkCollection(
            cs_plus,
            ['cs_cdf']
        )
        .map(function (image) {
            var clear_mask = image
                .select('cs_cdf')
                .gte(clear_threshold);

            return image
                .updateMask(clear_mask);
        });
};


/**
 * Creates a monthly NDVI composite.
 *
 * Each image represents exactly one calendar month.
 *
 * @param {number} month Month number from 1 to 12.
 * @return {ee.Image} Monthly NDVI image.
 */
var get_monthly_ndvi = function (month) {
  
    var start_date = ee.Date.fromYMD(year, month, 1);

    var end_date = start_date.advance(1, 'month');

    var collection = get_clear_collection(start_date, end_date);

    var empty_image = ee.Image
        .constant([ 0, 0 ])
        .rename([ 'B8', 'B4' ])
        .updateMask(
            ee.Image.constant(0)
        );

    var composite = ee.Image(
        ee.Algorithms.If(
            collection.size().gt(0),
            collection
                .select(['B8', 'B4' ])
                .median(),
            empty_image
        )
    );

    var ndvi = composite
        .normalizedDifference(['B8', 'B4'])
        .rename(
            'ndvi_' +
            String(month).padStart(2, '0')
        )
        .clip(aoi);

    return ndvi;
};


/**
 * Monthly NDVI images
 */
var monthly_images = months.map(function (month) {
    return get_monthly_ndvi(month);
});


/**
 * Monthly image availability
 */
var monthly_summary = ee.FeatureCollection(
    months.map(function (month) {
        var start_date = ee.Date.fromYMD( year, month, 1 );

        var end_date = start_date.advance( 1, 'month' );

        var image_count = get_clear_collection(start_date, end_date).size();

        return ee.Feature(
            null,
            {
                month: month,
                start_date: start_date.format('YYYY-MM-dd'),
                end_date: end_date
                    .advance(-1, 'day')
                    .format('YYYY-MM-dd'),
                image_count: image_count,
            }
        );
    })
);


/**
 * Multiband NDVI image
 *
 * Conventional representation:
 *
 *     ndvi_01
 *     ndvi_02
 *     ...
 *     ndvi_12
 */
var ndvi_stack = ee.Image.cat(
    monthly_images
);


/**
 * Number of valid months
 *
 * Indicates how many monthly NDVI values are available
 * for each pixel.
 */
var valid_month_count = ndvi_stack
    .mask()
    .reduce(ee.Reducer.sum())
    .rename('valid_month_count')
    .clip(aoi);


/**
 * Array Image
 *
 * Each pixel is converted from:
 *
 *     ndvi_01
 *     ndvi_02
 *     ...
 *     ndvi_12
 *
 * into:
 *
 *     [NDVI_01, NDVI_02, ..., NDVI_12]
 */
var ndvi_array = ndvi_stack.toArray();


/**
 * Arrays used for robust temporal calculations
 *
 * NDVI normally ranges from -1 to 1.
 *
 * Missing values are temporarily replaced with:
 *
 *     -2 for maximum calculations
 *      2 for minimum calculations
 *
 * Missing observations therefore cannot become
 * the annual maximum or minimum.
 */
var ndvi_array_max = ndvi_stack
    .unmask(-2)
    .toArray();

var ndvi_array_min = ndvi_stack
    .unmask(2)
    .toArray();


/**
 * Maximum annual NDVI
 */
var max_ndvi = ndvi_array_max
    .arrayReduce(
        ee.Reducer.max(),
        [0]
    )
    .arrayGet([0])
    .rename('max_ndvi');


/**
 * Minimum annual NDVI
 */
var min_ndvi = ndvi_array_min
    .arrayReduce(
        ee.Reducer.min(),
        [0]
    )
    .arrayGet([0])
    .rename('min_ndvi');


/**
 * Annual NDVI amplitude
 *
 * Large values indicate stronger seasonal change.
 */
var ndvi_amplitude = max_ndvi
    .subtract(min_ndvi)
    .rename('ndvi_amplitude');


/**
 * Seasonal validity mask
 *
 * A pixel must:
 *
 *     - Have enough valid months.
 *     - Show sufficient annual NDVI variation.
 */
var seasonal_mask = valid_month_count
    .gte(min_valid_months)
    .and(
        ndvi_amplitude.gte(min_amplitude)
    );


/**
 * Position of maximum annual NDVI
 *
 * Array positions:
 *
 *     0 = January
 *     1 = February
 *     ...
 *     11 = December
 */
var max_position = ndvi_array_max
    .arrayArgmax()
    .arrayGet([0])
    .rename('max_position');


/**
 * Month of maximum NDVI
 */
var peak_month = max_position
    .add(1)
    .rename('peak_month')
    .toByte()
    .updateMask(seasonal_mask);


/**
 * Relative phenological threshold
 *
 * threshold =
 *
 *     minimum NDVI +
 *     annual amplitude * phenology_fraction
 *
 * Example:
 *
 *     minimum NDVI = 0.20
 *     maximum NDVI = 0.80
 *     amplitude    = 0.60
 *
 *     fraction     = 0.25
 *
 *     threshold    = 0.20 + (0.60 * 0.25)
 *                  = 0.35
 */
var phenology_threshold = min_ndvi
    .add(
        ndvi_amplitude.multiply(
            phenology_fraction
        )
    )
    .rename('phenology_threshold');


/**
 * Active vegetation months
 *
 * Each image contains the month number only where
 * monthly NDVI is above the phenological threshold.
 *
 * All month images are explicitly cast to Int16
 * to create a homogeneous ImageCollection.
 */
var active_month_images = months.map(function (month) {
    var band_name =
        'ndvi_' +
        String(month).padStart(2, '0');

    var monthly_ndvi = ndvi_stack
        .select(band_name);

    var active = monthly_ndvi
        .gte(phenology_threshold);

    var month_image = ee.Image
        .constant(month)
        .rename('month')
        .toInt16();

    return month_image
        .updateMask(active)
        .updateMask(seasonal_mask)
        .clip(aoi);
});


/**
 * Active vegetation month collection
 */
var active_month_collection = ee.ImageCollection
    .fromImages(active_month_images)
    .map(function (image) {
        return image
            .select('month')
            .toInt16();
    });

print(
    'rojos → duración corta',
    'amarillos → duración intermedia',
    'verdes → duración larga'
);


/**
 * Start of active vegetation period
 *
 * First month where NDVI exceeds the relative threshold.
 */
var start_month = active_month_collection
    .min()
    .rename('start_month')
    .toByte();


/**
 * End of active vegetation period
 *
 * Last month where NDVI exceeds the relative threshold.
 */
var end_month = active_month_collection
    .max()
    .rename('end_month')
    .toByte();


/**
 * Approximate duration of active vegetation period
 *
 * Example:
 *
 *     Start = March
 *     End   = August
 *
 *     Duration = 8 - 3 + 1
 *              = 6 months
 */
var duration_months = end_month
    .toInt16()
    .subtract(
        start_month.toInt16()
    )
    .add(1)
    .rename('duration_months')
    .toByte();


/**
 * Annual true-color composite
 */
var annual_composite = get_clear_collection(
    year_start,
    year_end
)
    .median()
    .clip(aoi);


/**
 * Visualization
 */
var vis_rgb = {
    bands: [
        'B4',
        'B3',
        'B2',
    ],
    min: 0,
    max: 3000,
};

var vis_ndvi = {
    min: 0,
    max: 0.9,
    palette: [
        '#d73027',
        '#fee08b',
        '#d9ef8b',
        '#66bd63',
        '#1a9850',
        '#006837',
    ],
};

var vis_month = {
    min: 1,
    max: 12,
    palette: [
        '#313695',
        '#4575b4',
        '#74add1',
        '#abd9e9',
        '#e0f3f8',
        '#ffffbf',
        '#fee090',
        '#fdae61',
        '#f46d43',
        '#d73027',
        '#a50026',
        '#542788',
    ],
};

var vis_amplitude = {
    min: 0,
    max: 0.8,
    palette: [
        '#ffffcc',
        '#c2e699',
        '#78c679',
        '#31a354',
        '#006837',
    ],
};

var vis_duration = {
    min: 1,
    max: 12,
    palette: [
        '#a50026',
        '#d73027',
        '#f46d43',
        '#fdae61',
        '#fee08b',
        '#d9ef8b',
        '#a6d96a',
        '#66bd63',
        '#1a9850',
        '#063a7c',
    ],
};

var vis_valid_months = {
    min: 0,
    max: 12,
    palette: [
        '#d73027',
        '#fdae61',
        '#ffffbf',
        '#a6d96a',
        '#1a9850',
    ],
};


/**
 * Map
 */
Map.setOptions('SATELLITE');

Map.addLayer(
    annual_composite,
    vis_rgb,
    'Sentinel-2 annual composite',
    false
);

months.forEach(function (month, index) {
    var band_name =
        'ndvi_' +
        String(month).padStart(2, '0');

    Map.addLayer(
        ndvi_stack.select(band_name),
        vis_ndvi,
        'NDVI ' + month_names[index],
        false
    );
});

Map.addLayer(
    max_ndvi.updateMask(seasonal_mask),
    vis_ndvi,
    'Maximum annual NDVI',
    false
);

Map.addLayer(
    ndvi_amplitude.updateMask(seasonal_mask),
    vis_amplitude,
    'Annual NDVI amplitude',
    false
);

Map.addLayer(
    start_month,
    vis_month,
    'Start of vegetation activity',
    false
);

Map.addLayer(
    peak_month,
    vis_month,
    'Month of maximum NDVI',
    true
);

Map.addLayer(
    end_month,
    vis_month,
    'End of vegetation activity',
    false
);

Map.addLayer(
    duration_months,
    vis_duration,
    'Vegetation activity duration',
    false
);

Map.addLayer(
    valid_month_count,
    vis_valid_months,
    'Valid months',
    false
);

En la consola aparecen la ventana temporal y la leyenda de la capa de duración. En el mapa se cargan el compuesto anual, los doce NDVI mensuales y las capas de máximo, amplitud, inicio, mes de máximo, fin, duración y meses válidos; solo el mes de máximo queda visible al abrir. El script también arma monthly_summary, con el número de imágenes despejadas de cada mes: añade print(monthly_summary) si quieres revisarlo.

La lógica es idéntica. En Colab no hay herramienta de dibujo, así que el polígono se escribe como una lista de coordenadas [lon, lat] en aoi_coordinates. Los operadores lógicos van con mayúscula (.And()) porque and es una palabra reservada de Python.

ndvi_fenologia.py
"""
Google Earth Engine - Array Image example

Annual crop phenology using monthly Sentinel-2 NDVI.

Each pixel is represented as a temporal vector:

    [NDVI_Jan, NDVI_Feb, ..., NDVI_Dec]

The Array Image is used to derive the minimum and maximum annual NDVI,
the month of maximum NDVI and the annual NDVI amplitude. A relative
phenological threshold is then used to estimate the start, end and
approximate duration of the active season.

Important:
    Start and end months represent spectral phenology. They should not be
    interpreted directly as exact sowing or harvest dates.
"""

import ee
import geemap

ee.Authenticate()
ee.Initialize(project='tu-proyecto-gcp')


# Configuration
year = 2025
clear_threshold = 0.60
phenology_fraction = 0.25
min_amplitude = 0.15
min_valid_months = 8


# Area of interest: same polygon as "geometry" in the Code Editor,
# as [lon, lat] pairs.
aoi_coordinates = [
    # [lon, lat], [lon, lat], ...
]
aoi = ee.Geometry.Polygon([aoi_coordinates])


# Collections
s2 = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
cs_plus = ee.ImageCollection('GOOGLE/CLOUD_SCORE_PLUS/V1/S2_HARMONIZED')


# Months
months = list(range(1, 13))
month_names = [
    'January', 'February', 'March', 'April', 'May', 'June',
    'July', 'August', 'September', 'October', 'November', 'December',
]


# Annual temporal window
year_start = ee.Date.fromYMD(year, 1, 1)
year_end = year_start.advance(1, 'year')

print(
    'Temporal window:',
    year_start.format('YYYY-MM-dd').getInfo(),
    year_end.advance(-1, 'day').format('YYYY-MM-dd').getInfo(),
)


def get_clear_collection(start_date, end_date):
    """Returns a cloud-filtered Sentinel-2 collection."""

    def mask_clouds(image):
        clear_mask = image.select('cs_cdf').gte(clear_threshold)
        return image.updateMask(clear_mask)

    return (
        s2.filterBounds(aoi)
        .filterDate(start_date, end_date)
        .linkCollection(cs_plus, ['cs_cdf'])
        .map(mask_clouds)
    )


def get_monthly_ndvi(month):
    """Creates a monthly NDVI composite for one calendar month."""
    start_date = ee.Date.fromYMD(year, month, 1)
    end_date = start_date.advance(1, 'month')
    collection = get_clear_collection(start_date, end_date)

    empty_image = (
        ee.Image.constant([0, 0])
        .rename(['B8', 'B4'])
        .updateMask(ee.Image.constant(0))
    )

    composite = ee.Image(
        ee.Algorithms.If(
            collection.size().gt(0),
            collection.select(['B8', 'B4']).median(),
            empty_image,
        )
    )

    return (
        composite.normalizedDifference(['B8', 'B4'])
        .rename(f'ndvi_{month:02d}')
        .clip(aoi)
    )


# Monthly NDVI images
monthly_images = [get_monthly_ndvi(month) for month in months]


# Monthly image availability
def month_summary(month):
    start_date = ee.Date.fromYMD(year, month, 1)
    end_date = start_date.advance(1, 'month')
    image_count = get_clear_collection(start_date, end_date).size()

    return ee.Feature(None, {
        'month': month,
        'start_date': start_date.format('YYYY-MM-dd'),
        'end_date': end_date.advance(-1, 'day').format('YYYY-MM-dd'),
        'image_count': image_count,
    })


monthly_summary = ee.FeatureCollection([month_summary(m) for m in months])


# Multiband NDVI image: ndvi_01, ndvi_02, ..., ndvi_12
ndvi_stack = ee.Image.cat(monthly_images)


# Number of valid months per pixel
valid_month_count = (
    ndvi_stack.mask()
    .reduce(ee.Reducer.sum())
    .rename('valid_month_count')
    .clip(aoi)
)


# Array Image: each pixel becomes [NDVI_01, NDVI_02, ..., NDVI_12]
ndvi_array = ndvi_stack.toArray()


# Arrays used for robust temporal calculations. Missing values are
# replaced with -2 for maximum and 2 for minimum calculations, so they
# cannot become the annual maximum or minimum.
ndvi_array_max = ndvi_stack.unmask(-2).toArray()
ndvi_array_min = ndvi_stack.unmask(2).toArray()


# Maximum and minimum annual NDVI
max_ndvi = (
    ndvi_array_max.arrayReduce(ee.Reducer.max(), [0])
    .arrayGet([0])
    .rename('max_ndvi')
)

min_ndvi = (
    ndvi_array_min.arrayReduce(ee.Reducer.min(), [0])
    .arrayGet([0])
    .rename('min_ndvi')
)


# Annual NDVI amplitude: large values indicate stronger seasonal change
ndvi_amplitude = max_ndvi.subtract(min_ndvi).rename('ndvi_amplitude')


# Seasonal validity mask: enough valid months and enough annual variation
seasonal_mask = (
    valid_month_count.gte(min_valid_months)
    .And(ndvi_amplitude.gte(min_amplitude))
)


# Position of maximum annual NDVI (0 = January, ..., 11 = December)
max_position = (
    ndvi_array_max.arrayArgmax()
    .arrayGet([0])
    .rename('max_position')
)


# Month of maximum NDVI
peak_month = (
    max_position.add(1)
    .rename('peak_month')
    .toByte()
    .updateMask(seasonal_mask)
)


# Relative phenological threshold:
# minimum NDVI + annual amplitude * phenology_fraction
phenology_threshold = (
    min_ndvi.add(ndvi_amplitude.multiply(phenology_fraction))
    .rename('phenology_threshold')
)


# Active vegetation months: the month number only where monthly NDVI is
# above the threshold. Cast to Int16 to build a homogeneous collection.
def active_month(month):
    monthly_ndvi = ndvi_stack.select(f'ndvi_{month:02d}')
    active = monthly_ndvi.gte(phenology_threshold)
    month_image = ee.Image.constant(month).rename('month').toInt16()

    return (
        month_image.updateMask(active)
        .updateMask(seasonal_mask)
        .clip(aoi)
    )


active_month_collection = ee.ImageCollection.fromImages(
    [active_month(month) for month in months]
).map(lambda image: image.select('month').toInt16())

print(
    'rojos → duración corta',
    'amarillos → duración intermedia',
    'verdes → duración larga',
    sep='\n',
)


# Start and end of the active vegetation period: first and last month
# where NDVI exceeds the relative threshold
start_month = active_month_collection.min().rename('start_month').toByte()
end_month = active_month_collection.max().rename('end_month').toByte()


# Approximate duration of the active vegetation period
duration_months = (
    end_month.toInt16()
    .subtract(start_month.toInt16())
    .add(1)
    .rename('duration_months')
    .toByte()
)


# Annual true-color composite
annual_composite = get_clear_collection(year_start, year_end).median().clip(aoi)


# Visualization
vis_rgb = {'bands': ['B4', 'B3', 'B2'], 'min': 0, 'max': 3000}

vis_ndvi = {
    'min': 0,
    'max': 0.9,
    'palette': ['#d73027', '#fee08b', '#d9ef8b', '#66bd63', '#1a9850', '#006837'],
}

vis_month = {
    'min': 1,
    'max': 12,
    'palette': [
        '#313695', '#4575b4', '#74add1', '#abd9e9', '#e0f3f8', '#ffffbf',
        '#fee090', '#fdae61', '#f46d43', '#d73027', '#a50026', '#542788',
    ],
}

vis_amplitude = {
    'min': 0,
    'max': 0.8,
    'palette': ['#ffffcc', '#c2e699', '#78c679', '#31a354', '#006837'],
}

vis_duration = {
    'min': 1,
    'max': 12,
    'palette': [
        '#a50026', '#d73027', '#f46d43', '#fdae61', '#fee08b',
        '#d9ef8b', '#a6d96a', '#66bd63', '#1a9850', '#063a7c',
    ],
}

vis_valid_months = {
    'min': 0,
    'max': 12,
    'palette': ['#d73027', '#fdae61', '#ffffbf', '#a6d96a', '#1a9850'],
}


# Map
mapa = geemap.Map()
mapa.add_basemap('SATELLITE')
mapa.centerObject(aoi, 12)

mapa.addLayer(annual_composite, vis_rgb, 'Sentinel-2 annual composite', False)

for index, month in enumerate(months):
    mapa.addLayer(
        ndvi_stack.select(f'ndvi_{month:02d}'),
        vis_ndvi,
        f'NDVI {month_names[index]}',
        False,
    )

mapa.addLayer(max_ndvi.updateMask(seasonal_mask), vis_ndvi, 'Maximum annual NDVI', False)
mapa.addLayer(ndvi_amplitude.updateMask(seasonal_mask), vis_amplitude, 'Annual NDVI amplitude', False)
mapa.addLayer(start_month, vis_month, 'Start of vegetation activity', False)
mapa.addLayer(peak_month, vis_month, 'Month of maximum NDVI', True)
mapa.addLayer(end_month, vis_month, 'End of vegetation activity', False)
mapa.addLayer(duration_months, vis_duration, 'Vegetation activity duration', False)
mapa.addLayer(valid_month_count, vis_valid_months, 'Valid months', False)
mapa

Ciclos de cultivo al sur de Neiva en 2025

La capa de duración para 2025 separa con claridad las parcelas agrícolas del resto del paisaje en la zona agrícola al sur de Neiva, Huila, a orillas del río Magdalena. Las parcelas con ciclos cortos aparecen en rojo y naranja, las de ciclo intermedio en amarillo y las de ciclo largo en verde. En azul oscuro queda la vegetación que se mantiene activa casi todo el año, y los píxeles sin variación suficiente, como zonas urbanas, quedan sin color.

Mapa de duración de la temporada de vegetación: parcelas agrícolas en rojo, amarillo y verde sobre un fondo azul oscuro de vegetación activa todo el año, atravesado por un río.
Duración de la vegetación activa en 2025, en meses. Rojo: ciclo corto; amarillo: intermedio; verde: largo; azul oscuro: activa casi todo el año.
Periodo1 ene – 31 dic 2025
UbicaciónCultivos al sur de Neiva, Huila
Meses válidos mínimos8 de 12
Amplitud mínima0.15
Umbral fenológicomínimo + 25 % de la amplitud

Lo que la IA suele hacer mal

Si le pides este script a un asistente de IA, lo más probable es que funcione a medias. Estos son los fallos que más aparecen en el código generado, y cómo los evita el script de arriba:

  1. Convierte a array sin tratar los meses vacíos. toArray() enmascara el píxel entero si falta un mes, así que en zonas nubladas el mapa queda casi vacío. El script rellena los huecos antes de convertir.
  2. Rellena los meses vacíos con 0. Un 0 es un NDVI posible: el mes nublado pasa a ser el mínimo del año, la amplitud se infla y el inicio de la temporada se adelanta. El script usa −2 y 2, valores que el NDVI nunca alcanza.
  3. Olvida arrayGet([0]). Después de arrayReduce el resultado sigue siendo un array de un elemento, y restarlo o visualizarlo da errores de tipo.
  4. Mezcla tipos en la colección de meses. Si las imágenes de cada mes no tienen el mismo tipo de dato, min() y max() fallan. El script convierte todas a Int16.
  5. Usa COPERNICUS/S2_SR sin _HARMONIZED. Desde enero de 2022 Sentinel-2 cambió su procesamiento; sin la colección armonizada, una serie que cruce esa fecha muestra un salto de NDVI que no existe.
  6. Presenta el inicio y el fin como siembra y cosecha. Son fechas de cambio espectral. Para siembra y cosecha hace falta información de campo o calendarios de cultivo.

Cómo leer las capas

CapaQué muestra
Maximum annual NDVIEl verdor máximo del año; en cultivos, el momento de mayor desarrollo
Annual NDVI amplitudeCuánto cambia el píxel en el año; alto en cultivos, bajo en bosques y zonas urbanas
Month of maximum NDVIEl mes del pico de verdor, de 1 (enero) a 12 (diciembre)
Start / End of vegetation activityPrimer y último mes por encima del umbral fenológico
Vegetation activity durationMeses entre inicio y fin; rojo corto, amarillo intermedio, verde largo
Valid monthsMeses con dato despejado; si es bajo, el resultado es menos fiable

Como referencia para las capas de NDVI mensual y máximo (rangos orientativos; cómo interpretarlos sin caer en umbrales engañosos está en qué es el NDVI):

NDVILo que suele indicar
Menor que 0Agua, nubes residuales o nieve
0 a 0.2Suelo desnudo, zonas urbanas o parcelas recién cosechadas
0.2 a 0.5Pastos, cultivos jóvenes o vegetación dispersa
Mayor que 0.5Vegetación densa y activa, como bosques o cultivos en pleno desarrollo

Meses nublados, temporadas que cruzan el año y dobles cosechas

Muchos píxeles quedan sin resultado

Revisa la capa Valid months: en zonas muy nubladas es habitual que pocos píxeles lleguen a 8 meses válidos. Baja min_valid_months a 6 o el umbral de Cloud Score+ a 0.5, sabiendo que el resultado será menos robusto.

Temporadas que cruzan el cambio de año

El script busca el inicio y el fin dentro del año calendario. Un cultivo sembrado en noviembre y cosechado en febrero aparece con inicio en enero, fin en diciembre y 12 meses de duración. Si la temporada de tu zona cruza diciembre, desplaza la ventana de 12 meses para que empiece en un mes de baja actividad, por ejemplo de julio a junio.

Dos cosechas en el año

La duración va del primer al último mes activo, así que dos ciclos seguidos se suman en una sola temporada larga. Para contar solo los meses activos, usa active_month_collection.count() en lugar de fin − inicio + 1.

User memory limit exceeded o Computation timed out

Doce compuestos mensuales y las operaciones de array pesan más que un NDVI suelto. En áreas grandes, exporta las capas en lugar de visualizarlas a zoom bajo, o divide la zona en partes más pequeñas.

Preguntas frecuentes

¿Qué diferencia hay entre un NDVI puntual y una serie temporal de NDVI?

Un NDVI puntual muestra cuánto verdor hay en una fecha. Una serie temporal muestra cómo cambia a lo largo del año, y de esa curva salen la amplitud, el mes de máximo verdor y la duración de la temporada, que separan cultivos, pastos y bosques aunque tengan el mismo NDVI un día concreto.

¿Cada cuánto conviene calcular el NDVI para ver la fenología?

Un compuesto por mes es un buen punto de partida, porque reúne varias pasadas de Sentinel-2 y deja menos huecos por nubes. En cultivos de ciclo corto conviene bajar a quincenas, siempre que la capa Valid months muestre suficientes datos despejados.

¿Cómo se hace una serie temporal de NDVI en Google Earth Engine?

Se arma un compuesto sin nubes por periodo (aquí, uno por mes), se calcula el NDVI de cada uno y se apilan en una imagen de varias bandas. Con toArray() cada píxel pasa a tener su serie completa, lista para analizarla con operaciones de array.

¿Qué es una array image en Earth Engine?

Es una imagen en la que cada píxel guarda un vector o una matriz en lugar de un solo valor. Permite calcular sobre la serie completa de cada píxel, por ejemplo su máximo, su mínimo o la posición del máximo, con arrayReduce y arrayArgmax.

Fuentes

  1. Harmonized Sentinel-2 MSI: MultiSpectral Instrument, Level-2A (SR)Google Earth Engine, catálogo de datosBandas, longitudes de onda y procesamiento armonizado de Sentinel-2.
  2. Cloud Score+ S2_HARMONIZED V1Google Earth Engine, catálogo de datosLa puntuación de nubes de 0 a 1 que usa el script para enmascarar.
  3. Arrays and Array ImagesGoogle Earth Engine, documentaciónCómo funcionan toArray(), arrayReduce y arrayArgmax.
  4. Measuring Vegetation (NDVI & EVI)NASA Earth ObservatoryQué mide el NDVI y por qué.

Sigue aprendiendo

Mira esta técnica aplicada en el mapa de cobertura y uso del suelo de Venezuela. ¿Lo necesitas hecho y verificado para tu zona? Mira el servicio Mapas de incendios, agua y cobertura.

¿Tienes un script de Earth Engine que no termina o que necesitas llevar a Python?

Escribe a contacto@ventuari.dev con el enlace de tu script y recibirás una propuesta con alcance y precio.