Python

English version

Python

El paquete monitoreo es la base de todas las herramientas: el plugin de QGIS, el servidor MCP y el paquete de R llaman a sus funciones. Todas están en monitoreo.services, devuelven diccionarios y listas de Python (fáciles de pasar a pandas) y pueden escribir CSV, PNG y GeoTIFF.

Instalación y autenticación: Primeros pasos. Las opciones comunes (regiones, periodos, métricas, filtros) se explican en Conceptos.

Los ejemplos de esta página son reales y usan el Parque Nacional de Sierra Nevada (ID de WDPA 389011).

1. Conectar con Earth Engine

from monitoreo import services as mo

mo.init()                     # sin la variable EE_PROJECT: la cuenta compartida de MonitorEO
# mo.init('mi-proyecto-ee')   # o tu cuenta y tu proyecto (o mo.init() con EE_PROJECT)

u = mo.eecu_usage()  # tus EECU de este mes (con tu proyecto: mo.eecu_usage('Contributor') si es su nivel)
print(u['bar'], u['text'])
# [##----------------------------] 8% research EECUs: 8.0% used (2.00 of your 25 EECU-hours of this month with ...

Con la cuenta compartida son tus research EECUs: 25 EECU-h al mes, solo para investigación. Si te queda menos del 10 %, init() da un aviso (EECUWarning) y u['warning'] dice qué hacer; cuando se acaban, init() no conecta con la cuenta compartida hasta el mes siguiente (Conceptos).

2. Elegir la variable

mo.list_variables(sensor='Sentinel-2', text='NDVI')

[{'ref': 27, 'variable': 'NDVI', 'sensor': 'Sentinel-2', 'collection': 'COPERNICUS/S2_SR_HARMONIZED',
  'name': 'Normalized Difference Vegetation Index', 'theme': 'Vegetation', 'group': 'Carbon (Primary Production)',
  'resolution_m': 10, 'dates': '28/03/2017-now', 'unit': '[unitless]', 'rgb': True}]

mo.list_variables() da las 82 entradas del catálogo (cada variable con cada sensor, con el NBR), y mo.list_variables(theme='ecosystems') las de un tema. mo.describe_variable('NDSI', 'Sentinel-2') da su fórmula (N(B3, B11)), escala, filtros por píxel que se le aplican y qué filtros quitarían su señal ('measures': ['snow']: con el NDSI hay que dejar la nieve, quality={'snow': False}). En las dos, rgb dice si las imágenes de la variable tienen color verdadero, como mo.has_rgb('NDVI', 'Modis') (False). Las capas categóricas (las ecorregiones y los biomas, ECOREGION y BIOME de RESOLVE) tienen categorical y n_classes; describe_variable da sus clases (hasta 50; si son más, classes_note dice que se listen con list_classes):

mo.list_classes('BIOME', lang='es')[11]
# {'code': 12, 'name': 'Bosques, bosques abiertos y matorrales mediterráneos', 'colour': 'fe0000'}
mo.list_classes('ECOREGION', 'RESOLVE', text='Iberian')     # por palabras de su nombre, con su bioma (5)
# [{'code': 792, 'name': 'Iberian conifer forests', 'colour': 'cd8966', 'group': 12,
#   'group_name': 'Mediterranean Forests, Woodlands & Scrub'}, ...]

list_classes no usa Earth Engine: lee la tabla de clases del paquete (ver Conceptos).

3. Elegir la región

mo.find_protected_areas('Sierra Nevada')

[{'SITE_ID': 555722086, 'NAME': 'Sierra Nevada', 'DESIG': 'Special Areas of Conservation (Habitats Directive)', ...},
 ...
 {'SITE_ID': 389011, 'NAME': 'Sierra Nevada', 'DESIG': 'Parque Nacional', 'IUCN_CAT': 'II', 'ISO3': 'ESP',
  'REP_AREA': 858.83, 'STATUS_YR': 1999, ...},
 {'SITE_ID': 349467, 'NAME': 'Sierra Nevada', 'DESIG': 'Parque Natural', ...}, ...]

La región del parque nacional es {'wdpa_id': 389011}. Otras formas de dar la región (ecorregiones, biomas, rectángulos, puntos, GeoJSON, ficheros): Conceptos.

sn = {'wdpa_id': 389011}
zonas = {'file': 'zonas.gpkg', 'name_field': 'nombre'}          # una serie por polígono
punto = {'point': [-3.42, 37.13], 'buffer': 1500}               # círculo de 1,5 km de radio
bioma = {'biome_at': [-3.3, 37.1]}                              # el bioma de RESOLVE de un punto (o {'biome_id': 12})
mediterraneo_es = {'biome_id': 12, 'within': {'country': 'ESP'}}   # solo su parte dentro de otra región

Un bioma entero abarca continentes: dalo dentro de otra región (within, que vale con cualquier región) o calcula con MODIS o el clima y un paso anual. Para los píxeles de un bioma en un país, mejor el país con una condición de where sobre las clases de BIOME (abajo) o las zonas de BIOME (Conceptos, within).

Solo una parte de la región

where deja solo los píxeles cuyas variables cumplen unas condiciones, en todas las funciones con región (Conceptos):

arriba = [{'variable': 'ELEVATION', 'sensor': 'Copernicus', 'min': 1500}]
mo.time_series(sn, 'NDVI', 'Modis', '2015-01-01', '2024-12-31', step='annual', where=arriba)
humedos = [{'variable': 'TWI', 'sensor': 'MERIT Hydro', 'min': 10}]               # los fondos de valle húmedos

Una condición sobre una capa categórica (ECOREGION o BIOME de RESOLVE) deja los píxeles de unas clases (classes) o los quita (exclude), por sus códigos (list_classes), en lugar de min y max. Un píxel se queda si las clases que se dejan cubren la mitad de él o más:

mediterraneo = [{'variable': 'BIOME', 'sensor': 'RESOLVE', 'classes': [12]}]
mo.map_layer({'country': 'ESP'}, 'NDVI', 'Modis', '2023-06-01', '2023-08-31', where=mediterraneo)

