Hydroperiod Analysis with ndvi2gif

This notebook demonstrates the HydroperiodAnalyzer class for computing hydroperiod (flood days per hydrological cycle) entirely in Google Earth Engine — no data downloads required.

What is hydroperiod?

Hydroperiod represents the number of days a pixel remains flooded during a hydrological cycle (Sep 1 – Aug 31 by default). It is a key indicator for wetland ecology, species habitat assessment, and hydrological monitoring.

Methodology

The midpoint weighting method distributes 365 days among available scenes: each scene is assigned a temporal territory that spans from the midpoint with its predecessor to the midpoint with its successor. Weighted flood days are accumulated per pixel:

Scene dates:    Sep-01    Oct-15    Dec-20    Mar-10    Jun-05
Midpoints:           Oct-02    Nov-17    Feb-07    Apr-23
Territories:   [0, 31]  [31, 77]  [77, 160] [160, 235] [235, 365]

Products

Band	Type	Description
hydroperiod	int16	Accumulated weighted flood days
valid_days	int16	Accumulated valid (cloud-free) observation days
normalized	int16	Hydroperiod normalised to 365-day year
first_flood_doy	int16	First flood day of the cycle (0 = Sep 1st)
last_flood_doy	int16	Last flood day of the cycle

Study Area: Doñana National Park, Spain

One of Europe’s most important wetlands and a UNESCO World Heritage Site, Doñana has a highly dynamic flooding regime driven by seasonal rainfall. It is an ideal benchmark for hydroperiod analysis.

1. Setup

import warnings
warnings.filterwarnings('ignore', category=UserWarning, module='geemap.conversion')

import ee
import geemap
from ndvi2gif import NdviSeasonality, HydroperiodAnalyzer

print('✓ Modules imported')

# Authenticate only the first time
# ee.Authenticate()
ee.Initialize()
print('✓ Earth Engine initialised')

✓ Modules imported
✓ Earth Engine initialised

2. Region of Interest

We use Doñana’s marismas (marshes) as the study area. You can replace this with any ROI supported by NdviSeasonality: shapefile, GeoJSON, drawn geometry, eLTER DEIMS ID, or Sentinel-2 tile code.

# Doñana marshes bounding box
roi = ee.Geometry.Rectangle([-6.55, 36.85, -6.10, 37.20])

# --- Alternative ROIs ---
# roi = 'path/to/donana.shp'
# roi = 'deimsid:https://deims.org/...'
# roi = Map.draw_last_feature  # draw on the map below

Map = geemap.Map()
Map.centerObject(roi, zoom=10)
Map.addLayer(roi, {'color': 'blue'}, 'Study area')
Map

3. Single Hydrological Cycle

We start with one cycle (Sep 2022 – Aug 2023) using MNDWI (Modified Normalised Difference Water Index), which is generally the most robust water index for optical imagery.

3.1 Initialise NdviSeasonality and HydroperiodAnalyzer

# NdviSeasonality provides the satellite collection, ROI, and cloud filtering.
# start_year / end_year define which cycles will be available.
ns = NdviSeasonality(
    roi=roi,
    sat='S2',
    start_year=2022,
    end_year=2023,
    index='mndwi',
    cloud_filter=True,
    max_cloud_cover=10,   # scene-level: excludes scenes with >20% clouds
    scl_mask=False       # pixel-level: SCL cloud/shadow mask (set False to use legacy QA60)
)

# HydroperiodAnalyzer adds pixel-level SCL cloud/shadow masking on top (S2 only)
ha = HydroperiodAnalyzer(ns)
print(ha)

There we go again...
Applying cloud filter to Sentinel-2: max 10% cloud cover
Applying pixel-level cloud mask using QA60 band
Using MODIS Terra + Aqua LST (maximum coverage)
Using all Sentinel-1 orbits (ascending + descending).
Applying S1 ARD preprocessing:
  - Speckle filter: REFINED_LEE
  - Terrain correction: True
  - Terrain model: VOLUME
