Ambiente y hábitats#
¿Cuáles son las condiciones ambientales generales en los lugares donde se realizaron las observaciones de mosquitos?
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import seaborn as sns
import ee
import geemap
pd.set_option("display.max_columns", None)
ee.Authenticate()
# Escribe aquí el ID de tu proyecto, entre comillas
ee.Initialize(project="emerge-lessons")
mosquito = gpd.read_file(
"https://github.com/geo-di-lab/emerge-lessons/raw/"
"refs/heads/main/docs/data/globe_mosquito.zip"
)
Para cada uno de los conjuntos de datos de Google Earth Engine que utilizamos en la lección anterior, agregaremos los valores ambientales como columnas adicionales al conjunto de datos de mosquitos.
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)
def get_data_for_point(point):
"""
Agrega a una observación de mosquitos los valores ambientales
correspondientes a su ubicación y fecha.
"""
# Extraer la fecha de medición del punto
point_date = ee.Date(point.get("MeasuredDate"))
day_before = point_date.advance(-1, "day")
month_before = point_date.advance(-30, "day")
# Temperatura de la superficie terrestre
lst = (
ee.ImageCollection("MODIS/061/MOD11A1")
.filterDate(day_before, point_date)
.select(
[
"LST_Day_1km",
"QC_Day",
"LST_Night_1km",
"QC_Night",
"Clear_day_cov",
"Clear_night_cov"
]
)
.median()
# Obtener los valores de LST en la ubicación del punto
.reduceRegion(
reducer=ee.Reducer.median(),
geometry=point.geometry(),
scale=1000
)
)
point = point.set(lst)
# Sentinel-2: mediana del mes anterior
# Se utilizará para calcular el NDVI y el NDWI
sentinel2 = (
ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
.filterDate(month_before, point_date)
.map(mask_s2_clouds)
.median()
# Obtener los valores medianos de las bandas en el punto
.reduceRegion(
reducer=ee.Reducer.median(),
geometry=point.geometry(),
scale=10
)
)
point = point.set(sentinel2)
# Precipitación
rain = (
ee.ImageCollection("UCSB-CHG/CHIRPS/DAILY")
.filterDate(day_before, point_date)
.select("precipitation")
.sum()
# Obtener el valor en la ubicación del punto
.reduceRegion(
reducer=ee.Reducer.median(),
geometry=point.geometry(),
scale=5566
)
)
point = point.set(rain)
return point
Como el conjunto de datos completo requiere mucha memoria, utilizaremos solamente las observaciones de 2024. Este subconjunto contiene 2,039 observaciones, en comparación con las 43,009 observaciones totales recopiladas entre 2018 y 2024.
mosquito_2024 = mosquito[
mosquito["MeasuredDate"].str.startswith("2024")
]
# Aplicar la función a todos los puntos de 2024
points_fc = geemap.geopandas_to_ee(mosquito_2024)
result_fc = points_fc.map(get_data_for_point)
data = geemap.ee_to_gdf(result_fc)
# Calcular el NDVI
data["NDVI"] = (
(data["B8"] - data["B4"])
/ (data["B8"] + data["B4"])
)
# Calcular el NDWI
data["NDWI"] = (
(data["B3"] - data["B8"])
/ (data["B3"] + data["B8"])
)
data
Si el código anterior produce errores, algo que puede ocurrir debido a las limitaciones de memoria, ejecuta la siguiente celda para cargar directamente una copia de los datos procesados disponible en GitHub.
data = gpd.read_file(
"https://github.com/geo-di-lab/emerge-lessons/raw/"
"refs/heads/main/docs/data/globe_mosquito_env_2024.zip"
)
# Eliminar los valores atípicos de los conteos positivos
# mediante el método del rango intercuartílico (IQR)
positive_larvae = data[
data["LarvaeCountProcessed"] > 0
]
q1 = positive_larvae["LarvaeCountProcessed"].quantile(0.25)
q3 = positive_larvae["LarvaeCountProcessed"].quantile(0.75)
iqr = q3 - q1
threshold = 1.5
lower_bound = q1 - threshold * iqr
upper_bound = q3 + threshold * iqr
# Conservar las observaciones sin larvas para poder compararlas.
# Aplicar el filtro de valores atípicos solo a los conteos positivos.
data_cleaned = data[
(data["LarvaeCountProcessed"] == 0)
| (
(data["LarvaeCountProcessed"] >= lower_bound)
& (data["LarvaeCountProcessed"] <= upper_bound)
)
]
Gráficos#
fig, (ax1, ax2) = plt.subplots(2, figsize=(8, 7))
ax1.hist(
data_cleaned[
data_cleaned["LarvaeCountProcessed"] > 0
]["NDVI"],
bins=20
)
ax1.set_title("Con larvas")
ax1.set_xlabel("NDVI")
ax1.set_ylabel("Frecuencia")
ax2.hist(
data_cleaned[
data_cleaned["LarvaeCountProcessed"] == 0
]["NDVI"],
bins=20
)
ax2.set_title("Sin larvas")
ax2.set_xlabel("NDVI")
ax2.set_ylabel("Frecuencia")
plt.suptitle(
"Distribución de la vegetación (NDVI) "
"en las observaciones de mosquitos"
)
plt.tight_layout()
plt.show()
fig, (ax1, ax2) = plt.subplots(2, figsize=(8, 7))
ax1.hist(
data_cleaned[
data_cleaned["LarvaeCountProcessed"] > 0
]["precipitation"],
bins=200
)
ax1.set_title("Con larvas")
ax1.set_xlim([1, 20])
ax1.set_ylim([0, 50])
ax1.set_xlabel("Precipitación (mm/día)")
ax1.set_ylabel("Frecuencia")
ax2.hist(
data_cleaned[
data_cleaned["LarvaeCountProcessed"] == 0
]["precipitation"],
bins=200
)
ax2.set_title("Sin larvas")
ax2.set_xlim([1, 20])
ax2.set_ylim([0, 50])
ax2.set_xlabel("Precipitación (mm/día)")
ax2.set_ylabel("Frecuencia")
plt.suptitle(
"Distribución de la precipitación "
"en las observaciones de mosquitos"
)
plt.tight_layout()
plt.show()
Correlación#
from scipy.stats import spearmanr
Realizaremos una correlación de rangos de Spearman mediante la biblioteca SciPy de Python.
Esta prueba genera un coeficiente de correlación entre -1 y 1 que representa la fuerza y dirección de la relación. Un valor cercano a 0 indica que no existe una relación monotónica clara. Una correlación positiva significa que, en general, una variable aumenta cuando la otra aumenta. Una correlación negativa significa que una variable tiende a disminuir cuando la otra aumenta.
La correlación de Spearman se diferencia de la correlación de Pearson porque Pearson mide relaciones lineales, mientras que Spearman también puede identificar relaciones monotónicas que no sean lineales.
columns = [
"LarvaeCountProcessed",
"precipitation",
"LST_Day_1km",
"LST_Night_1km",
"NDVI",
"NDWI",
"MeasurementLongitude",
"MeasurementLatitude"
]
correlation_matrix = data_cleaned[
columns
].corr(method="spearman")
sns.heatmap(
correlation_matrix,
annot=True,
cmap="coolwarm",
fmt=".2f",
linewidths=0.5
)
plt.title("Correlaciones de Spearman entre las variables")
plt.tight_layout()
plt.show()
def interpret_spearmanr(variable, alpha=0.05):
"""
Calcula e interpreta la correlación de Spearman entre una
variable ambiental y el conteo procesado de larvas.
"""
correlation_coefficient, p_value = spearmanr(
data_cleaned[variable],
data_cleaned["LarvaeCountProcessed"],
nan_policy="omit"
)
print(
"Coeficiente de correlación de Spearman: "
f"{correlation_coefficient:.3f}"
)
print(f"Valor p: {p_value}")
if p_value < alpha:
if correlation_coefficient > 0:
direction = "positiva"
elif correlation_coefficient < 0:
direction = "negativa"
else:
direction = "nula"
print(
"El resultado es estadísticamente significativo, "
f"con una correlación {direction}."
)
else:
print("El resultado no es estadísticamente significativo.")
print(f"Número de observaciones: {len(data_cleaned)}")
print("\n---NDVI y conteo de larvas de mosquitos---")
interpret_spearmanr("NDVI")
print("\n---NDWI y conteo de larvas de mosquitos---")
interpret_spearmanr("NDWI")
print("\n---Precipitación y conteo de larvas de mosquitos---")
interpret_spearmanr("precipitation")
print(
"\n---Temperatura nocturna y conteo de larvas de mosquitos---"
)
interpret_spearmanr("LST_Night_1km")
print(
"\n---Temperatura diurna y conteo de larvas de mosquitos---"
)
interpret_spearmanr("LST_Day_1km")
A continuación, crearemos gráficos que muestran el ajuste de un modelo de regresión lineal mediante la biblioteca Seaborn de Python.
Puedes pensar en este modelo como una línea que intenta representar la relación entre las variables ambientales y el conteo de larvas de mosquitos.
positive_counts = data_cleaned[
data_cleaned["LarvaeCountProcessed"] > 0
]
fig, (ax1, ax2, ax3) = plt.subplots(
1,
3,
figsize=(15, 4)
)
sns.regplot(
x="NDVI",
y="LarvaeCountProcessed",
data=positive_counts,
line_kws={"color": "darkblue"},
ax=ax1
)
ax1.set_title("Vegetación (NDVI)")
ax1.set_xlabel("NDVI")
ax1.set_ylabel("Conteo de larvas")
sns.regplot(
x="precipitation",
y="LarvaeCountProcessed",
data=positive_counts,
line_kws={"color": "darkblue"},
ax=ax2
)
ax2.set_title("Precipitación")
ax2.set_xlabel("Precipitación (mm/día)")
ax2.set_ylabel("Conteo de larvas")
sns.regplot(
x="LST_Night_1km",
y="LarvaeCountProcessed",
data=positive_counts,
line_kws={"color": "darkblue"},
ax=ax3
)
ax3.set_title("Temperatura nocturna de la superficie terrestre")
ax3.set_xlabel("LST nocturna")
ax3.set_ylabel("Conteo de larvas")
plt.suptitle(
"Condiciones ambientales en las observaciones de mosquitos"
)
plt.tight_layout()
plt.show()