Una condición también puede quitar algunos usos del suelo según un mapa de usos de un año (mo.landcover_maps() lista los mapas, sus años y las clases de cada mapa en cada uso; Conceptos):

sin_agua = [{'landcover': 'CORINE', 'year': 2018, 'exclude': ['water', 'built']}]
mo.time_series(sn, 'NDVI', 'Sentinel-2', '2018-01-01', '2024-12-31', step='annual', where=sin_agua)

4. Serie temporal

La media mensual del NDVI de Sentinel-2 en el parque, de 2022 a 2023, con CSV y gráfica:

ts = mo.time_series(sn, 'NDVI', 'Sentinel-2', '2022-01-01', '2023-12-31', step='monthly',
                    metric='mean', stat='mean', csv='ndvi_sierra_nevada.csv', chart='ndvi_sierra_nevada.png')
ts['rows'][:2]

[{'region': 'Sierra Nevada', 'date': '2022-01-01', 'value': 0.3537, 'n_pixels': 346520, 'n_images': 30},
 {'region': 'Sierra Nevada', 'date': '2022-02-01', 'value': 0.4036, 'n_pixels': 298880, 'n_images': 29}]

ts también trae title, unit, scale (el píxel usado: 50 m, porque el parque es grande; ver Conceptos), csv y chart:

Con pandas:

import pandas as pd
df = pd.DataFrame(ts['rows'])
df['date'] = pd.to_datetime(df['date'])
df.pivot(index='date', columns='region', values='value')      # una columna por región

Variantes:

mo.time_series(sn, 'LST', 'Modis', '2010-01-01', '2023-12-31', step='annual', metric='max')   # máximo anual
mo.time_series(sn, 'NDVI', 'Landsat', '1985-01-01', '2023-12-31', step='annual',
               season=['06-01', '09-30'])                                           # solo los veranos
mo.time_series(sn, 'NDSI', 'Modis', '2022-10-01', '2023-06-30', step='16-day',
               quality={'snow': False})                                             # la nieve no se quita
mo.time_series(sn, 'SCF', 'Landsat', '1990-10-01', '2023-05-31', step='annual',
               season=['10-01', '05-31'])                      # frecuencia de nieve de cada temporada (con la media)
mo.time_series(punto, 'NDVI', 'Sentinel-2', '2023-01-01', '2023-12-31', step='daily')  # cada imagen

Por los días del año

Una línea por año sobre los días del año, para comparar los años (ver Conceptos):

mo.time_series({'wdpa_id': 389011}, 'NDSI', 'Modis', '2019-01-01', '2024-12-31', step='16-day', doy=True,
               chart='snow.png')

O el perfil estacional, una línea por variable: con interannual, la métrica entre años de cada paso del año (mean, median, min, max, std, cv o pN), y con other_variable (y other_sensor, por defecto sensor), una segunda variable del catálogo en el eje derecho (ver Conceptos). La mediana de los años del NDSI y el NDVI mensuales de Sentinel-2 de 2018 a 2024, cada eje del mínimo al máximo de su variable (lo que se usa por defecto):

pr = mo.time_series(sn, 'NDSI', 'Sentinel-2', '2018-01-01', '2024-12-31', step='monthly', doy=True,
                    interannual='median', other_variable='NDVI',
                    quality={'snow': 'auto'},          # la nieve se queda en el NDSI, que la mide, no en el NDVI
                    csv='ndsi_ndvi.csv', chart='ndsi_ndvi.png')
pr['variables']                 # ['NDSI Sentinel-2', 'NDVI Sentinel-2']: el eje izquierdo y el derecho
pr['ylim'], pr['other_ylim']    # (None, None): cada eje según los valores de su variable
pd.DataFrame(pr['profile']).pivot(index=['doy', 'label'], columns='variable', values='value')  # una columna cada una

pr['profile'] es lo que muestran el gráfico y el CSV: una fila por región, variable y paso del año, con region, variable, part (el mes 1-12, el periodo de 16 días 0, 1..., o el día MM-DD), doy (dónde se dibuja), label (Jan...), value y n_years (los años con valor). pr trae también units y scales (de cada variable, por su nombre), ylabel y other_ylabel, ylim y other_ylim (los ejes dados; None: según los datos), interannual, y rows, los valores de cada año de las dos variables, cada uno con su variable. Los ejes se fijan con ylim=(0, 0.8) y other_ylim=(-0.2, 1); con una sola variable (interannual sin other_variable), hay un solo eje.

5. Mapa

La mediana del NDVI del verano de 2023, como teselas de Earth Engine y una vista previa en PNG:

capa = mo.map_layer(sn, 'NDVI', 'Sentinel-2', '2023-06-01', '2023-08-31', metric='median',
                    thumbnail='ndvi_verano_2023.png')
capa['url']       # 'https://earthengine.googleapis.com/v1/projects/.../maps/.../tiles/{z}/{x}/{y}'
capa['min'], capa['max'], capa['bounds']   # rampa (percentiles 5 y 95) y extensión [oeste, sur, este, norte]

La rampa de colores es la de MonitorEO para el NDVI: del azul (agua, roca o nieve, NDVI cercano a 0) al verde oscuro (vegetación densa).

La elevación, la pendiente, el índice topográfico de humedad y las variables de WorldClim son capas fijas, sin fechas: van sin periodo ni métrica, p. ej. mo.map_layer(sn, 'SLOPE', 'Copernicus') o mo.map_layer(sn, 'TWI', 'MERIT Hydro') (ver Conceptos y TWI).

Las ecorregiones y los biomas (ECOREGION y BIOME, sensor RESOLVE) son capas fijas categóricas: sus valores son códigos de clases. Su mapa tiene el color de cada clase y, en lugar de la rampa (min, max y palette son None), la leyenda de las clases de la región (Conceptos):

b = mo.map_layer({'country': 'ESP'}, 'BIOME', 'RESOLVE')
b['classes']   # [..., {'value': 12, 'label': '12 Mediterranean Forests, Woodlands & Scrub', 'color': '#fe0000'}, ...]

La URL vale unas horas en cualquier visor de teselas XYZ. Con folium:

