Cómo hacer un mosaico Landsat sin nubes en Google Earth Engine

/**
 * Automatic monthly Landsat mosaic
 * Collection 2, Level 2
 *
 * WRS-2:
 *     Path: 2
 *     Row: 52
 *
 * The script automatically selects the clearest Landsat
 * scene available for each month and creates an annual
 * cloud-free median mosaic.
 *
 * Satellites:
 *     - Landsat 8
 *     - Landsat 9
 */


/**
 * Configuration
 */
var path = 2;

var row = 52;

var year = 2025;

// Initial scene-level cloud filter.
// A relatively permissive threshold is useful because the
// best scene will later be selected independently for each month.
var max_cloud_cover = 70;

var collections = [
    'LANDSAT/LC08/C02/T1_L2',
    'LANDSAT/LC09/C02/T1_L2',
];

var input_bands = [
    'SR_B2',
    'SR_B3',
    'SR_B4',
    'SR_B5',
    'SR_B6',
    'SR_B7',
];

var output_bands = [
    'blue',
    'green',
    'red',
    'nir',
    'swir1',
    'swir2',
];


/**
 * Dates
 */
var start_date = ee.Date.fromYMD(year, 1, 1);

var end_date = start_date.advance(1, 'year');   // filterDate excludes the end date


/**
 * Candidate Landsat images
 */
var raw_collection = get_collection(collections, start_date, end_date);

print(
    'Available image count',
    raw_collection.size()
);

print(
    'Available images',
    raw_collection.aggregate_array(
        'LANDSAT_PRODUCT_ID'
    )
);

print(
    'Available CLOUD_COVER',
    raw_collection.aggregate_array(
        'CLOUD_COVER'
    )
);


/**
 * Automatic monthly selection
 */
var selected_raw = select_best_monthly_images(
    raw_collection
);

print(
    'Automatically selected image count',
    selected_raw.size()
);

print(
    'Automatically selected images',
    selected_raw.aggregate_array(
        'LANDSAT_PRODUCT_ID'
    )
);

print(
    'Selected months',
    selected_raw.aggregate_array('selection_month')
);

print(
    'Selected CLOUD_COVER',
    selected_raw.aggregate_array('CLOUD_COVER')
);


/**
 * Apply pixel-level QA masking
 *
 * Only the automatically selected monthly scenes are processed.
 */
var selection = selected_raw
    .map(prepare_image)
    .sort('system:time_start');


/**
 * Scene geometry and projection
 *
 * A median composite has no native projection: Earth Engine
 * gives it WGS84 at 1 degree. The UTM projection of the scene
 * is kept so the export does not fall back to EPSG:4326.
 */
var first_image = ee.Image(raw_collection.first());

var zone = first_image.geometry();

var projection = first_image
    .select('SR_B4')
    .projection();


/**
 * Annual median mosaic
 */
var mosaic = selection
    .median()
    .setDefaultProjection(projection)
    .clip(zone);


/**
 * Cloud-free observation count
 *
 * Number of valid observations contributing to each
 * pixel of the annual mosaic (0 where none is left).
 */
var observations = selection
    .select('red')
    .count()
    .unmask(0)
    .rename('observations')
    .setDefaultProjection(projection)
    .clip(zone);


/**
 * Pixels without observations
 */
var no_observations = observations
    .eq(0)
    .selfMask()
    .rename('no_observations');


/**
 * Low-observation pixels
 *
 * Areas supported by only one valid observation
 * during the year.
 */
var low_observations = observations
    .eq(1)
    .selfMask()
    .rename('low_observations');


/**
 * Visualization
 */
var vis_true_color = {
    bands: ['red', 'green', 'blue'],
    min: 0,
    max: 0.25,
};

var vis_false_color = {
    bands: ['swir1', 'nir', 'red'],
    min: [0, 0, 0],
    max: [0.35, 0.45, 0.20],
};

var vis_observations = {
    min: 1,
    max: 12,
    palette: [
        'd73027',
        'fc8d59',
        'fee08b',
        'ffffbf',
        'd9ef8b',
        '91cf60',
        '1a9850',
        '006837',
    ],
};

var vis_no_observations = {
    palette: ['ff0000'],
};