Sentinel-2 collection configured with cloud filtering
HydroperiodAnalyzer(sat='S2', hyd_year=2022/2023, hyd_start=09-01)

3.2 Inspect water masks (optional)

Before running the full computation you can examine the binary water masks and their temporal weights.

The scl_mask=True parameter (default) applies a pixel-level cloud/shadow mask using the Sentinel-2 SCL band, which is more accurate than the QA60 band used by NdviSeasonality.

masks = ha.get_water_masks(
    index='mndwi',
    threshold=0.0,
    hyd_year=2022,
    # scl_mask is handled by NdviSeasonality (scl_mask=True by default for S2)
)

n_scenes = masks.size().getInfo()
print(f'Scenes in 2022-2023 cycle: {n_scenes}')

# Verify that weights sum to 365
total_weight = sum(masks.aggregate_array('weight').getInfo())
print(f'Sum of weights: {total_weight:.1f} days  (should be ~365)')

# Inspect the first 5 scenes
props = masks.limit(5).aggregate_array('weight').getInfo()
starts = masks.limit(5).aggregate_array('start_doy').getInfo()
ends = masks.limit(5).aggregate_array('end_doy').getInfo()

print('\nFirst 5 scenes (weight, start_doy, end_doy):')
for w, s, e in zip(props, starts, ends):
    print(f'  weight={w:.1f} days  [{s:.0f} → {e:.0f}]')

Scenes in 2022-2023 cycle: 37
Sum of weights: 365.0 days  (should be ~365)

First 5 scenes (weight, start_doy, end_doy):
  weight=4.5 days  [0 → 4]
  weight=10.0 days  [4 → 14]
  weight=10.0 days  [14 → 24]
  weight=7.5 days  [24 → 32]
  weight=7.5 days  [32 → 40]

# Visualise a single water mask on the map
first_mask = ee.Image(masks.first())

water_vis = {'min': 0, 'max': 1, 'palette': ['tan', '0000ff']}

Map1 = geemap.Map()
Map1.centerObject(roi, zoom=10)
Map1.addLayer(first_mask, water_vis, 'Water mask (first scene)')
Map1.addLayer(roi, {'color': 'red'}, 'ROI')
Map1

3.3 Compute hydroperiod

result = ha.compute_hydroperiod(
    index='mndwi',
    threshold=0.0,        # pixels with MNDWI > 0 → water
    hyd_year=2022,        # Sep 2022 – Aug 2023
    normalize=True,
    compute_first_last=True,
)

print('Bands:', result.bandNames().getInfo())
print('Properties:', result.propertyNames().getInfo())

Bands: ['hydroperiod', 'valid_days', 'normalized', 'first_flood_doy', 'last_flood_doy']
Properties: ['hyd_year_end', 'system:footprint', 'index', 'threshold', 'hyd_year_start', 'system:bands', 'system:band_names']

3.4 Visualise results

Map2 = geemap.Map()
Map2.centerObject(roi, zoom=10)

# Normalised hydroperiod (0–365 days): white = never flooded, dark blue = always flooded
hyd_vis = {
    'bands': ['normalized'],
    'min': 0, 'max': 365,
    'palette': ['ffffff', 'cce5f5', '91bfdb', '4393c3', '2166ac', '053061'],
}
Map2.addLayer(result, hyd_vis, 'Hydroperiod normalised (days)')

# Valid days (data coverage)
valid_vis = {'bands': ['valid_days'], 'min': 0, 'max': 365, 'palette': ['white', 'green']}
Map2.addLayer(result, valid_vis, 'Valid days', shown=False)

# First flood day (day of hydrological year, 0 = Sep 1st)
first_vis = {'bands': ['first_flood_doy'], 'min': 0, 'max': 180, 'palette': ['blue', 'cyan', 'yellow']}
Map2.addLayer(result, first_vis, 'First flood DOY', shown=False)