import folium
m = folium.Map(location=[37.07, -3.1], zoom_start=10)
folium.TileLayer(capa['url'], attr='Google Earth Engine', name=capa['name']).add_to(m)
m

Para trabajar con la imagen de Earth Engine (geemap, más cálculos con ee), composite_image la devuelve con la región:

imagen, regiones = mo.composite_image(sn, 'NDVI', 'Sentinel-2', '2023-06-01', '2023-08-31', metric='median')
# imagen es un ee.Image con la banda 'value'

Agregaciones entre años (ver Conceptos):

mo.map_layer(sn, 'LST', 'Modis', '2014-01-01', '2023-12-31', metric='max', temporal='annual', interannual='median')

Color verdadero de las imágenes

rgb_map toma los mismos argumentos que map_layer y da el color verdadero de las mismas imágenes y píxeles: su reflectancia de superficie roja, verde y azul, con la misma métrica, con las tres bandas estiradas juntas entre el menor percentil 2 y el mayor percentil 98 de las bandas en la región (o vmin y vmax). Ayuda a leer el mapa de un índice: zonas quemadas, suelo desnudo, nieve, agua o nubes que dejaron los filtros (ver Conceptos).

rgb = mo.rgb_map(sn, 'NDVI', 'Sentinel-2', '2023-06-01', '2023-08-31', metric='median',
                 thumbnail='rgb_verano_2023.png')
rgb['bands'], rgb['min'], rgb['max']   # ['red', 'green', 'blue'] y el rango de reflectancia de las tres bandas

Con folium, añade su TileLayer antes que el del NDVI para tener el color verdadero debajo. En GeoTIFF, download_rgb toma los argumentos de download_image y escribe tres bandas float32, red, green y blue (NoData −9999), en la misma rejilla que el GeoTIFF de la variable:

r = mo.download_rgb(punto, 'NDVI', 'Sentinel-2', '2023-06-01', '2023-08-31', path='ndvi_verano_rgb.tif')
r['bands'], r['min'], r['max']         # su 'min' y 'max': el rango para ver las tres bandas juntas

Solo lo tienen las variables de Sentinel-2 y Landsat y la frecuencia de nubes de MODIS (mo.has_rgb(variable, sensor)), con las métricas median, mean, min, max o pN (con max y min, los colores de la imagen en la que la variable es más alta o más baja en cada píxel); las demás se rechazan con un ValueError antes de calcular nada.

Animación (timelapse)

Un fotograma por año (o mes, 16 días o día con imágenes) con su fecha, los mismos colores en todos y los píxeles sin datos en gris claro (ver Conceptos):

a = mo.timelapse({'wdpa_id': 389011}, 'NDVI', 'Modis', '2001-01-01', '2024-12-31', step='annual', path='ndvi.gif')
a['frames'][:2]   # [{'date': '2001-01-01', 'end': '2002-01-01', 'label': '2001', 'images': 23}, ...]
a['min'], a['max'], a['readme']
# un vídeo MP4 (con ffmpeg), dos fotogramas por segundo, un fotograma por mes de los veranos
mo.timelapse({'wdpa_id': 389011}, 'NDSI', 'Sentinel-2', '2019-01-01', '2024-12-31', step='monthly',
             season=['06-01', '09-30'], fps=2, path='ndsi.mp4')

Figura para un informe o un artículo

mo.map_figure({'wdpa_id': 389011}, 'NDVI', 'Modis', '2023-06-01', '2023-08-31', path='ndvi.pdf')  # or .png, .svg

El mapa con su retícula, barra de escala, flecha del norte, escala de colores y créditos (ver Conceptos). Una capa categórica (ECOREGION, BIOME) lleva la leyenda de sus clases en la región en lugar de la escala de colores (como mucho 11, y después cuántas más), y la respuesta trae classes en lugar de min, max y palette.

6. Estadísticas de un mapa por polígono

mo.region_stats(sn, 'LST', 'Modis', '2023-06-01', '2023-08-31', metric='mean')

{'rows': [{'region': 'Sierra Nevada', 'mean': 32.21, 'std': 1.98, 'cv': 6.14, 'min': 27.27, 'p5': 29.41,
           'p50': 32.09, 'p95': 34.90, 'max': 40.86, 'n_pixels': 1051}],
 'stats': ['mean', 'std', 'cv', 'min', 'p5', 'p50', 'p95', 'max'],
 'scale': 1000.0, 'unit': '[ºC]', 'title': 'LST Modis | 2023-06-01 / 2023-08-31 (mean)'}

stats elige los estadísticos: cualquiera de mean, median, min, max, std (desviación típica), cv (coeficiente de variación, 100 × desviación típica / media, en %) y percentiles pN, para cualquier variable, también las capas fijas (mo.region_stats(sn, 'TWI', 'MERIT Hydro')); por defecto, mean, std, cv, min, p5, p50, p95 y max. Una capa categórica (ECOREGION, BIOME) solo admite mode (la clase mayoritaria, por defecto) y count: mo.region_stats(sn, 'BIOME', 'RESOLVE') da en cada fila mode (el código de la clase) y mode_label (su código y su nombre); la superficie de cada clase es area_by_class. Con temporal e interannual, el mapa es el de un valor entre años, p. ej. cuánto cambia la precipitación anual de un año a otro (ver Conceptos):

mo.region_stats(sn, 'PRECIP', 'TerraClimate', '2001-01-01', '2020-12-31', metric='sum', temporal='annual',
                interannual='cv', stats=['mean', 'median', 'p90'])   # en %

Superficie por clase

Las hectáreas y el % de cada clase (ver Conceptos): de un mapa de usos del suelo de un año, del mapa de una variable cortado en unos valores, de su tendencia en cada píxel, o de las clases de una capa categórica:

a = mo.area_by_class({'wdpa_id': 389011}, landcover='CORINE', year=2018, csv='corine.csv', chart='corine.png')
next(r for r in a['rows'] if r['class'] == 311)  # su fila (rows: cada clase del mapa, por código, también con 0 ha):
# {'region': ..., 'map': 'CORINE Land Cover', 'year': 2018, 'class': 311, 'label': '311 Broad-leaved forest', 'landcover': 'trees', 'ha': ..., 'percent': ...}
a['by_landcover']   # las hectáreas de cada uso del suelo general en cada polígono
# sus cambios de 2006 a 2018: cada par de clases, la superficie que cambió y un diagrama de Sankey
c = mo.area_by_class({'wdpa_id': 389011}, landcover='CORINE', year=2006, other_year=2018, chart='sankey.png')
c['changed_ha'], c['unchanged_ha'], c['by_landcover'][:3]
mo.area_by_class({'wdpa_id': 389011}, 'ELEVATION', 'Copernicus', breaks=[1500, 2000, 2500, 3000])
mo.area_by_class({'wdpa_id': 389011}, 'NDVI', 'Modis', '2001-01-01', '2024-12-31', trend=True)
# las ecorregiones del parque (sin breaks ni trend): mode 'classes', cada una con su bioma (group, group_label)
e = mo.area_by_class({'wdpa_id': 389011}, 'ECOREGION', 'RESOLVE', csv='ecorregiones.csv', chart='ecorregiones.png')
e['by_group']       # las hectáreas de cada bioma en cada polígono

Las clases de una capa categórica son las suyas, las que tienen superficie en la región, con sus colores; con ella no caben breaks ni trend.

Gradientes: franjas de altitud y perfiles

# el NDVI de cada franja de altitud de 500 m (ver Conceptos, Gradientes)
z = mo.region_stats({'wdpa_id': 389011}, 'NDVI', 'Modis', '2023-06-01', '2023-08-31', stats=['p25', 'p50', 'p75'],
                    zones={'variable': 'ELEVATION', 'sensor': 'Copernicus', 'breaks': [1500, 2000, 2500, 3000]},
                    csv='bands.csv', chart='bands.png')
# o las franjas del índice topográfico de humedad: {'variable': 'TWI', 'sensor': 'MERIT Hydro', 'breaks': [8, 12]}
# el NDVI de cada bioma de España: las clases de una capa categórica son las zonas, sin breaks
b = mo.region_stats({'country': 'ESP'}, 'NDVI', 'Modis', '2023-06-01', '2023-08-31', stats=['p25', 'p50', 'p75'],
                    zones={'variable': 'BIOME', 'sensor': 'RESOLVE'}, csv='biomas.csv')
# la LST a lo largo de una línea que sube a Sierra Nevada, con la elevación
p = mo.transect([[-3.60, 37.15], [-3.45, 37.10], [-3.39, 37.06]], 'LST', 'Modis', '2023-07-01', '2023-08-31',
                chart='profile.png')
p['rows'][0]   # {'distance_m': 500.0, 'lon': ..., 'lat': ..., 'value': 43.5, 'elevation': 684.0}

Con las zonas de una capa categórica, cada clase con píxeles en cada polígono es una fila (zone, su código, y zone_label, su código y su nombre; zone_from y zone_to vacíos) y, con chart, el gráfico de cada clase; 'classes': [12, 13] deja solo esas. Una capa categórica no tiene perfil a lo largo de una línea (transect da un ValueError).

7. Descargar un mapa en GeoTIFF

El percentil 95 del NDVI de 2023 en un círculo de 1,5 km, a 10 m:

r = mo.download_image(punto, 'NDVI', 'Sentinel-2', '2023-01-01', '2023-12-31', metric='p95', path='ndvi_p95.tif')

{'file': 'ndvi_p95.tif', 'grid': {'crs': 'EPSG:32630', 'resolution': 10.0, 'width': 300, 'height': 301, ...},
 'resolution': 10.0, 'native_resolution': 10.0, 'images': 74, 'unit': '[unitless]', 'min': 0.31, 'max': 1.0, ...}

El GeoTIFF es float32 con NoData −9999, en la zona UTM del centro de la región (crs='auto'; en una región de más de 40° de longitud, como un bioma entero, EPSG:6933, de igual área: ver Conceptos) u otro sistema (crs='EPSG:25830'), con el píxel de los datos o el que se pida (resolution=20). Earth Engine lo calcula por teselas en paralelo: para regiones de unos millones de píxeles. Con max_pixels=4_000_000, una región mayor se descarga con píxeles 2, 3... veces más grandes; a resolución completa, mejor exportarlo a Google Drive. El color verdadero de las mismas imágenes, en la misma rejilla: download_rgb (Color verdadero de las imágenes).

El GeoTIFF de una capa fija va sin periodo ni métrica: mo.download_image(sn, 'TWI', 'MERIT Hydro', path='twi.tif'). El de una capa categórica guarda los códigos de las clases (en píxeles más del doble de grandes que su píxel nominal, 1 km, la clase mayoritaria); la respuesta trae classes, la leyenda de las clases del fichero, en lugar de min, max y palette, y su README y sus metadatos (classes) dicen el nombre de cada código:

e = mo.download_image(sn, 'ECOREGION', 'RESOLVE', path='ecorregiones.tif')
e['classes']   # las clases del fichero: [{'value' (su código), 'label' (código y nombre), 'color'}, ...]

8. Tendencias estacionales

¿Ha cambiado el NDVI de cada mes entre 2005 y 2024? Con MODIS (20 años):

t = mo.trend_table(sn, 'NDVI', 'Modis', '2005-01-01', '2024-12-31', step='monthly',
                   csv='tendencias_ndvi.csv', chart='tendencias_ndvi.png')
[(r['label'], round(r['slope'], 4), round(r['p_value'], 3), r['significant']) for r in t['rows']]

[('January', -0.0003, 0.697, False), ('February', 0.0007, 0.65, False), ('March', 0.0015, 0.119, False),
 ('April', 0.0004, 0.795, False), ('May', 0.0009, 0.516, False), ('June', 0.0012, 0.243, False),
 ('July', 0.0013, 0.27, False), ('August', 0.0019, 0.005, True), ('September', 0.002, 0.023, True),
 ('October', 0.001, 0.399, False), ('November', 0.0011, 0.436, False), ('December', 0.0005, 0.516, False)]