var vis_low_observations = {
    palette: ['ff9900'],
};


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

Map.centerObject(zone, 9);

Map.addLayer(
    mosaic,
    vis_false_color,
    'False color (SWIR1, NIR, red)',
    true
);

Map.addLayer(
    mosaic,
    vis_true_color,
    'True color',
    false
);

Map.addLayer(
    observations,
    vis_observations,
    'Cloud-free observations',
    false
);

Map.addLayer(
    low_observations,
    vis_low_observations,
    'Only one valid observation',
    false
);

Map.addLayer(
    no_observations,
    vis_no_observations,
    'Pixels without observations',
    false
);


/**
 * Export
 *
 * Surface reflectance is multiplied by 10000 and
 * exported as signed 16-bit integers.
 *
 * No CRS is hardcoded: the export takes the CRS of the first
 * band, which is the scene's UTM projection set above.
 */
Export.image.toDrive({
    image: mosaic
        .multiply(10000)
        .round()
        .toInt16(),
    description: 'landsat_monthly_mosaic_' + path + '_' + row + '_' + year,
    folder: 'earthengine',
    region: zone,
    scale: 30,
    maxPixels: 1e10,
});


/**
 * Export observation count
 *
 * Optional quality layer showing how many cloud-free
 * observations support each pixel.
 */
Export.image.toDrive({
    image: observations.toUint8(),
    description: 'landsat_observations_' + path + '_' + row + '_' + year,
    folder: 'earthengine',
    region: zone,
    scale: 30,
    maxPixels: 1e10,
});


/**
 * Prepare Landsat image
 *
 * Applies Collection 2 Level 2 surface reflectance scaling
 * and masks invalid pixels using QA_PIXEL and QA_RADSAT.
 *
 * Masked conditions:
 *     - fill
 *     - dilated cloud
 *     - cirrus
 *     - cloud
 *     - cloud shadow
 *     - snow
 *     - radiometric saturation or terrain occlusion
 *
 * @param {ee.Image} image Landsat Collection 2 Level 2 image.
 * @returns {ee.Image} Cloud-free surface reflectance image.
 */
function prepare_image(image) {

    var qa = image.select('QA_PIXEL');

    var clear_mask = qa
        .bitwiseAnd(1 << 0)
        .eq(0)
        .and(qa.bitwiseAnd(1 << 1).eq(0))
        .and(qa.bitwiseAnd(1 << 2).eq(0))
        .and(qa.bitwiseAnd(1 << 3).eq(0))
        .and(qa.bitwiseAnd(1 << 4).eq(0))
        .and(qa.bitwiseAnd(1 << 5).eq(0));

    var saturation_mask = image
        .select('QA_RADSAT')
        .eq(0);

    var reflectance = image
        .select(input_bands, output_bands)
        .multiply(0.0000275)
        .add(-0.2)
        .updateMask(clear_mask)
        .updateMask(saturation_mask);

    // copyProperties returns an ee.Element: cast it back to ee.Image
    return ee.Image(
        reflectance.copyProperties(
            image,
            [
                'system:time_start',
                'CLOUD_COVER',
                'CLOUD_COVER_LAND',
                'WRS_PATH',
                'WRS_ROW',
                'SPACECRAFT_ID',
                'LANDSAT_PRODUCT_ID',
            ]
        )
    ).set('image_id', image.id());
}


/**
 * Load Landsat collections
 *
 * Loads all candidate Landsat scenes for the selected
 * WRS-2 path/row and year.
 *
 * @param {Array<string>} collection_ids Landsat collection IDs.
 * @param {ee.Date} start_date Start date.
 * @param {ee.Date} end_date End date.
 * @returns {ee.ImageCollection} Combined raw Landsat collection.
 */
function get_collection(collection_ids, start_date, end_date) {

    var collection = ee.ImageCollection([]);

    collection_ids.forEach(function (collection_id) {

        var landsat = ee.ImageCollection(collection_id)
            .filter(ee.Filter.eq('WRS_PATH', path))
            .filter(ee.Filter.eq('WRS_ROW', row))
            .filterDate(start_date, end_date)
            .filter(
                ee.Filter.lte('CLOUD_COVER', max_cloud_cover)
            );

        collection = collection.merge(landsat);
    });

    return collection
        .sort('system:time_start');
}


