Python package

Versión en español

Python package

The monitoreo package is the base of all the tools: the QGIS plugin, the MCP server and the R package call its functions. They are all in monitoreo.services, return Python dictionaries and lists (easy to pass to pandas) and can write CSV, PNG and GeoTIFF files.

Installation and authentication: Getting started. The common options (regions, periods, metrics, filters) are explained in Concepts.

The examples on this page are real and use the Sierra Nevada National Park (WDPA ID 389011).

1. Connect to Earth Engine

from monitoreo import services as mo

mo.init()                     # without the EE_PROJECT variable: the shared account of MonitorEO
# mo.init('my-ee-project')    # or your account and your project (or mo.init() with EE_PROJECT)

u = mo.eecu_usage()  # your EECUs of this month (with your project: mo.eecu_usage('Contributor') if that is its tier)
print(u['bar'], u['text'])
# [##----------------------------] 8% research EECUs: 8.0% used (2.00 of your 25 EECU-hours of this month with ...

With the shared account they are your research EECUs: 25 EECU-h per month, for research only. If less than 10 % is left, init() gives a warning (EECUWarning) and u['warning'] says what to do; when they run out, init() does not connect to the shared account until the next month (Concepts).

2. Choose the 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', 'group': 'Carbon (Primary Production)',
  'resolution_m': 10, 'dates': '28/03/2017-now', 'unit': '[unitless]'}]

mo.list_variables() gives the 34 of the catalog (with the NBR). mo.describe_variable('NDSI', 'Sentinel-2') gives its formula (N(B3, B11)), scale, the per-pixel filters applied to it and which filters would remove its signal ('measures': ['snow']: with the NDSI the snow must be kept, quality={'snow': False}).

3. Choose the region

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', ...}, ...]

The region of the national park is {'wdpa_id': 389011}. Other ways to give the region (ecoregions, rectangles, points, GeoJSON, files): Concepts.

sn = {'wdpa_id': 389011}
zones = {'file': 'zonas.gpkg', 'name_field': 'nombre'}          # one series per polygon
point = {'point': [-3.42, 37.13], 'buffer': 1500}               # circle of 1.5 km radius

4. Time series

The monthly mean of the Sentinel-2 NDVI in the park, from 2022 to 2023, with CSV and chart:

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 also has title, unit, scale (the pixel used: 50 m, because the park is large; see Concepts), csv and chart:

With 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')      # one column per region

Variants:

mo.time_series(sn, 'LST', 'Modis', '2010-01-01', '2023-12-31', step='annual', metric='max')   # annual maximum
mo.time_series(sn, 'NDVI', 'Landsat', '1985-01-01', '2023-12-31', step='annual',
               season=['06-01', '09-30'])                                           # only the summers
mo.time_series(sn, 'NDSI', 'Modis', '2022-10-01', '2023-06-30', step='16-day',
               quality={'snow': False})                                             # the snow is not removed
mo.time_series(point, 'NDVI', 'Sentinel-2', '2023-01-01', '2023-12-31', step='daily')  # every image

5. Map

The median of the summer 2023 NDVI, as Earth Engine tiles and a PNG preview:

layer = mo.map_layer(sn, 'NDVI', 'Sentinel-2', '2023-06-01', '2023-08-31', metric='median',
                     thumbnail='ndvi_summer_2023.png')
layer['url']       # 'https://earthengine.googleapis.com/v1/projects/.../maps/.../tiles/{z}/{x}/{y}'
layer['min'], layer['max'], layer['bounds']   # ramp (percentiles 5 and 95) and extent [west, south, east, north]

The colour ramp is the one of MonitorEO for the NDVI: from blue (water, rock or snow, NDVI close to 0) to dark green (dense vegetation).

The URL is valid for a few hours in any XYZ tile viewer. With folium:

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

To work with the Earth Engine image (geemap, more computations with ee), composite_image returns it with the region:

image, regions = mo.composite_image(sn, 'NDVI', 'Sentinel-2', '2023-06-01', '2023-08-31', metric='median')
# image is an ee.Image with the band 'value'

Aggregations across years (see Concepts):

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

6. Statistics of a map by polygon

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

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

7. Download a map as GeoTIFF

The 95th percentile of the 2023 NDVI in a circle of 1.5 km, at 10 m:

r = mo.download_image(point, '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, ...}

The GeoTIFF is float32 with NoData −9999, in UTM (crs='auto') or another system (crs='EPSG:25830'), with the pixel of the data or the one requested (resolution=20). Earth Engine computes it in tiles, in parallel: for regions of a few million pixels. With max_pixels=4_000_000, a larger region is downloaded with pixels 2, 3... times larger; at full resolution, it is better to export it to Google Drive.

8. Seasonal trends

Has the NDVI of each month changed between 2005 and 2024? With MODIS (20 years):

t = mo.trend_table(sn, 'NDVI', 'Modis', '2005-01-01', '2024-12-31', step='monthly',
                   csv='trends_ndvi.csv', chart='trends_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)]

The mean NDVI of the park rises in every month except January, but the rise is only significant (p < 0.05) in the late summer: about +0.002 units per year in August and September. Each row also has n (years with data), tau, z and intercept; t['yearly'] has the annual values of each month. The chart shows those values and the Sen line of each month:

In each pixel. The map of Sen's slope of the annual NDVI, only where the trend is significant:

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

Green: increase; red: decrease; white, the pixels without a significant trend. With values, as GeoTIFF:

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

For one month or one 16-day block: step='monthly', part=8 (August) or step='16-day', part=14 (days 209-224). output chooses the band of the tile map: 'slope', 'tau', 'z', 'p_value', 'intercept' or 'n'.

9. Compare two periods

The NDVI of each month in 2020-2024 against 2005-2009:

c = mo.compare_series(sn, 'NDVI', 'Modis', '2005-01-01', '2009-12-31',     # reference period
                      '2020-01-01', '2024-12-31',                          # the other period
                      step='monthly', csv='comparison_ndvi.csv', chart='comparison_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'] has the values of each year of the two periods. The map of the difference (the median of 2020-2024 minus that of 2005-2009):

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

Burn severity

The August 2023 fire in the Teide National Park, with Sentinel-2 (see Concepts):

s = mo.burn_severity({'wdpa_id': 389013}, 'Sentinel-2', '2023-05-15', '2023-08-14',   # before
                     '2023-09-15', '2023-11-30', path='severity.tif', csv='severity.csv')   # after
s['burned_ha']                                  # {'El Teide': 956.4}: hectares with 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')   # the RBR, continuous
r['stats']                                      # [{'region', 'mean', 'median', 'p10', 'p90', 'min', 'max', 'ha'}]

With index='RBR' it gives the RBR instead of the dNBR and its classes (see Concepts).

10. Read the GeoTIFF files

With GDAL (it comes in the conda environment) or rasterio:

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

11. Export to Google Drive

For large maps at full resolution or very long series. Earth Engine computes them with no time limit and leaves them in a folder of your Google Drive:

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

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

The tasks can also be seen at https://code.earthengine.google.com/tasks.

12. Several variables for a model: recipes

download_variables downloads several variables and metrics on aligned grids (same CRS, same edges, NoData −9999) over a study area and, with a coarser pixel, over the ecoregions it touches. It is what eHabitat+ uses to prepare its input variables. The recipe is a YAML (or a dictionary):

crs: auto                                   # UTM of the centre of the area, or e.g. EPSG:25830
study: {resolution: 10, buffer: 500}        # study area: rectangle + 500 m, direct download
ecoregion: {resolution: 100, drive: true}   # optional: the ecoregions it touches, via 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]}                       # one file per metric
  - {index: LST, sensor: Modis, temporal: annual, metric: max, interannual: median}
  - {name: slope, image: COPERNICUS/DEM/GLO30_2024_1, band: DEM, terrain: slope}  # any Earth Engine image

First, a trial run (dry_run) to see the size of the grids and the number of images:

p = mo.download_variables('recipe.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_slope.tif']

And then the download:

mo.download_variables('recipe.yaml', 'study_area.gpkg', 'variables', log=print)

• The variables of the study area end up in variables/ and those of the ecoregions in variables/ecoregion/, with their polygons in ecoregions.gpkg.
• monitoreo.json keeps how each file was made: when you repeat it, only what has changed is downloaded (overwrite=True downloads everything).
• The grids with drive: true are exported to Google Drive (with no size limit) and imported from Google Drive for desktop when they finish. If it is interrupted, when you repeat it the exports in progress are awaited instead of launching them again.

13. Other Earth Engine data

All the functions accept, instead of a variable of the catalog, an Earth Engine collection or image:

mo.time_series(sn, collection='IDAHO_EPSCOR/TERRACLIMATE', band='pr', start='2000-01-01', end='2023-12-31',
               step='annual', metric='sum')                                   # annual precipitation
mo.region_stats(sn, image='COPERNICUS/DEM/GLO30_2024_1', band='DEM', terrain='slope')   # slope of the terrain

Reference

Function	What it does	Returns
init(project, shared)	Connects to Earth Engine: with your project or, without it (or with shared=True), with the shared account of MonitorEO, which does not export to Google Drive	the project
eecu_usage(limit, refresh, lang)	Your EECUs spent this month (research EECUs with the shared account; with your account, the quota of your project)	used_percent, left_percent, used_hours, limit_hours, warning, text, bar
about(lang)	Description, acknowledgements, logos and links	dict
list_variables(sensor, text)	Variables of the catalog	list
describe_variable(variable, sensor)	Everything about a variable	dict
find_protected_areas(text)	Protected areas of WDPA by name	list
time_series(region, variable, sensor, start, end, step, metric, stat, ...)	Time series	rows, csv, chart, scale, unit, title
map_layer(region, variable, sensor, start, end, metric, temporal, ...)	Map tiles	url, min, max, palette, bounds, thumbnail
thumbnail(...)	Preview of the map	PNG (bytes)
composite_image(...)	The Earth Engine image of the map	(ee.Image, ee.FeatureCollection)
region_stats(...)	Statistics of the map by polygon	rows, scale, unit
download_image(..., path)	Map as GeoTIFF	file, grid, resolution, ...
trend_table(..., step, alpha)	Trends of the region	rows, yearly, csv, chart
trend_map(..., step, part, output, alpha)	Tiles of the trend map	like map_layer
download_trend(..., step, part, path)	Trends by pixel as GeoTIFF	file, bands, ...
compare_series(..., other_start, other_end, step)	Comparison of periods	profile, rows, csv, chart
compare_map(...) / download_comparison(..., path)	Map of the difference	like map_layer / download_image
burn_severity(region, sensor, start, end, other_start, other_end, path, csv, index)	Burn severity (dNBR, USGS classes; or, with index='RBR', the RBR)	classes, burned_ha (or stats of the RBR), file
export_image(...) / time_series_to_drive(...)	Exports to Google Drive	operation, file, folder
task_status(operation)	State of an export	state, progress, error
download_variables(recipe, area, folder, ...)	Recipe of variables as GeoTIFF	summary

Each function has its full documentation: help(mo.time_series).