El NDVI medio del parque sube en todos los meses salvo enero, pero la subida solo es significativa (p < 0,05) a finales del verano: unas +0,002 unidades por año en agosto y septiembre. Cada fila trae también n (años con datos), tau, z e intercept; t['yearly'] tiene los valores anuales de cada mes. La gráfica muestra esos valores y la recta de Sen de cada mes:

En cada píxel. El mapa de la pendiente de Sen del NDVI anual, solo donde la tendencia es significativa:

mo.trend_map(sn, 'NDVI', 'Modis', '2005-01-01', '2024-12-31', step='annual', output='slope', alpha=0.05,
             thumbnail='tendencia_ndvi_anual.png')

Verde: aumento; rojo: descenso; en blanco, los píxeles sin tendencia significativa. Con valores, en GeoTIFF:

d = mo.download_trend(sn, 'NDVI', 'Modis', '2005-01-01', '2024-12-31', step='annual', path='tendencia_ndvi_anual.tif')
d['bands']    # ['slope', 'tau', 'z', 'p_value', 'intercept', 'n']
d['unit']     # '[unitless] / year'

Para un mes o un bloque de 16 días: step='monthly', part=8 (agosto) o step='16-day', part=14 (días 209-224). output elige la banda del mapa de teselas: 'slope', 'tau', 'z', 'p_value', 'intercept' o 'n'.

9. Comparar dos periodos

El NDVI de cada mes en 2020-2024 frente a 2005-2009:

c = mo.compare_series(sn, 'NDVI', 'Modis', '2005-01-01', '2009-12-31',     # periodo de referencia
                      '2020-01-01', '2024-12-31',                          # el otro periodo
                      step='monthly', csv='comparacion_ndvi.csv', chart='comparacion_ndvi.png')
c['profile'][:2]

[{'region': 'Sierra Nevada', 'part': 1, 'key': 'm01', 'label': 'January', 'reference': 0.4954, 'other': 0.4953,
  'difference': -0.0002, 'n_reference': 5, 'n_other': 5},
 {'region': 'Sierra Nevada', 'part': 2, 'key': 'm02', 'label': 'February', 'reference': 0.4681, 'other': 0.4829,
  'difference': 0.0149, 'n_reference': 5, 'n_other': 5}]

c['rows'] tiene los valores de cada año de los dos periodos. El mapa de la diferencia (la mediana de 2020-2024 menos la de 2005-2009):

mo.compare_map(sn, 'NDVI', 'Modis', '2005-01-01', '2009-12-31', '2020-01-01', '2024-12-31',
               metric='median', thumbnail='diferencia_ndvi.png')
mo.download_comparison(sn, 'NDVI', 'Modis', '2005-01-01', '2009-12-31', '2020-01-01', '2024-12-31',
                       metric='median', path='diferencia_ndvi.tif')     # bandas difference, reference, other

Severidad de un incendio

El incendio de agosto de 2023 en el Parque Nacional del Teide, con Sentinel-2 (ver Conceptos):

s = mo.burn_severity({'wdpa_id': 389013}, 'Sentinel-2', '2023-05-15', '2023-08-14',   # antes
                     '2023-09-15', '2023-11-30', path='severidad.tif', csv='severidad.csv')   # después
s['burned_ha']                                  # {'El Teide': 956.4}: hectáreas con dNBR >= 0,1
[(r['label'], round(r['ha'])) for r in s['classes'] if r['class'] >= 1]
# [('Low severity', 559), ('Moderate-low severity', 318), ('Moderate-high severity', 79), ('High severity', 0)]
r = mo.burn_severity({'wdpa_id': 389013}, 'Sentinel-2', '2023-05-15', '2023-08-14',
                     '2023-09-15', '2023-11-30', index='RBR', path='rbr.tif')   # el RBR, continuo
r['stats']                                      # [{'region', 'mean', 'median', 'p10', 'p90', 'min', 'max', 'ha'}]
d = mo.burn_severity({'wdpa_id': 389013}, 'Sentinel-2', '2023-05-15', '2023-08-14',
                     '2024-05-15', '2024-08-14', index='RdNBR', path='rdnbr.tif')   # un año después
d['classes']                                    # [{'region', 'class', 'label', 'rdnbr_min', 'rdnbr_max', 'ha'}]
d['burned_ha']                                  # hectáreas que cambiaron: RdNBR >= 0,069
d['stats']                                      # media, mediana, p10, p90, mín. y máx. del RdNBR
m = mo.burn_severity({'wdpa_id': 389013}, 'Modis', '2023-05-15', '2023-08-14',   # MODIS Aqua, 500 m
                     '2024-05-15', '2024-08-14', index='RdNBR', collection='MODIS/061/MYD09GA')

Con index='RBR' da el RBR en lugar del dNBR y sus clases (ver Conceptos). Con index='RdNBR', el dNBR relativo de Miller y Thode (2007) y sus clases de severidad (class 0 a 3: Unchanged, Low severity, Moderate severity y High severity), la superficie que cambió (burned_ha, RdNBR ≥ 0,069), las estadísticas del RdNBR (stats) y un GeoTIFF con las bandas rdnbr, severity, dnbr, nbr_before y nbr_after (ver Conceptos). Sus clases se calibraron con el periodo después un año tras el incendio, en la misma época que el periodo antes. Solo el RdNBR se calcula también con MODIS (sensor='Modis', 500 m, clases aproximadas): Terra (MODIS/061/MOD09GA) por defecto, o Aqua con collection='MODIS/061/MYD09GA'.

Con perimeter, también el perímetro de la superficie quemada en un fichero de polígonos (.gpkg, .geojson o .shp), para usarlo como región de otro análisis (ver Conceptos):

s = mo.burn_severity({'wdpa_id': 389013}, 'Sentinel-2', '2023-05-15', '2023-08-14', '2023-09-15', '2023-11-30',
                     perimeter='perimetro.gpkg')  # dNBR >= 0,1, sin parches ni islas de menos de 1 ha
s['perimeter_ha']                                 # hectáreas de cada polígono de la región
quemado = {'file': 'perimetro.gpkg', 'name_field': 'region'}
ndvi = mo.time_series(quemado, 'NDVI', 'Sentinel-2', '2023-09-01', '2025-08-31')  # la recuperación, dentro