/**
 * Select best image of each month
 *
 * For each month, selects the candidate scene with the
 * lowest CLOUD_COVER metadata value.
 *
 * Only months with images are visited, so months without
 * candidates are skipped without conditionals or nulls.
 *
 * Selection happens before pixel-level QA masking to keep
 * computation lightweight.
 *
 * @param {ee.ImageCollection} collection Raw Landsat collection.
 * @returns {ee.ImageCollection} Best scene from each available month.
 */
function select_best_monthly_images(collection) {

    var with_month = collection.map(function (image) {
        return image.set(
            'selection_month',
            ee.Date(image.get('system:time_start')).get('month')
        );
    });

    var months = with_month
        .aggregate_array('selection_month')
        .distinct();

    var images = months.map(function (month) {
        return with_month
            .filter(ee.Filter.eq('selection_month', month))
            .sort('CLOUD_COVER')
            .first();
    });

    return ee.ImageCollection
        .fromImages(images)
        .sort('system:time_start');
}

"""
Automatic monthly Landsat mosaic
Collection 2, Level 2

WRS-2:
    Path: 2
    Row: 52

The script automatically selects the clearest Landsat
scene available for each month and creates an annual
cloud-free median mosaic.

Satellites:
    - Landsat 8
    - Landsat 9
"""

import ee
import geemap

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


# Configuration
path = 2
row = 52
year = 2025

# Initial scene-level cloud filter.
# A relatively permissive threshold is useful because the
# best scene will later be selected independently for each month.
max_cloud_cover = 70

collections = [
    'LANDSAT/LC08/C02/T1_L2',
    'LANDSAT/LC09/C02/T1_L2',
]

input_bands = ['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7']
output_bands = ['blue', 'green', 'red', 'nir', 'swir1', 'swir2']


def prepare_image(image):
    """
    Applies Collection 2 Level 2 surface reflectance scaling
    and masks invalid pixels using QA_PIXEL and QA_RADSAT.

    Masked conditions: fill, dilated cloud, cirrus, cloud,
    cloud shadow, snow, radiometric saturation or terrain occlusion.
    """
    qa = image.select('QA_PIXEL')

    clear_mask = (
        qa.bitwiseAnd(1 << 0).eq(0)
        .And(qa.bitwiseAnd(1 << 1).eq(0))
        .And(qa.bitwiseAnd(1 << 2).eq(0))
        .And(qa.bitwiseAnd(1 << 3).eq(0))
        .And(qa.bitwiseAnd(1 << 4).eq(0))
        .And(qa.bitwiseAnd(1 << 5).eq(0))
    )

    saturation_mask = image.select('QA_RADSAT').eq(0)

    reflectance = (
        image.select(input_bands, output_bands)
        .multiply(0.0000275)
        .add(-0.2)
        .updateMask(clear_mask)
        .updateMask(saturation_mask)
    )

    # copyProperties returns an ee.Element: cast it back to ee.Image
    return ee.Image(
        reflectance.copyProperties(
            image,
            [
                'system:time_start',
                'CLOUD_COVER',
                'CLOUD_COVER_LAND',
                'WRS_PATH',
                'WRS_ROW',
                'SPACECRAFT_ID',
                'LANDSAT_PRODUCT_ID',
            ],
        )
    ).set('image_id', image.id())


def get_collection(collection_ids, start_date, end_date):
    """Loads all candidate Landsat scenes for the WRS-2 path/row and year."""
    collection = ee.ImageCollection([])
    for collection_id in collection_ids:
        landsat = (
            ee.ImageCollection(collection_id)
            .filter(ee.Filter.eq('WRS_PATH', path))
            .filter(ee.Filter.eq('WRS_ROW', row))
            .filterDate(start_date, end_date)
            .filter(ee.Filter.lte('CLOUD_COVER', max_cloud_cover))
        )
        collection = collection.merge(landsat)
    return collection.sort('system:time_start')


