¿Qué relación hay entre deforestación y quemas en Colombia?

/**
 * Fire after deforestation
 *
 * Confirmed GLAD-L alerts crossed with MapBiomas Fire Colombia.
 *
 * Study area:
 *     San Vicente del Caguán, Caquetá, Colombia.
 *
 * Products:
 *     - Confirmed GLAD alerts (January-April)
 *     - Alerts followed by fire
 *     - Monthly alert and burned area statistics
 *     - Fire occurrence before and after alerts
 *
 * Sources:
 *     GLAD-L
 *     MapBiomas Fuego Colombia Collection 1
 */


/**
 * Configuration
 */
var municipality = 'San Vicente Del Caguan';

var country = 'COL';

var alert_year = 2025;

var last_alert_month = 4;

var alerts_path = 'projects/glad/alert/UpdResult';

var fire_folder = 'projects/mapbiomas-public/assets/colombia/fire/collection1/';

var regions_path = 'FAO/GAUL/2025/level2';

var alert_id = String(alert_year).slice(-2);


/**
 * Study area
 */
var region = ee.FeatureCollection(regions_path)
    .filter(ee.Filter.eq('ISO3_CODE', country))
    .filter(ee.Filter.eq('GAUL2_NAME', municipality));

var geometry = region.geometry();


/**
 * Input datasets
 *
 * GLAD values:
 *     conf = 2 (probable)
 *     conf = 3 (confirmed)
 *
 * alertDate represents the day of the year.
 *
 * Images without the requested annual bands are excluded.
 */
var alerts = ee.ImageCollection(alerts_path)
    .filterBounds(geometry)
    .filter(
        ee.Filter.listContains(
            'system:band_names',
            'conf' + alert_id
        )
    );

var glad = alerts
    .select(
        ['conf' + alert_id, 'alertDate' + alert_id],
        ['conf', 'date']
    )
    .mosaic();

var fire = ee.Image(
    fire_folder + 'mapbiomas_colombia_fire_collection1_monthly_burned_v1'
)
    .select('burned_monthly_' + alert_year);


/**
 * Area calculation
 *
 * Uses the projection of a GLAD tile
 * and the actual area of each pixel.
 */
var projection = ee.Image(alerts.first())
    .select('conf' + alert_id)
    .projection();

var hectares = ee.Image.pixelArea().divide(10000);


/**
 * Calculate the sum of pixel values within the study area.
 *
 * @param {ee.Image} image - Image containing area values.
 * @returns {ee.Dictionary} Sum of values by band.
 */
function sum_hectares(image) {
    return image.reduceRegion({
        reducer: ee.Reducer.sum(),
        geometry: geometry,
        crs: projection,
        maxPixels: 1e13,
        tileScale: 4,
    });
}


/**
 * Alert month
 *
 * Converts the day of year into a month (1-12).
 */
var month_starts = [];

for (var month_number = 2; month_number <= 12; month_number++) {
    var start = new Date(
        Date.UTC(alert_year, month_number - 1, 1)
    );

    var day = (
        start - Date.UTC(alert_year, 0, 1)
    ) / 86400000 + 1;

    month_starts.push(day);
}

var alert_month = ee.Image(1);

month_starts.forEach(function (day) {
    alert_month = alert_month.add(
        glad.select('date').gte(day)
    );
});


/**
 * Fire month
 *
 * Values:
 *     1-12  = Landsat burn month
 *     13-24 = MODIS complement (month + 12)
 *
 * MODIS values are normalized to months 1-12.
 * Pixels without fire are assigned 0.
 */
var fire_month = fire
    .updateMask(fire.gte(1))
    .subtract(1)
    .mod(12)
    .add(1)
    .unmask(0);


/**
 * Alerts and fire
 *
 * confirmed:
 *     GLAD confirmed alerts from January to April.
 *
 * burned_after:
 *     Fire in the alert month or later.
 *
 * same_month:
 *     Fire in the same month as the alert.
 *
 * burned_before:
 *     Fire before the alert month.
 */
var confirmed = glad.select('conf')
    .eq(3)
    .and(alert_month.lte(last_alert_month))
    .unmask(0);

var burned_after = confirmed
    .and(fire_month.gte(alert_month));

var same_month = confirmed
    .and(fire_month.eq(alert_month));

var burned_before = confirmed
    .and(fire_month.gt(0))
    .and(fire_month.lt(alert_month));


/**
 * Area statistics
 */
var totals = sum_hectares(
    ee.Image.cat([
        confirmed.rename('alerts'),
        burned_after.rename('burned_after'),
        same_month.rename('same_month'),
        burned_before.rename('burned_before'),
    ])
        .multiply(hectares)
);