burned_min cambia el valor desde el que un píxel está quemado (obligatorio con index='RBR', que no tiene clases) y min_ha, el parche y la isla más pequeños (min_ha=0: todos los píxeles).

10. Leer los GeoTIFF

Con GDAL (viene en el entorno de conda) o rasterio:

from osgeo import gdal
ds = gdal.Open('tendencia_ndvi_anual.tif')
nombres = [ds.GetRasterBand(i + 1).GetDescription() for i in range(ds.RasterCount)]   # slope, tau, ...
pendiente = ds.GetRasterBand(1).ReadAsArray()        # NoData: -9999

11. Exportar a Google Drive

Para mapas grandes a resolución completa o series muy largas. Earth Engine los calcula sin límite de tiempo y los deja en una carpeta de tu Google Drive:

tarea = mo.export_image(sn, 'NDVI', 'Sentinel-2', '2023-01-01', '2023-12-31', metric='median',
                        resolution=10, folder='MonitorEO')
mo.task_status(tarea['operation'])     # {'state': 'RUNNING', 'progress': 0.4, ...} hasta 'SUCCEEDED'

mo.time_series_to_drive(sn, 'NDVI', 'Sentinel-2', '2017-04-01', '2023-12-31', step='daily', folder='MonitorEO')

Las tareas también se ven en https://code.earthengine.google.com/tasks.

12. Varias variables para un modelo: recetas

download_variables descarga varias variables y métricas en rejillas alineadas (mismo CRS, mismos bordes, NoData −9999) sobre un área de estudio y, con un píxel más grueso, sobre las ecorregiones que toca. Es lo que usa eHabitat+ para preparar sus variables de entrada. La receta es un YAML (o un diccionario):

crs: auto                                   # UTM del centro del área (EPSG:6933 si abarca más de 40° de longitud), o p. ej. EPSG:25830
study: {resolution: 10, buffer: 500}        # área de estudio: rectángulo + 500 m, descarga directa
ecoregion: {resolution: 100, drive: true}   # opcional: ecorregiones que toca, vía Google Drive
period: ['2021-01-01', '2023-12-31']
quality: {clouds: true, water: true, snow: true, aerosols: true, cloud_score: 0.6}
variables:
  - {index: NDVI, sensor: Sentinel-2, metric: [p95, p5]}                       # un fichero por métrica
  - {index: LST, sensor: Modis, temporal: annual, metric: max, interannual: median}
  - {name: pendiente, image: COPERNICUS/DEM/GLO30_2024_1, band: DEM, terrain: slope}  # cualquier imagen de Earth Engine

Primero, un ensayo (dry_run) para ver el tamaño de las rejillas y el número de imágenes:

p = mo.download_variables('receta.yaml', {'wdpa_id': 389011}, 'variables', dry_run=True)
p['grids']['study']       # {'crs': 'EPSG:32630', 'resolution': 10.0, 'width': 7962, 'height': 3102,
                          #  'pixels': 24698124, 'mb_per_variable': 94.22, ...}
[f['file'] for f in p['files']]
# ['1_NDVI_S2_p95.tif', '2_NDVI_S2_p5.tif', '3_LST_MODIS_annual-max_median.tif', '4_pendiente.tif']

Y después la descarga:

mo.download_variables('receta.yaml', 'area_estudio.gpkg', 'variables', log=print)

• Las variables del área de estudio quedan en variables/ y las de las ecorregiones en variables/ecoregion/, con sus polígonos en ecoregions.gpkg.
• monitoreo.json guarda con qué se hizo cada fichero: al repetirlo solo se descarga lo que haya cambiado (overwrite=True lo descarga todo).
• Las rejillas con drive: true se exportan a Google Drive (sin límite de tamaño) y se importan desde Google Drive para escritorio al terminar. Si se interrumpe, al repetirlo se esperan las exportaciones en curso en vez de lanzarlas otra vez.

13. Otros datos de Earth Engine, por su ID

Todas las funciones que calculan aceptan, en lugar de variable y sensor, cualquier imagen o colección de imágenes de Earth Engine por su ID: collection (imágenes con fechas) o image (una imagen fija, o una colección de teselas sin fechas), con band (una, o varias combinadas en cada píxel con band_reducer) y, si hace falta, band_scale, band_offset y unit (en UDUNITS); con image y una banda, también terrain. Qué son, qué análisis tiene cada uno y cómo se citan: Conceptos.

Primero, qué es el ID (lee sus metadatos en Earth Engine, sin calcular):

d = mo.describe_asset('IDAHO_EPSCOR/TERRACLIMATE')
d['type'], d['use'], d['cadence']      # ('image_collection', 'collection', '1 month')
d['dates']                             # ['1958-01-01', '2024-12-01'], como dice su ficha
next(b for b in d['bands'] if b['band'] == 'pr')
# {'band': 'pr', 'description': 'Precipitation accumulation', 'units': 'mm', 'unit': 'mm', 'resolution_m': 4638.3, ...}
d['arguments']      # los que calculan su primera banda: {'collection': 'IDAHO_EPSCOR/TERRACLIMATE', 'band': ..., ...}
d['references']     # cómo citarlo (Abatzoglou et al., 2018, con su DOI); d['credits'], los que piden sus condiciones

La precipitación de cada mes de TerraClimate y la temperatura media anual de WorldClim (1960-1990), una imagen fija guardada en décimas de grado:

pr = dict(collection='IDAHO_EPSCOR/TERRACLIMATE', band='pr', unit='mm month-1')    # mm de cada mes
ts = mo.time_series(sn, **pr, start='1991-01-01', end='2020-12-31', step='annual', metric='sum',
                    csv='precip_terraclimate.csv', chart='precip_terraclimate.png')
ts['unit']          # 'mm': la suma de los 12 meses de cada año es el total anual
mo.trend_table(sn, **pr, start='1991-01-01', end='2020-12-31', step='annual', metric='sum')   # ¿ha cambiado?