def select_best_monthly_images(collection):
    """
    For each month, selects the candidate scene with the lowest
    CLOUD_COVER. Only months with images are visited.
    """
    with_month = collection.map(
        lambda image: image.set(
            'selection_month',
            ee.Date(image.get('system:time_start')).get('month'),
        )
    )
    months = with_month.aggregate_array('selection_month').distinct()
    images = months.map(
        lambda month: with_month
        .filter(ee.Filter.eq('selection_month', month))
        .sort('CLOUD_COVER')
        .first()
    )
    return ee.ImageCollection.fromImages(images).sort('system:time_start')


# Dates
start_date = ee.Date.fromYMD(year, 1, 1)
end_date = start_date.advance(1, 'year')   # filterDate excludes the end date


# Candidate Landsat images
raw_collection = get_collection(collections, start_date, end_date)

print('Available image count:', raw_collection.size().getInfo())
print('Available images:', raw_collection.aggregate_array('LANDSAT_PRODUCT_ID').getInfo())
print('Available CLOUD_COVER:', raw_collection.aggregate_array('CLOUD_COVER').getInfo())


# Automatic monthly selection
selected_raw = select_best_monthly_images(raw_collection)

print('Automatically selected image count:', selected_raw.size().getInfo())
print('Automatically selected images:', selected_raw.aggregate_array('LANDSAT_PRODUCT_ID').getInfo())
print('Selected months:', selected_raw.aggregate_array('selection_month').getInfo())
print('Selected CLOUD_COVER:', selected_raw.aggregate_array('CLOUD_COVER').getInfo())


# Pixel-level QA masking, only on the selected monthly scenes
selection = selected_raw.map(prepare_image).sort('system:time_start')


# Scene geometry and projection.
# A median composite has no native projection: Earth Engine gives it
# WGS84 at 1 degree. The UTM projection of the scene is kept so the
# export does not fall back to EPSG:4326.
first_image = ee.Image(raw_collection.first())
zone = first_image.geometry()
projection = first_image.select('SR_B4').projection()


# Annual median mosaic
mosaic = selection.median().setDefaultProjection(projection).clip(zone)


# Cloud-free observations per pixel (0 where none is left)
observations = (
    selection.select('red')
    .count()
    .unmask(0)
    .rename('observations')
    .setDefaultProjection(projection)
    .clip(zone)
)
no_observations = observations.eq(0).selfMask().rename('no_observations')
low_observations = observations.eq(1).selfMask().rename('low_observations')


# Export: reflectance x 10000 as signed 16-bit integers. No CRS is
# hardcoded: the export takes the CRS of the first band (the scene's UTM)
task_mosaic = ee.batch.Export.image.toDrive(
    image=mosaic.multiply(10000).round().toInt16(),
    description=f'landsat_monthly_mosaic_{path}_{row}_{year}',
    folder='earthengine',
    region=zone,
    scale=30,
    maxPixels=1e10,
)
task_mosaic.start()

# Optional quality layer: cloud-free observations per pixel
task_observations = ee.batch.Export.image.toDrive(
    image=observations.toUint8(),
    description=f'landsat_observations_{path}_{row}_{year}',
    folder='earthengine',
    region=zone,
    scale=30,
    maxPixels=1e10,
)
task_observations.start()


# Map
vis_true_color = {'bands': ['red', 'green', 'blue'], 'min': 0, 'max': 0.25}
vis_false_color = {'bands': ['swir1', 'nir', 'red'], 'min': [0, 0, 0], 'max': [0.35, 0.45, 0.20]}
vis_observations = {
    'min': 1,
    'max': 12,
    'palette': ['d73027', 'fc8d59', 'fee08b', 'ffffbf', 'd9ef8b', '91cf60', '1a9850', '006837'],
}

mapa = geemap.Map()
mapa.add_basemap('SATELLITE')
mapa.centerObject(zone, 9)
mapa.addLayer(mosaic, vis_false_color, 'False color (SWIR1, NIR, red)', True)
mapa.addLayer(mosaic, vis_true_color, 'True color', False)
mapa.addLayer(observations, vis_observations, 'Cloud-free observations', False)
mapa.addLayer(low_observations, {'palette': ['ff9900']}, 'Only one valid observation', False)
mapa.addLayer(no_observations, {'palette': ['ff0000']}, 'Pixels without observations', False)
mapa