var total_alerts = ee.Number(totals.get('alerts'));

var total_burned_after = ee.Number(
    totals.get('burned_after')
);

var burned_percentage = ee.Algorithms.If(
    total_alerts.gt(0),
    total_burned_after
        .divide(total_alerts.max(1))
        .multiply(100),
    0
);


/**
 * Results
 */
print(
    'Alertas GLAD confirmadas, enero-abril ' + alert_year + ' (ha)',
    totals.get('alerts')
);

print(
    'Área quemada en el mes de alerta o posteriormente (ha)',
    totals.get('burned_after')
);

print(
    'Área quemada en el mismo mes de la alerta (ha)',
    totals.get('same_month')
);

print(
    'Área quemada antes del mes de la alerta (ha)',
    totals.get('burned_before')
);

print(
    'Porcentaje de alertas con quema posterior o simultánea (%)',
    burned_percentage
);


/**
 * Monthly statistics
 */
var month_labels = [
    'Ene', 'Feb', 'Mar', 'Abr',
    'May', 'Jun', 'Jul', 'Ago',
    'Sep', 'Oct', 'Nov', 'Dic',
];

var month_bands = [];

for (var month_index = 1; month_index <= 12; month_index++) {
    var month_id = (
        month_index < 10 ? '0' : ''
    ) + month_index;

    month_bands.push(
        confirmed
            .and(alert_month.eq(month_index))
            .rename('alert_' + month_id)
    );

    month_bands.push(
        burned_after
            .and(fire_month.eq(month_index))
            .rename('fire_' + month_id)
    );
}

var by_month = sum_hectares(
    ee.Image.cat(month_bands)
        .multiply(hectares)
);

var month_table = ee.FeatureCollection(
    month_labels.map(function (label, index) {
        var month_id = (
            index < 9 ? '0' : ''
        ) + (index + 1);

        return ee.Feature(null, {
            month: label,
            alerts: by_month.get('alert_' + month_id),
            fire: by_month.get('fire_' + month_id),
        });
    })
);


/**
 * Monthly chart
 */
var month_chart = ui.Chart.feature.byFeature(
    month_table,
    'month',
    ['alerts', 'fire']
)
    .setSeriesNames([
        'Alertas GLAD confirmadas',
        'Área quemada después de la alerta',
    ])
    .setChartType('ColumnChart')
    .setOptions({
        title: 'Alertas de deforestación y áreas quemadas | ' + alert_year,
        hAxis: {
            title: 'Mes',
        },
        vAxis: {
            title: 'Superficie (ha)',
        },
        legend: {
            position: 'bottom',
        },
        colors: [
            '#12202b',
            '#a1561b',
        ],
    });

print(month_chart);


/**
 * Visualization parameters
 */
var alert_vis = {
    palette: ['#449ae0'],
};

var burned_vis = {
    palette: ['#ff3b00'],
};

var boundary_vis = {
    color: '#ffffff',
    fillColor: '#00000000',
    width: 2,
};


/**
 * Map visualization
 */
Map.setOptions('HYBRID');

Map.centerObject(region, 9);

Map.addLayer(
    confirmed.selfMask().clip(geometry),
    alert_vis,
    'Alertas GLAD confirmadas | Ene-Abr ' + alert_year,
    true
);

Map.addLayer(
    burned_after.selfMask().clip(geometry),
    burned_vis,
    'Alertas con quema posterior o simultánea | ' + alert_year,
    true
);

Map.addLayer(
    region.style(boundary_vis),
    {},
    municipality,
    true
);


/**
 * Map legend
 */
var legend = ui.Panel({
    style: {
        position: 'bottom-left',
        padding: '10px 14px',
        backgroundColor: 'rgba(255,255,255,0.95)',
    },
});

var legend_title = ui.Label({
    value: 'Deforestación y fuego | ' + alert_year,
    style: {
        fontWeight: 'bold',
        fontSize: '13px',
        margin: '0 0 8px 0',
    },
});

var legend_items = [
    {
        label: 'Alertas GLAD confirmadas',
        color: '#449ae0',
    },
    {
        label: 'Alertas con quema posterior o simultánea',
        color: '#ff3b00',
    },
];


/**
 * Create a legend item.
 *
 * @param {string} label - Legend description.
 * @param {string} color - Category color.
 * @returns {ui.Panel} Legend row.
 */