# Last flood day
last_vis = {'bands': ['last_flood_doy'], 'min': 0, 'max': 365, 'palette': ['yellow', 'orange', 'red']}
Map2.addLayer(result, last_vis, 'Last flood DOY', shown=False)

Map2.addLayer(roi, {'color': 'red'}, 'ROI')
Map2.add_colorbar(hyd_vis, label='Flood days', orientation='horizontal')
Map2

4. Trying Different Water Indices and Thresholds

The choice of index and threshold affects how conservative the water mapping is. Available indices: 'ndwi', 'mndwi', 'awei', 'aweinsh', 'wi2015'.

Map3 = geemap.Map()
Map3.centerObject(roi, zoom=10)

hyd_vis = {
    'min': 0, 'max': 365,
    'palette': ['ffffff', 'cce5f5', '91bfdb', '4393c3', '2166ac', '053061'],
}

# MNDWI — default, works well for open water
r_mndwi = ha.compute_hydroperiod(index='mndwi', threshold=0.0, hyd_year=2022)
Map3.addLayer(r_mndwi.select('normalized'), hyd_vis, 'MNDWI thr=0.0')

# NDWI — classic McFeeters index
r_ndwi = ha.compute_hydroperiod(index='ndwi', threshold=0.0, hyd_year=2022)
Map3.addLayer(r_ndwi.select('normalized'), hyd_vis, 'NDWI thr=0.0', shown=False)

# AWEI — includes shadow correction, useful in areas with mixed land cover
r_awei = ha.compute_hydroperiod(index='awei', threshold=0.0, hyd_year=2022)
Map3.addLayer(r_awei.select('normalized'), hyd_vis, 'AWEI thr=0.0', shown=False)

# Stricter threshold with MNDWI to map only permanent/deep water
r_strict = ha.compute_hydroperiod(index='mndwi', threshold=0.3, hyd_year=2022)
Map3.addLayer(r_strict.select('normalized'), hyd_vis, 'MNDWI thr=0.3 (strict)', shown=False)

Map3.addLayer(roi, {'color': 'red'}, 'ROI')
Map3

5. Multiple Hydrological Cycles

compute_all_cycles() iterates automatically over all years in the NdviSeasonality date range and returns a dict {year: ee.Image}.

# Expand the date range to cover several hydrological cycles
ns_multi = NdviSeasonality(
    roi=roi,
    sat='S2',
    start_year=2019,
    end_year=2023,
    index='mndwi',
    cloud_filter=True,
)

ha_multi = HydroperiodAnalyzer(ns_multi)

# Compute all cycles: 2019-20, 2020-21, 2021-22, 2022-23
cycles = ha_multi.compute_all_cycles(index='mndwi', threshold=0.0)

print('Computed cycles:')
for yr, img in cycles.items():
    print(f'  {yr}/{yr+1}')

# Visualise all cycles — toggle layers to compare
Map4 = geemap.Map()
Map4.centerObject(roi, zoom=10)

hyd_vis = {
    'min': 0, 'max': 365,
    'palette': ['ffffff', 'cce5f5', '91bfdb', '4393c3', '2166ac', '053061'],
}

for yr, img in cycles.items():
    Map4.addLayer(img.select('normalized'), hyd_vis, f'Hydroperiod {yr}/{yr+1}')

Map4.addLayer(roi, {'color': 'red'}, 'ROI')
Map4

5.1 Compare two cycles — inter-annual variability

Subtract two cycles to see which areas were wetter or drier in a given year.

# Δ: 2022-23 minus 2019-20
diff = cycles[2022].select('normalized').subtract(cycles[2019].select('normalized'))

diff_vis = {
    'min': -150, 'max': 150,
    'palette': ['d73027', 'fc8d59', 'fee090', 'ffffff', 'e0f3f8', '91bfdb', '4575b4'],
}

