Merge branch 'main' into git_attributes_hot_fix

docs/changelog.rst

@@ -31,6 +31,7 @@ Fixed
 - Allow for string format in zoom_level path, e.g. `{zoom_level:02d}` (#851)
 - Fixed incorrect renaming of single variable raster datasets (#883)
 - Provide better error message for 0D geometry arrays in GeoDataset (#885)
+- Fixed index error when the number of peaks varies between stations in get_hydrographs method (#933)
 
 Deprecated
 ----------

hydromt/stats/design_events.py

@@ -111,7 +111,7 @@ def hydrograph_1d(
     n = ts.size
     d0 = int(np.floor(wdw_size / 2))
     d1 = wdw_size - d0
-    for i in range(n0):
+    for i in range(min(n0, idxs.size)):
         idx = idxs[seq[i]]
         idx0 = idx - d0
         idx1 = idx + d1

tests/test_stats_design_events.py

@@ -7,11 +7,13 @@
 def test_get_peak_hydrographs(ts_extremes):
     ts_extremes = ts_extremes.isel(time=slice(365 * 10))  # smaller sample
     da_peaks = extremes.get_peaks(
-        ts_extremes, ev_type="BM", period="182.625D"
+        ts_extremes,
+        ev_type="BM",
+        period="182.625D",  # this returns 20 peaks
     )  # default: ev_type='BM', period='year'
     da = design_events.get_peak_hydrographs(
         ts_extremes, da_peaks, wdw_size=7, n_peaks=20, normalize=False
-    )
+    ).load()
     assert da.time.shape[0] == 7  # wdw_size=7
     assert (da.argmax("time") == 3).all()  # peak at time=3
 
@@ -24,6 +26,13 @@ def test_get_peak_hydrographs(ts_extremes):
     # Testing if normalize values are set to 1
     da = design_events.get_peak_hydrographs(
         ts_extremes, da_peaks, wdw_size=7, normalize=True
-    )
+    ).load()
     damax_1 = da.sel(stations=1, time=0)
     assert (damax_1 == 1).all()
+
+    # test when number of peaks varies between stations
+    # set last 4 peaks of station 2 (index 1) to nan
+    da_peaks[dict(stations=1, time=slice(365 * 8, -1))] = np.nan
+    da = design_events.get_peak_hydrographs(ts_extremes, da_peaks, wdw_size=7).load()
+    assert da.isel(time=3, stations=1).count("peak") == 16
+    assert da.isel(time=3, stations=0).count("peak") == 20