function create_legend_item(label, color) {
    var color_box = ui.Label({
        style: {
            backgroundColor: color,
            padding: '8px',
            margin: '0 8px 4px 0',
        },
    });

    var description = ui.Label({
        value: label,
        style: {
            fontSize: '12px',
            margin: '0 0 4px 0',
        },
    });

    return ui.Panel({
        widgets: [
            color_box,
            description,
        ],
        layout: ui.Panel.Layout.flow('horizontal'),
    });
}


/**
 * Build legend
 */
var legend_widgets = [legend_title];

legend_items.forEach(function (item) {
    legend_widgets.push(
        create_legend_item(item.label, item.color)
    );
});

legend = ui.Panel({
    widgets: legend_widgets,
    style: {
        position: 'bottom-left',
        padding: '10px 14px',
        backgroundColor: 'rgba(255,255,255,0.95)',
    },
});

Map.add(legend);

"""
Fire after deforestation

Confirmed GLAD-L alerts crossed with MapBiomas Fire Colombia.

Study area:
    San Vicente del Caguán, Caquetá, Colombia.

Products:
    - Confirmed GLAD alerts (January-April)
    - Alerts followed by fire
    - Monthly alert and burned area statistics
    - Fire occurrence before and after alerts

Sources:
    GLAD-L
    MapBiomas Fuego Colombia Collection 1
"""
from datetime import date

import ee
import geemap
import matplotlib.pyplot as plt

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


# Configuration
municipality = 'San Vicente Del Caguan'

country = 'COL'

alert_year = 2025

last_alert_month = 4

alerts_path = 'projects/glad/alert/UpdResult'

fire_folder = 'projects/mapbiomas-public/assets/colombia/fire/collection1/'

regions_path = 'FAO/GAUL/2025/level2'

alert_id = str(alert_year)[-2:]


# Study area
region = (
    ee.FeatureCollection(regions_path)
    .filter(ee.Filter.eq('ISO3_CODE', country))
    .filter(ee.Filter.eq('GAUL2_NAME', municipality))
)

geometry = region.geometry()


# Input datasets
#
# GLAD values:
#     conf = 2 (probable)
#     conf = 3 (confirmed)
#
# alertDate represents the day of the year.
#
# Images without the requested annual bands are excluded.
alerts = (
    ee.ImageCollection(alerts_path)
    .filterBounds(geometry)
    .filter(
        ee.Filter.listContains(
            'system:band_names',
            f'conf{alert_id}'
        )
    )
)

glad = (
    alerts
    .select(
        [f'conf{alert_id}', f'alertDate{alert_id}'],
        ['conf', 'date']
    )
    .mosaic()
)

fire = (
    ee.Image(
        fire_folder + 'mapbiomas_colombia_fire_collection1_monthly_burned_v1'
    )
    .select(f'burned_monthly_{alert_year}')
)


# Area calculation
#
# Uses the projection of a GLAD tile
# and the actual area of each pixel.
projection = (
    ee.Image(alerts.first())
    .select(f'conf{alert_id}')
    .projection()
)

hectares = ee.Image.pixelArea().divide(10000)


def sum_hectares(image):
    """
    Calculate the sum of pixel values within the study area.

    Args:
        image (ee.Image): Image containing area values.

    Returns:
        ee.Dictionary: Sum of values by band.
    """
    return image.reduceRegion(
        reducer=ee.Reducer.sum(),
        geometry=geometry,
        crs=projection,
        maxPixels=1e13,
        tileScale=4,
    )


# Alert month
#
# Converts the day of year into a month (1-12).
month_starts = []

for month_number in range(2, 13):
    day = (
        date(alert_year, month_number, 1) - date(alert_year, 1, 1)
    ).days + 1

    month_starts.append(day)

alert_month = ee.Image(1)

for day in month_starts:
    alert_month = alert_month.add(
        glad.select('date').gte(day)
    )


# Fire month
#
# Values:
#     1-12  = Landsat burn month
#     13-24 = MODIS complement (month + 12)
#
# MODIS values are normalized to months 1-12.
# Pixels without fire are assigned 0.
fire_month = (
    fire
    .updateMask(fire.gte(1))
    .subtract(1)
    .mod(12)
    .add(1)
    .unmask(0)
)


# Alerts and fire
#
# confirmed:
#     GLAD confirmed alerts from January to April.
#
# burned_after:
#     Fire in the alert month or later.
#
# same_month:
#     Fire in the same month as the alert.
#
# burned_before:
#     Fire before the alert month.
confirmed = (
    glad.select('conf')
    .eq(3)
    .And(alert_month.lte(last_alert_month))
    .unmask(0)
)