Map5 = geemap.Map()
Map5.centerObject(roi, zoom=10)
Map5.addLayer(diff, diff_vis, 'Δ Hydroperiod 2022-23 vs 2019-20 (days)')
Map5.addLayer(roi, {'color': 'red'}, 'ROI')
Map5.add_colorbar(diff_vis, label='Δ flood days  (blue = wetter in 2022-23)', orientation='horizontal')
Map5

5.2 Mean hydroperiod and anomalies

compute_anomalies() calculates the mean hydroperiod across a reference period and the anomaly of each individual cycle relative to that mean:

• Positive anomaly → wetter than average year

• Negative anomaly → drier than average year

The reference parameter controls what “average” means:

reference=	Baseline
'period' (default)	Mean of the cycles in the analysis window (start_year → end_year)
'historical'	Mean of the full satellite archive (S2 from 2017, Landsat from 1984, …)

# Compute mean hydroperiod + anomaly per cycle
# Pass the already-computed cycles dict to avoid recomputing
anomaly_result = ha_multi.compute_anomalies(cycles=cycles)

mean_hyd = anomaly_result['mean']       # ee.Image — mean_hydroperiod band
anomalies = anomaly_result['anomalies'] # {yr: ee.Image} — anomaly band per cycle

print('Mean hydroperiod band:', mean_hyd.bandNames().getInfo())
print('Anomaly cycles:', list(anomalies.keys()))

Map5b = geemap.Map()
Map5b.centerObject(roi, zoom=10)

# Mean hydroperiod
hyd_vis = {
    'min': 0, 'max': 365,
    'palette': ['ffffff', 'cce5f5', '91bfdb', '4393c3', '2166ac', '053061'],
}
Map5b.addLayer(mean_hyd, hyd_vis, 'Mean hydroperiod (all cycles)')

# Anomalies — diverging palette: red = drier, blue = wetter
anom_vis = {
    'min': -150, 'max': 150,
    'palette': ['d73027', 'fc8d59', 'fee090', 'ffffff', 'e0f3f8', '91bfdb', '4575b4'],
}
for yr, anom in anomalies.items():
    Map5b.addLayer(anom, anom_vis, f'Anomaly {yr}/{yr+1}', shown=False)

Map5b.addLayer(roi, {'color': 'red'}, 'ROI')
Map5b.add_colorbar(anom_vis, label='Anomaly (days)  blue = wetter than average', orientation='horizontal')
Map5b

5.2b Climatological baseline — reference='historical'

Using the full S2 archive (2017 – present) as the baseline instead of just the 2019–2023 window. GEE builds the extra computation graph lazily, so the call is instantaneous; actual computation happens on render or export.

# Climatological baseline: compare each 2019–2023 cycle against the
# full S2 archive mean (2017 – last complete hydrological year).
# Pass the already-computed `cycles` dict so those images are reused.
result_hist = ha_multi.compute_anomalies(
    cycles=cycles,
    reference='historical',
)

mean_hist = result_hist['mean'] # Use this variable to export the mean historic hydroperiod
anomalies_hist = result_hist['anomalies'] # Use this variable to export the mean historic anomalies

# Confirm baseline metadata
print('Reference :', mean_hist.get('reference').getInfo())
print('Baseline  :', mean_hist.get('ref_start_year').getInfo(),
      '→', mean_hist.get('ref_end_year').getInfo())
print('N cycles  :', mean_hist.get('ref_n_cycles').getInfo())

Map5c = geemap.Map()
Map5c.centerObject(roi, zoom=10)

hyd_vis = {
    'min': 0, 'max': 365,
    'palette': ['ffffff', 'cce5f5', '91bfdb', '4393c3', '2166ac', '053061'],
}
anom_vis = {
    'min': -150, 'max': 150,
    'palette': ['d73027', 'fc8d59', 'fee090', 'ffffff', 'e0f3f8', '91bfdb', '4575b4'],
}

# Mean computed from the full S2 archive
Map5c.addLayer(mean_hist, hyd_vis, 'Mean hydroperiod (S2 archive, historical)')