bio01 = dict(image='WORLDCLIM/V1/BIO', band='bio01', band_scale=0.1, unit='degC')
mo.region_stats(sn, **bio01)                       # media, mínimo, percentiles y máximo, en degC
mo.download_image(sn, **bio01, path='bio01.tif')   # sin periodo ni métrica: una imagen fija

Varias bandas en un valor. band también es una lista de bandas (['aet', 'def']) o un texto con comas ('aet, def'), con band_reducer: 'sum' (su suma) o 'mean' (su media) en cada píxel. Cada banda se multiplica por band_scale y se le suma band_offset, y después se combinan, solo donde todas tienen datos, en los píxeles de la primera banda. La unidad es la de cada banda, y por tanto la de su suma o su media: una suma de bandas no pasa de una tasa a un total (la métrica temporal sum sí). La evapotranspiración real más el déficit hídrico climático de TerraClimate es la evapotranspiración potencial; y la media de la temperatura máxima del mes más cálido y la mínima del mes más frío de WorldClim (Conceptos):

pet = dict(collection='IDAHO_EPSCOR/TERRACLIMATE', band=['aet', 'def'], band_reducer='sum', band_scale=0.1,
           unit='mm month-1')
ts = mo.time_series(sn, **pet, start='1991-01-01', end='2020-12-31', step='annual', metric='sum')
ts['title'], ts['unit']   # 'IDAHO_EPSCOR/TERRACLIMATE sum of aet, def ...', 'mm': el total de cada año

bio = dict(image='WORLDCLIM/V1/BIO', band=['bio05', 'bio06'], band_reducer='mean', band_scale=0.1, unit='degC')
mo.region_stats(sn, **bio)

band_reducer hace falta con dos bandas o más y no se admite con una; terrain solo va con una banda. Los resultados se llaman por el ID y la combinación (IDAHO_EPSCOR/TERRACLIMATE sum of aet, def), sus ficheros por el final del ID, la combinación y las bandas (TERRACLIMATE_sum_aet_def; con más de tres bandas, la combinación, cuántas son y la primera y la última, como sum12_jan-dec), y el README y el informe dicen the sum of the bands aet and def on each pixel; values × 0.1 + 0 (each band, before the sum).

Los conjuntos de datos del catálogo de la comunidad de GEE, por sus temas (Conceptos), de la copia de su índice que trae el paquete (sin Earth Engine):

mo.community_themes('es')[2]
# {'theme': 'Elevation and Bathymetry', 'name': 'Elevación y batimetría', 'datasets': 15}
[d['id'] for d in mo.community_datasets(theme='Elevación y batimetría')]     # o 'Elevation and Bathymetry'
# ['projects/sat-io/open-datasets/ASTER/GDEM', ..., 'projects/sat-io/open-datasets/FABDEM', ...]
mo.community_datasets(text='fabdem')
# [{'id': 'projects/sat-io/open-datasets/FABDEM', 'title': 'FABDEM (Forest And Buildings removed Copernicus 30m DEM)',
#   'type': 'image_collection', 'theme': 'Elevation and Bathymetry', 'provider': 'Hawker, L., Uhe, P., ...',
#   'license': 'Creative Commons Attribution Non Commercial Share Alike 4.0 International', 'tags': '...',
#   'docs': 'https://gee-community-catalog.org/projects/fabdem/'}]
mo.describe_asset('projects/sat-io/open-datasets/FABDEM')['use']   # 'image': teselas sin fechas, en un mosaico

fabdem = dict(image='projects/sat-io/open-datasets/FABDEM', band='b1')
mo.region_stats(sn, **fabdem, unit='m')            # la elevación del suelo, sin árboles ni edificios
mo.download_image(sn, **fabdem, terrain='slope', path='fabdem_pendiente.tif')   # su pendiente, en grados (degree)

• Los errores dan un ValueError que dice qué hacer: un ID que esta cuenta no puede leer, un ID dado como variable, una imagen o una colección de teselas sin fechas dada como collection, una serie temporal, una tendencia o una comparación de una image, sin band, varias bandas sin band_reducer (o band_reducer con una sola), una banda repetida, terrain con varias bandas, una unit que no está en UDUNITS, o una unit distinta de degree con terrain (Problemas frecuentes).
• rgb_map, download_rgb y burn_severity son solo para las variables del catálogo.
• En una receta de download_variables se indican con las mismas claves, con scale y offset en lugar de band_scale y band_offset: {collection: IDAHO_EPSCOR/TERRACLIMATE, band: pr, unit: mm month-1, temporal: annual, metric: sum}; varias bandas, en una lista con band_reducer en esa variable: {name: bio_temp, image: WORLDCLIM/V1/BIO, band: [bio05, bio06], band_reducer: mean, scale: 0.1, unit: degC}. Una unit vacía (unit: '') es una unidad no indicada, como en las funciones, no una adimensional (unit: '1').

Busca en el catálogo de datos de Earth Engine por palabras en inglés, fechas y zona (ver Conceptos):

mo.find_datasets('soil moisture', start='2020-01-01', end='2020-12-31', bbox=[-3.6, 36.9, -2.6, 37.3], limit=5)
# [{'id': 'NASA/SMAP/SPL3SMP_E/005', 'title': ..., 'type': 'image_collection', 'start': ..., 'gsd': ..., 'page': ...}, ...]

14. En un cuaderno, sobre los mapas de geemap

En Jupyter o Google Colab, monitoreo.notebook pone los mapas de MonitorEO sobre los mapas interactivos de geemap (Wu, 2020), con la escala de colores de la variable y su unidad, el contorno de la región y los créditos de los datos. Instala MonitorEO con el extra notebook (geemap) o, si ya tienes MonitorEO, solo geemap (pip install geemap):

pip install "monitoreo[notebook] @ https://github.com/monitoreo-obsnev/monitoreo/archive/refs/heads/main.zip"

import monitoreo.services as mo
from monitoreo import notebook
mo.init()

m = notebook.Map()
cite = notebook.add_map(m, {'wdpa_id': 389011}, 'NDVI', 'Sentinel-2', '2023-06-01', '2023-08-31')
m  # el mapa, con su barra de colores, el contorno del parque y los créditos