burned_after = confirmed.And(fire_month.gte(alert_month))

same_month = confirmed.And(fire_month.eq(alert_month))

burned_before = (
    confirmed
    .And(fire_month.gt(0))
    .And(fire_month.lt(alert_month))
)


# Area statistics
totals = sum_hectares(
    ee.Image.cat([
        confirmed.rename('alerts'),
        burned_after.rename('burned_after'),
        same_month.rename('same_month'),
        burned_before.rename('burned_before'),
    ])
    .multiply(hectares)
).getInfo()

total_alerts = totals['alerts']

total_burned_after = totals['burned_after']

burned_percentage = (
    total_burned_after / total_alerts * 100
    if total_alerts > 0
    else 0
)


# Results
print(
    f'Alertas GLAD confirmadas, enero-abril {alert_year} (ha):',
    round(totals['alerts'], 1)
)

print(
    'Área quemada en el mes de alerta o posteriormente (ha):',
    round(totals['burned_after'], 1)
)

print(
    'Área quemada en el mismo mes de la alerta (ha):',
    round(totals['same_month'], 1)
)

print(
    'Área quemada antes del mes de la alerta (ha):',
    round(totals['burned_before'], 1)
)

print(
    'Porcentaje de alertas con quema posterior o simultánea (%):',
    round(burned_percentage, 1)
)


# Monthly statistics
month_labels = [
    'Ene', 'Feb', 'Mar', 'Abr',
    'May', 'Jun', 'Jul', 'Ago',
    'Sep', 'Oct', 'Nov', 'Dic',
]

month_bands = []

for month_index in range(1, 13):
    month_bands.append(
        confirmed
        .And(alert_month.eq(month_index))
        .rename(f'alert_{month_index:02d}')
    )

    month_bands.append(
        burned_after
        .And(fire_month.eq(month_index))
        .rename(f'fire_{month_index:02d}')
    )

by_month = sum_hectares(
    ee.Image.cat(month_bands)
    .multiply(hectares)
).getInfo()

month_alerts = [
    by_month[f'alert_{index:02d}'] for index in range(1, 13)
]

month_fire = [
    by_month[f'fire_{index:02d}'] for index in range(1, 13)
]


# Monthly chart
positions = range(len(month_labels))

bar_width = 0.4

figure, axis = plt.subplots(figsize=(10, 5))

axis.bar(
    [position - bar_width / 2 for position in positions],
    month_alerts,
    width=bar_width,
    color='#12202b',
    label='Alertas GLAD confirmadas',
)

axis.bar(
    [position + bar_width / 2 for position in positions],
    month_fire,
    width=bar_width,
    color='#a1561b',
    label='Área quemada después de la alerta',
)

axis.set_title(
    f'Alertas de deforestación y áreas quemadas | {alert_year}',
    loc='left',
    fontweight='bold',
)

axis.set_xticks(list(positions))
axis.set_xticklabels(month_labels)
axis.set_xlabel('Mes')
axis.set_ylabel('Superficie (ha)')
axis.grid(axis='y', color='#e3e7ea')
axis.set_axisbelow(True)
axis.spines[['top', 'right']].set_visible(False)

axis.legend(
    loc='upper center',
    bbox_to_anchor=(0.5, -0.15),
    ncol=2,
    frameon=False,
)

figure.tight_layout()

plt.show()


# Visualization parameters
alert_vis = {
    'palette': ['#449ae0'],
}

burned_vis = {
    'palette': ['#ff3b00'],
}

boundary_vis = {
    'color': '#ffffff',
    'fillColor': '#00000000',
    'width': 2,
}


# Map legend
legend_items = [
    {
        'label': 'Alertas GLAD confirmadas',
        'color': '#449ae0',
    },
    {
        'label': 'Alertas con quema posterior o simultánea',
        'color': '#ff3b00',
    },
]


# Map visualization
Map = geemap.Map()

Map.add_basemap('HYBRID')

Map.centerObject(region, 9)

Map.addLayer(
    confirmed.selfMask().clip(geometry),
    alert_vis,
    f'Alertas GLAD confirmadas | Ene-Abr {alert_year}',
    True
)

Map.addLayer(
    burned_after.selfMask().clip(geometry),
    burned_vis,
    f'Alertas con quema posterior o simultánea | {alert_year}',
    True
)

Map.addLayer(
    region.style(**boundary_vis),
    {},
    municipality,
    True
)

Map.add_legend(
    title=f'Deforestación y fuego | {alert_year}',
    legend_dict={item['label']: item['color'] for item in legend_items},
    position='bottomleft',
)

Map