# Anomalies of each 2019-2023 cycle vs the historical baseline
for yr, anom in anomalies_hist.items():
    Map5c.addLayer(anom, anom_vis, f'Anomaly {yr}/{yr+1} (vs historical)', shown=False)

Map5c.addLayer(roi, {'color': 'red'}, 'ROI')
Map5c.add_colorbar(anom_vis, label='Anomaly (days)  blue = wetter than S2 archive mean', orientation='horizontal')
Map5c

6. Temporal Representativity Index (IRT)

The IRT measures how evenly distributed the cloud-free observations are across the hydrological year:

• IRT = 1: observations perfectly spread across all months

• IRT = 0: all observations concentrated in a single period

A low IRT warns that hydroperiod estimates may be unreliable.

# Global IRT (one number per cycle) — computed client-side from scene dates
for yr in range(2019, 2023):
    irt = ha_multi.compute_irt_global(hyd_year=yr)
    print(f'  IRT {yr}/{yr+1}: {irt:.3f}')

# Per-pixel IRT — highlights areas with poor temporal coverage
ha_multi.get_water_masks(index='mndwi', threshold=0.0, hyd_year=2022)
irt_img = ha_multi.compute_irt_image(hyd_year=2022)

irt_vis = {
    'bands': ['irt'],
    'min': 0, 'max': 1,
    'palette': ['d7191c', 'fdae61', 'ffffbf', 'abd9e9', '2c7bb6'],
}

Map6 = geemap.Map()
Map6.centerObject(roi, zoom=10)
Map6.addLayer(irt_img, irt_vis, 'Per-pixel IRT 2022-23')
Map6.addLayer(roi, {'color': 'red'}, 'ROI')
Map6.add_colorbar(irt_vis, label='IRT (1 = perfect coverage)', orientation='horizontal')
Map6

7. Using Landsat for Longer Time Series

Landsat provides data since 1984, enabling multi-decadal hydroperiod analysis. Useful for detecting long-term trends in wetland extent.

ns_ls = NdviSeasonality(
    roi=roi,
    sat='Landsat',
    start_year=2000,
    end_year=2005,
    index='mndwi',
    cloud_filter=True,
)

ha_ls = HydroperiodAnalyzer(ns_ls)

# scl_mask is ignored for Landsat (SCL is Sentinel-2 specific)
ls_cycles = ha_ls.compute_all_cycles(index='mndwi', threshold=0.0)

Map7 = geemap.Map()
Map7.centerObject(roi, zoom=10)

hyd_vis = {
    'min': 0, 'max': 365,
    'palette': ['ffffff', 'cce5f5', '91bfdb', '4393c3', '2166ac', '053061'],
}
for yr, img in ls_cycles.items():
    Map7.addLayer(img.select('normalized'), hyd_vis, f'Landsat hydroperiod {yr}/{yr+1}')

Map7.addLayer(roi, {'color': 'red'}, 'ROI')
Map7

8. Export Results

8.1 Export to Google Drive

# Export the 2022-23 result to Google Drive
task = ha.export_to_drive(
    folder='donana_hydroperiod',
    description='donana_hydroperiod_2022_2023',
    scale=10,          # 10m — Sentinel-2 native resolution
    crs='EPSG:25829',  # UTM zone 29N — appropriate for SW Spain
)
# Monitor at: https://code.earthengine.google.com/tasks

# Export all cycles
for yr, img in cycles.items():
    ha_multi.export_to_drive(
        image=img,
        folder='donana_hydroperiod',
        description=f'donana_hydroperiod_{yr}_{yr+1}',
        scale=10,
        crs='EPSG:25829',
    )