# los dos periodos lado a lado (un mapa partido), con los mismos colores
m2, cite2 = notebook.compare({'wdpa_id': 389011}, 'NDVI', 'Modis', '2005-06-01', '2005-08-31', '2023-06-01', '2023-08-31')
m2

# el mapa de cada año con un deslizador de tiempo
m3 = notebook.Map()
notebook.time_slider(m3, {'wdpa_id': 389011}, 'NDVI', 'Modis', '2019-01-01', '2023-12-31', step='annual')
m3

Cada función toma los argumentos de map_layer (filtros de píxel, where, datos de Earth Engine por su ID...) y devuelve lo que hay que citar con el mapa (references, credits). El mapa de una capa categórica (ECOREGION, BIOME) lleva, en lugar de la barra de colores, la leyenda de sus clases en la región (y add_map devuelve sus classes). mo.composite_image() da la imagen de Earth Engine de un mapa, para otros usos con geemap o el Code Editor.

Wu, Q. (2020). geemap: A Python package for interactive mapping with Google Earth Engine. Journal of Open Source Software, 5(51), 2305. https://doi.org/10.21105/joss.02305

Referencia

Función	Qué hace	Devuelve
init(project, shared)	Conecta con Earth Engine: con tu proyecto o, sin él (o con shared=True), con la cuenta compartida de MonitorEO, que no exporta a Google Drive	el proyecto
eecu_usage(limit, refresh, lang)	Tus EECU gastados este mes (research EECUs con la cuenta compartida; con tu cuenta, la cuota de tu proyecto)	used_percent, left_percent, used_hours, limit_hours, warning, text, bar
about(lang)	Descripción, agradecimientos, logotipos y enlaces	dict
list_variables(sensor, text, lang, theme)	Variables del catálogo (theme='ecosystems': las de un tema)	lista
describe_variable(variable, sensor)	Todo sobre una variable (de una capa categórica, también sus clases)	dict
list_classes(variable, sensor, text, limit, lang)	Las clases de una capa categórica (ECOREGION, BIOME de RESOLVE): código, nombre, color y, de las ecorregiones, su bioma (sin Earth Engine)	lista
describe_asset(asset)	Qué es una imagen o colección de Earth Engine, por su ID: imagen o colección, bandas, unidades, escala, fechas, licencia y cómo citarla (lee sus metadatos en Earth Engine)	type, use, bands, arguments, references, ...
community_themes(lang), community_datasets(theme, text, kind, limit)	Temas y conjuntos de datos del catálogo de la comunidad de GEE (sin Earth Engine)	listas
find_protected_areas(text, designation=None)	Áreas protegidas de WDPA por nombre (designation='ramsar': solo sitios Ramsar)	lista
find_admin_units(text, level=1, country=None)	Unidades administrativas de FAO GAUL 2025 por nombre, con su código	lista
time_series(region, variable, sensor, start, end, step, metric, stat, ..., doy, interannual, other_variable, other_sensor, ylim, other_ylim)	Serie temporal; con doy, por los días del año: una línea por año o, con interannual u other_variable, el perfil estacional de una o dos variables, en dos ejes y	rows, csv, chart, scale, unit, title; con el perfil, también profile, variables, units, ylim, other_ylim
map_layer(region, variable, sensor, start, end, metric, temporal, ...)	Teselas del mapa	url, min, max, palette (de una capa categórica, classes, la leyenda de sus clases en la región), bounds, thumbnail (None, y el motivo en thumbnail_error, si Earth Engine no consigue hacer la vista previa PNG en 5 minutos: las teselas son el mapa)
rgb_map(...), has_rgb(variable, sensor)	Teselas del color verdadero de las mismas imágenes y píxeles que map_layer; si una variable lo tiene	url, bands, min, max, bounds, thumbnail; True o False
thumbnail(...)	Vista previa del mapa	PNG (bytes)
composite_image(...)	La imagen de Earth Engine del mapa	(ee.Image, ee.FeatureCollection)
region_stats(..., stats, zones)	Estadísticas del mapa por polígono (stats: cualquiera de mean, median, min, max, std, cv, pN; de una capa categórica, mode y count), o de cada franja o clase de otra capa (zones)	rows, stats, scale, unit
area_by_class(region, variable, sensor, ..., landcover, year, breaks, trend)	Hectáreas y % de cada clase: de un mapa de usos del suelo, de una variable cortada, de su tendencia o de una capa categórica	rows, mode, classes, csv, chart; by_landcover o by_group
download_image(..., path)	Mapa en GeoTIFF	file, grid, resolution, ... (de una capa categórica, classes)
download_rgb(..., path)	Color verdadero en GeoTIFF (bandas red, green, blue), en la rejilla de download_image	file, grid, bands, min, max, ...
trend_table(..., step, alpha)	Tendencias de la región	rows, yearly, csv, chart
trend_map(..., step, part, output, alpha)	Teselas del mapa de tendencias	como map_layer
download_trend(..., step, part, path)	Tendencias por píxel en GeoTIFF	file, bands, ...
compare_series(..., other_start, other_end, step)	Comparación de periodos	profile, rows, csv, chart
compare_map(...) / download_comparison(..., path)	Mapa de la diferencia	como map_layer / download_image
burn_severity(region, sensor, start, end, other_start, other_end, path, csv, index, collection)	Severidad de un incendio (dNBR, clases del USGS; con index='RBR', el RBR; con index='RdNBR', el RdNBR y las clases de Miller y Thode, también con MODIS)	classes, burned_ha (o stats del RBR; el RdNBR, los tres), file
export_image(...) / time_series_to_drive(...)	Exportaciones a Google Drive	operation, file, folder
task_status(operation)	Estado de una exportación	state, progress, error
download_variables(recipe, area, folder, ...)	Receta de variables en GeoTIFF	resumen

Las funciones que calculan aceptan también, en lugar de variable y sensor, collection o image, con band, band_reducer, band_scale, band_offset, unit y terrain (13. Otros datos de Earth Engine). Cada función tiene su documentación completa, también de estos argumentos: help(mo.time_series).