Temperatura, precipitación y vegetación#
Este cuaderno muestra cómo representar datos ambientales en mapas mediante Google Earth Engine y crear gráficos de tendencias a lo largo del tiempo.
import folium
import ee
import geemap
import geopandas as gpd
ee.Authenticate()
# Escribe aquí el ID de tu proyecto, entre comillas
ee.Initialize(project="emerge-lessons")
Visualicemos los datos de Google Earth Engine#
Primero, elige una fecha inicial y una fecha final. Estas fechas se utilizarán para filtrar los datos. Después, crea un mapa vacío.
start_date = "2025-01-01"
end_date = "2025-01-31"
Podemos crear un mapa vacío centrado en Florida con una sola instrucción:
map = folium.Map(
location=[28.263363, -83.497652],
tiles="Cartodb dark_matter",
zoom_start=7
)
A continuación, definimos una función, adaptada de este tutorial, para agregar datos de Google Earth Engine a un mapa y mostrarlos de manera interactiva.
def add_ee_layer(self, ee_image_object, vis_params, name):
"""
Agrega a un mapa de Folium una capa de teselas
procedente de una imagen de Earth Engine.
"""
map_id_dict = ee.Image(ee_image_object).getMapId(vis_params)
folium.raster_layers.TileLayer(
tiles=map_id_dict["tile_fetcher"].url_format,
attr=(
'Map Data © '
'<a href="https://earthengine.google.com/">'
"Google Earth Engine</a>"
),
name=name,
overlay=True,
control=True
).add_to(self)
folium.Map.add_ee_layer = add_ee_layer
También cargaremos el límite de Florida. El archivo original de límites estatales se descargó de la Oficina del Censo de Estados Unidos y se filtró para conservar solamente Florida. Utilizaremos este límite para recortar los datos y concentrarnos en el estado.
fl = gpd.read_file(
"https://github.com/geo-di-lab/emerge-lessons/raw/refs/heads/main/docs/data/florida_boundary.geojson"
)[["geometry"]]
region = geemap.geopandas_to_ee(fl)
Temperatura de la superficie terrestre (LST)#
MOD11A1.061 Terra Land Surface Temperature and Emissivity Daily Global 1km
De este conjunto de datos utilizaremos LST_Day_1km, que representa la temperatura diurna de la superficie terrestre y se mide en kelvin.
lst = (
ee.ImageCollection("MODIS/061/MOD11A1")
.filterDate(start_date, end_date)
.select("LST_Day_1km")
.median()
.clip(region)
)
lst_vis = {
"min": 13000.0,
"max": 16500.0,
"palette": [
"040274", "040281", "0502a3", "0502b8", "0502ce", "0502e6",
"0602ff", "235cb1", "307ef3", "269db1", "30c8e2", "32d3ef",
"3be285", "3ff38f", "86e26f", "3ae237", "b5e22e", "d6e21f",
"fff705", "ffd611", "ffb613", "ff8b13", "ff6e08", "ff500d",
"ff0000", "de0101", "c21301", "a71001", "911003"
]
}
map.add_ee_layer(
lst,
lst_vis,
"Temperatura de la superficie terrestre"
)
display(map)
Imagen satelital#
Harmonized Sentinel-2 MSI: MultiSpectral Instrument, Level-2A (SR)
def mask_s2_clouds(image):
"""
Oculta las nubes de una imagen Sentinel-2 mediante la banda QA.
Parámetros:
image (ee.Image): imagen de Sentinel-2.
Devuelve:
ee.Image: imagen de Sentinel-2 con las nubes ocultas.
"""
qa = image.select("QA60")
# Los bits 10 y 11 representan nubes y cirros, respectivamente
cloud_bit_mask = 1 << 10
cirrus_bit_mask = 1 << 11
# Ambos indicadores deben ser cero para representar condiciones despejadas
mask = (
qa.bitwiseAnd(cloud_bit_mask)
.eq(0)
.And(
qa.bitwiseAnd(cirrus_bit_mask).eq(0)
)
)
return image.updateMask(mask).divide(10000)
rgb = (
ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
.filterDate(start_date, end_date)
.filter(
ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 20)
)
.map(mask_s2_clouds)
.median()
.clip(region)
)
rgb_vis = {
"min": 0.0,
"max": 0.3,
"bands": ["B4", "B3", "B2"]
}
map.add_ee_layer(rgb, rgb_vis, "Sentinel-2 RGB")
display(map)
La capa puede tardar un poco en cargarse. Es posible que observes espacios vacíos entre las imágenes satelitales. Esto ocurre porque la capa está formada por varias imágenes independientes que se unen para crear un mosaico. Las imágenes con más de un 20 % de nubosidad se excluyeron mediante el filtro.
Índice de Vegetación de Diferencia Normalizada (NDVI)#
Podemos utilizar la imagen de Sentinel-2 para calcular el Índice de Vegetación de Diferencia Normalizada (NDVI). Este valor, que varía entre -1 y 1, representa los niveles de vegetación. El NDVI se calcula a partir de las bandas de la imagen satelital:
\(NDVI = \frac{NIR - Red}{NIR + Red}\)
ndvi = rgb.normalizedDifference(
["B8", "B4"]
).rename("NDVI")
Sentinel-2 contiene 12 bandas. La banda 8 corresponde al infrarrojo cercano, o NIR, y la banda 4 corresponde al rojo. En la celda anterior, la función .normalizedDifference calcula:
\(NDVI = \frac{B8 - B4}{B8 + B4}\)
ndvi_vis = {
"min": -1,
"max": 1,
"palette": ["blue", "white", "green"]
}
map.add_ee_layer(ndvi, ndvi_vis, "NDVI")
display(map)
Índice de Agua de Diferencia Normalizada (NDWI)#
El Índice de Agua de Diferencia Normalizada (NDWI) mide la presencia de agua. Se calcula de la siguiente manera:
\(NDWI = \frac{B3 - B8}{B3 + B8}\)
ndwi = rgb.normalizedDifference(
["B3", "B8"]
).rename("NDWI")
ndwi_vis = {
"min": -1,
"max": 1,
"palette": ["white", "blue"]
}
map.add_ee_layer(ndwi, ndwi_vis, "NDWI")
display(map)
Precipitación#
CHIRPS Daily: Climate Hazards Center InfraRed Precipitation With Station Data (Version 2.0 Final)
rain = (
ee.ImageCollection("UCSB-CHG/CHIRPS/DAILY")
.filterDate(start_date, end_date)
.select("precipitation")
.sum()
.clip(region)
)
rain_vis = {
"min": 0,
"max": 100,
"palette": ["white", "blue"]
}
map.add_ee_layer(rain, rain_vis, "Precipitación")
display(map)
Mapa final#
Mediante folium.LayerControl, podemos agregar una opción para seleccionar los datos que queremos mostrar. En el mapa interactivo, haz clic en el menú de la esquina superior derecha para elegir las capas.
Ten en cuenta que la última capa agregada se muestra encima de las demás. Por esta razón, cuando todas las capas estén activadas, Precipitación aparecerá como la capa superior.
folium.LayerControl(
collapsed=False
).add_to(map)
display(map)
Datos a lo largo del tiempo#
import pandas as pd
import geopandas as gpd
import matplotlib
import matplotlib.pyplot as plt
from datetime import datetime
import numpy as np
from pylab import *
pd.set_option("display.max_columns", None)
mosquito = gpd.read_file(
"https://github.com/geo-di-lab/emerge-lessons/raw/refs/heads/main/docs/data/globe_mosquito.zip"
)
# Definir un conjunto de coordenadas en Florida
longitude = -80.57107
latitude = 25.48361
# Obtener imágenes satelitales históricas y recientes cerca del punto
coords = [longitude, latitude]
point = ee.Geometry.Point(coords)
sentinel2 = (
ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
.map(mask_s2_clouds)
.getRegion(point, 500)
.getInfo()
)
# Crear una tabla
sentinel2 = pd.DataFrame(
sentinel2[1:],
columns=sentinel2[0]
)
sentinel2.head()
# Calcular el NDVI y el NDWI
sentinel2["time"] = pd.to_datetime(
sentinel2["id"].str[0:15],
format="%Y%m%dT%H%M%S"
)
sentinel2["NDVI"] = (
(sentinel2["B8"] - sentinel2["B4"])
/ (sentinel2["B8"] + sentinel2["B4"])
)
sentinel2["NDWI"] = (
(sentinel2["B3"] - sentinel2["B8"])
/ (sentinel2["B3"] + sentinel2["B8"])
)
# Interpolar para sustituir los valores faltantes
sentinel2["NDVI"] = sentinel2["NDVI"].interpolate()
# Crear un gráfico del NDVI a lo largo del tiempo
plt.figure(figsize=(10, 6))
plt.plot(sentinel2["time"], sentinel2["NDVI"])
plt.title(
"Índice de Vegetación de Diferencia Normalizada (NDVI)"
)
plt.xlabel("Fecha")
plt.ylabel("NDVI")
plt.tight_layout()
plt.show()
Podemos observar que el NDVI cambia cada año con las estaciones, aumentando y disminuyendo a lo largo del tiempo. Por esta razón, puede ser útil representar los datos de otra manera y mostrar el NDVI de cada mes, en lugar de utilizar una sola línea que abarque varios años.
# Agregar columnas para el año y el mes
sentinel2["year"] = sentinel2["time"].dt.year
sentinel2["month"] = sentinel2["time"].dt.month
# Obtener los años únicos
years = sorted(sentinel2["year"].dropna().unique())
# Configurar el gráfico
fig, ax = plt.subplots(figsize=(10, 6))
# Dibujar una línea por año con colores diferentes
cmap = plt.get_cmap("YlGnBu")
colors = [
cmap(i / len(years))
for i in range(len(years))
]
for i, year in enumerate(years):
year_data = sentinel2[
sentinel2["year"] == year
]
monthly_avg = (
year_data.groupby("month")["NDVI"].mean()
)
plt.plot(
monthly_avg.index,
monthly_avg.values,
label=str(year),
color=colors[i]
)
# Agregar una barra de color
sm = matplotlib.cm.ScalarMappable(
cmap=cmap,
norm=matplotlib.colors.Normalize(
vmin=min(years),
vmax=max(years)
)
)
sm.set_array([])
cbar = plt.colorbar(
sm,
orientation="vertical",
pad=0.02,
ax=ax
)
cbar.set_label("Año")
cbar.set_ticks([min(years), max(years)])
cbar.set_ticklabels(
[str(min(years)), str(max(years))]
)
# Personalizar el gráfico
plt.title("NDVI mensual (2018–2025)")
plt.xlabel("Mes")
plt.ylabel(
"Índice de Vegetación de Diferencia Normalizada (NDVI)"
)
plt.xticks(range(1, 13))
plt.tight_layout()
plt.show()
# Crear una función para representar un conjunto de datos
# a lo largo del tiempo
def plot_over_time(longitude, latitude, name):
coords = [longitude, latitude]
point = ee.Geometry.Point(coords)
if name in ("NDVI", "NDWI"):
image = (
ee.ImageCollection(
"COPERNICUS/S2_SR_HARMONIZED"
)
.map(mask_s2_clouds)
.getRegion(point, scale=10)
.getInfo()
)
# Crear una tabla
image = pd.DataFrame(
image[1:],
columns=image[0]
)
image[["B8", "B4", "B3"]] = (
image[["B8", "B4", "B3"]]
.interpolate()
)
image["time"] = pd.to_datetime(
image["id"].str[0:15],
format="%Y%m%dT%H%M%S"
)
image["NDVI"] = (
(image["B8"] - image["B4"])
/ (image["B8"] + image["B4"])
)
image["NDWI"] = (
(image["B3"] - image["B8"])
/ (image["B3"] + image["B8"])
)
ylabel = name
elif name == "LST":
image = (
ee.ImageCollection("MODIS/061/MOD11A1")
.getRegion(point, scale=1000)
.getInfo()
)
# Crear una tabla
image = pd.DataFrame(
image[1:],
columns=image[0]
)
image["time"] = pd.to_datetime(
image["id"],
format="%Y_%m_%d"
)
# Convertir kelvin a grados Fahrenheit
# Los datos incluyen un factor de escala de 0.02
scale_value = 0.02
image["LST"] = (
(
image["LST_Day_1km"].interpolate()
* scale_value
- 273.15
)
* 1.8
+ 32
)
ylabel = (
"Temperatura de la superficie terrestre "
"(grados Fahrenheit)"
)
elif name == "Precipitation":
image = (
ee.ImageCollection(
"UCSB-CHG/CHIRPS/DAILY"
)
.select("precipitation")
.getRegion(point, scale=5566)
.getInfo()
)
# Crear una tabla
image = pd.DataFrame(
image[1:],
columns=image[0]
)
image["time"] = pd.to_datetime(
image["id"],
format="%Y%m%d"
)
image["Precipitation"] = (
image["precipitation"].interpolate()
)
ylabel = "Precipitación (mm/día)"
else:
raise ValueError(
"El nombre debe ser NDVI, NDWI, LST "
"o Precipitation."
)
# Agregar columnas para el año y el mes
image["year"] = image["time"].dt.year
image["month"] = image["time"].dt.month
# Obtener los años únicos
years = sorted(image["year"].dropna().unique())
if not years:
raise ValueError(
"No se encontraron datos con fechas válidas."
)
# Crear los gráficos
fig, (ax1, ax2) = plt.subplots(
2,
figsize=(10, 8)
)
ax1.plot(image["time"], image[name])
ax1.set_title(
f"{name} a lo largo del tiempo"
)
ax1.set_xlabel("Fecha")
ax1.set_ylabel(ylabel)
# Dibujar una línea por año con colores diferentes
cmap = plt.get_cmap("YlGnBu")
colors = [
cmap(i / len(years))
for i in range(len(years))
]
for i, year in enumerate(years):
year_data = image[
image["year"] == year
]
monthly_avg = (
year_data.groupby("month")[name].mean()
)
ax2.plot(
monthly_avg.index,
monthly_avg.values,
label=str(year),
color=colors[i]
)
# Agregar una barra de color
sm = matplotlib.cm.ScalarMappable(
cmap=cmap,
norm=matplotlib.colors.Normalize(
vmin=min(years),
vmax=max(years)
)
)
sm.set_array([])
cbar = fig.colorbar(
sm,
orientation="vertical",
pad=0.02,
ax=ax2
)
cbar.set_label("Año")
cbar.set_ticks([min(years), max(years)])
cbar.set_ticklabels(
[str(min(years)), str(max(years))]
)
# Personalizar el gráfico
ax2.set_title(
f"{name} mensual "
f"({min(years)}–{max(years)})"
)
ax2.set_xlabel("Mes")
ax2.set_ylabel(ylabel)
ax2.set_xticks(range(1, 13))
plt.tight_layout()
plt.show()
plot_over_time(
-80.469649,
26.030573,
"LST"
)
plot_over_time(
-80.469649,
26.030573,
"Precipitation"
)
plot_over_time(
-80.469649,
26.030573,
"NDVI"
)
plot_over_time(
-80.469649,
26.030573,
"NDWI"
)