# Export mean hydroperiod and anomalies
ha_multi.export_to_drive(
    image=mean_hyd,
    folder='donana_hydroperiod',
    description='donana_hydroperiod_mean_2019_2023',
    scale=10,
    crs='EPSG:25829',
)
for yr, anom in anomalies.items():
    ha_multi.export_to_drive(
        image=anom,
        folder='donana_hydroperiod',
        description=f'donana_hydroperiod_anomaly_{yr}_{yr+1}',
        scale=10,
        crs='EPSG:25829',
    )

8.2 Export to Earth Engine Asset

# Replace with your EE asset path
asset_id = 'projects/your-project/assets/donana_hydroperiod_2022_2023'

task = ha.export_to_asset(
    asset_id=asset_id,
    scale=10,
    crs='EPSG:25829',
)

9. Custom Hydrological Year Start

For regions with a different rainy season (e.g. tropical areas starting in April), pass a custom hydrological_year_start to HydroperiodAnalyzer.

# Example: cycle starting April 1st (e.g. for tropical/monsoon regions)
ha_tropical = HydroperiodAnalyzer(
    ns,
    hydrological_year_start=(4, 1),  # April 1st
)

result_tropical = ha_tropical.compute_hydroperiod(
    index='mndwi',
    threshold=0.0,
    hyd_year=2022,  # Apr 2022 – Mar 2023
)

print('Properties:', result_tropical.propertyNames().getInfo())

Summary

Task	Method
Binary water masks with temporal weights	ha.get_water_masks(index, threshold, hyd_year)
Single hydrological cycle	ha.compute_hydroperiod(hyd_year=YYYY)
All cycles in date range	ha.compute_all_cycles()
Mean + anomalies (analysis window)	ha.compute_anomalies(cycles, reference='period')
Mean + anomalies (satellite archive)	ha.compute_anomalies(cycles, reference='historical')
Global IRT (one number)	ha.compute_irt_global()
Per-pixel IRT (ee.Image)	ha.compute_irt_image()
Export to Drive	ha.export_to_drive(folder=...)
Export to Asset	ha.export_to_asset(asset_id=...)
Pixel-level cloud/shadow mask (S2)	Pass scl_mask=True to NdviSeasonality (default)

Supported water indices: ndwi, mndwi, awei, aweinsh, wi2015

Supported satellites: S2 (10 m, with SCL masking), Landsat (30 m), MODIS (500 m)

Historical archive starts: S2 → 2017, Landsat → 1984, MODIS → 2000, Sentinel-1 → 2014

References

Related Python library

• phydroperiod (García Díaz & Bustamante Díaz, 2026) — standalone Python library for hydroperiod computation from water masks. Available on PyPI: pip install phydroperiod. Repository: https://pypi.org/project/phydroperiod/. Use phydroperiod when you already have water masks as local raster files or NumPy arrays; use the HydroperiodAnalyzer in ndvi2gif when you want to run the full pipeline (satellite collection → water masks → hydroperiod) server-side on Google Earth Engine.

Methodology and context

• Bustamante, J., Aragonés, D., Afán, I., Luque, C.J., Pérez-Vázquez, A., Castellanos, E.M., Díaz-Delgado, R. (2016). Predictive models of floristic composition based on flooding frequency in Mediterranean wetlands. Remote Sensing, 8(9), 776. https://doi.org/10.3390/rs8090776

Water indices

• McFeeters, S.K. (1996). The use of the Normalized Difference Water Index (NDWI) in the delineation of open water features. International Journal of Remote Sensing, 17(7), 1425–1432.

• Xu, H. (2006). Modification of normalised difference water index (NDWI) to enhance open water features in remotely sensed imagery. International Journal of Remote Sensing, 27(14), 3025–3033.

• Feyisa, G.L., Meilby, H., Fensholt, R., Proud, S.R. (2014). Automated Water Extraction Index: A new technique for surface water mapping using Landsat imagery. Remote Sensing of Environment, 140, 23–35.

• Fisher, A., Flood, N., Danaher, T. (2016). Comparing Landsat water index methods for automated water classification in eastern Australia. Remote Sensing of Environment, 175, 167–182.