Skip to content

Meteorological data

The meteorological data module handles both station-based and gridded (raster) meteorological forcing for MOBIDIC simulations. It also includes tools for generating synthetic design storm hyetographs from IDF (Intensity-Duration-Frequency) parameters.

Overview

MOBIDICpy supports three modes of meteorological forcing, each suited to different workflows:

Station-based forcing (MeteoData)

Use when you have time-series from weather stations that need spatial interpolation.

  • Station time series: Time-series from weather stations
  • Spatial interpolation: Performed during simulation (IDW or Nearest Neighbor)
  • Format support: MATLAB .mat files and CF-compliant NetCDF
  • Optional grid export: Can save interpolated grids as NetCDF for reuse

Raster-based forcing (MeteoRaster)

Use when you have pre-interpolated gridded data or want to reuse previously interpolated results.

  • Gridded data: Reads pre-interpolated NetCDF files with dimensions (time, y, x)
  • Performance: Faster simulation by skipping real-time interpolation
  • Grid validation: Automatically validates that rasters match the model grid
  • Best for: Production runs, scenario analysis, or when using external gridded datasets (e.g., reanalysis products)

Design storm hyetographs (HyetographGenerator)

Use for design storm simulations with synthetic precipitation from IDF curves.

  • IDF-based generation: Constructs precipitation from spatially distributed IDF parameters (a, n, k)
  • Chicago method: Currently implements the Chicago decreasing hyetograph (after-peak curve)
  • Areal reduction: Applies areal reduction factor (ARF) coefficient for spatial averaging
  • Auto-resampling: IDF parameter rasters are automatically resampled to match the model grid
  • Best for: Design flood simulations, flood mapping, infrastructure sizing

IDF formula

The total precipitation depth for a given duration is computed using the Depth-Duration-Frequency (DDF) relationship:

\[DDF(t) = k_a \cdot k \cdot a \cdot t^n\]

where:

  • \(DDF(t)\) is the cumulative precipitation depth (mm) for duration \(t\)
  • \(k_a\) is the areal reduction factor (ARF) coefficient
  • \(k\) is the return period factor (spatially distributed raster)
  • \(a\) is the IDF scale parameter (spatially distributed raster)
  • \(n\) is the IDF exponent parameter (spatially distributed raster)
  • \(t\) is duration in hours

Workflow recommendation:

  1. Initial run: Use station-based forcing with output_forcing_data.meteo_data = True
  2. Subsequent runs: Use the exported raster forcing for faster performance
  3. Different interpolation: Re-run with station forcing if you need to change interpolation method

API Reference

Input data classes

MeteoData

Container for station-based meteorological data with spatial interpolation capability.

Container for meteorological station data.

This class holds meteorological observations from multiple stations, organized by variable type (precipitation, temperature, wind, etc.).

Attributes:

Name Type Description
stations

Dictionary mapping variable names to lists of station dictionaries
variables

List of available meteorological variables
start_date

First timestamp in the dataset
end_date

Last timestamp in the dataset
Source code in mobidic/preprocessing/meteo_preprocessing.py 
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
class MeteoData:
    """Container for meteorological station data.

    This class holds meteorological observations from multiple stations,
    organized by variable type (precipitation, temperature, wind, etc.).

    Attributes:
        stations: Dictionary mapping variable names to lists of station dictionaries
        variables: List of available meteorological variables
        start_date: First timestamp in the dataset
        end_date: Last timestamp in the dataset
    """

    def __init__(self, stations: dict[str, list[dict[str, Any]]]):
        """Initialize MeteoData container.

        Args:
            stations: Dictionary with variable names as keys ('precipitation',
                'temperature_min', 'temperature_max', 'humidity', 'wind_speed',
                'radiation') and lists of station dictionaries as values.
                Each station dict contains: 'code', 'x', 'y', 'elevation',
                'time', 'data', 'name' (optional).
        """
        self.stations = stations
        self.variables = list(stations.keys())

        # Extract date range from all stations
        all_times = []
        for var_stations in stations.values():
            for station in var_stations:
                if len(station["time"]) > 0:
                    all_times.extend(station["time"])

        if all_times:
            self.start_date = pd.Timestamp(min(all_times))
            self.end_date = pd.Timestamp(max(all_times))
        else:
            self.start_date = None
            self.end_date = None

    def __repr__(self):
        n_stations = {var: len(stations) for var, stations in self.stations.items()}
        return (
            f"MeteoData(variables={self.variables}, "
            f"n_stations={n_stations}, "
            f"period={self.start_date} to {self.end_date})"
        )

    def to_netcdf(
        self,
        output_path: str | Path,
        compression_level: int = 4,
        add_metadata: dict[str, Any] | None = None,
    ) -> None:
        """Save meteorological data to CF-compliant NetCDF file.

        Args:
            output_path: Path to output NetCDF file
            compression_level: Compression level for NetCDF (0-9, default 4)
            add_metadata: Additional global attributes to add to NetCDF file

        Examples:
            >>> meteo = MeteoData.from_mat("meteodata.mat")
            >>> meteo.to_netcdf("meteodata.nc", add_metadata={"basin": "Arno"})
        """
        output_path = Path(output_path)
        logger.info(f"Saving meteorological data to NetCDF: {output_path}")

        # Create xarray Dataset with station data organized by variable
        # Group related variables that should share station dimensions
        # temperature_min and temperature_max share the same stations
        var_groups = {
            "temperature": ["temperature_min", "temperature_max"],
            "precipitation": ["precipitation"],
            "humidity": ["humidity"],
            "wind_speed": ["wind_speed"],
            "radiation": ["radiation"],
        }

        datasets = {}

        for group_name, var_list in var_groups.items():
            # Check which variables in this group are available
            available_vars = [v for v in var_list if v in self.stations and len(self.stations[v]) > 0]

            if len(available_vars) == 0:
                continue

            # Use the first available variable to define the station dimension
            # For temperature, this ensures min and max share the same stations
            primary_var = available_vars[0]
            var_stations = self.stations[primary_var]

            # Extract station metadata from primary variable
            n_stations = len(var_stations)
            station_codes = np.array([s["code"] for s in var_stations])
            station_x = np.array([s["x"] for s in var_stations])
            station_y = np.array([s["y"] for s in var_stations])
            station_elevation = np.array([s["elevation"] for s in var_stations])
            station_names = np.array([s.get("name", "") for s in var_stations], dtype=str)

            # Collect all unique timestamps across all variables in this group
            all_times = []
            for var_name in available_vars:
                for station in self.stations[var_name]:
                    all_times.extend(station["time"])
            unique_times = pd.DatetimeIndex(sorted(set(all_times)))
            n_times = len(unique_times)

            # Station dimension name based on group (e.g., "station_temperature")
            station_dim = f"station_{group_name}"

            # Create data variables for all variables in this group
            data_vars = {}
            for var_name in available_vars:
                var_stations = self.stations[var_name]

                # Create data array (time x station)
                data_array = np.full((n_times, n_stations), np.nan, dtype=np.float32)

                # Fill data array
                for i, station in enumerate(var_stations):
                    if len(station["time"]) > 0:
                        # Find indices of this station's times in the unified time array
                        time_indices = unique_times.get_indexer(pd.DatetimeIndex(station["time"]))
                        valid_mask = time_indices >= 0
                        data_array[time_indices[valid_mask], i] = station["data"][valid_mask]

                data_vars[var_name] = (
                    ["time", station_dim],
                    data_array,
                    {
                        "long_name": _get_variable_longname(var_name),
                        "units": _get_variable_units(var_name),
                        "missing_value": np.nan,
                    },
                )

            # Create xarray Dataset for this group
            ds = xr.Dataset(
                data_vars,
                coords={
                    "time": (["time"], unique_times, {"long_name": "time", "axis": "T"}),
                    station_dim: ([station_dim], np.arange(n_stations), {"long_name": "station index"}),
                    f"station_code_{group_name}": (
                        [station_dim],
                        station_codes,
                        {"long_name": f"{group_name} station code"},
                    ),
                    f"x_{group_name}": (
                        [station_dim],
                        station_x,
                        {"long_name": f"{group_name} station x coordinate (easting)", "units": "m"},
                    ),
                    f"y_{group_name}": (
                        [station_dim],
                        station_y,
                        {"long_name": f"{group_name} station y coordinate (northing)", "units": "m"},
                    ),
                    f"elevation_{group_name}": (
                        [station_dim],
                        station_elevation,
                        {"long_name": f"{group_name} station elevation", "units": "m"},
                    ),
                    f"station_name_{group_name}": (
                        [station_dim],
                        station_names,
                        {"long_name": f"{group_name} station name"},
                    ),
                },
            )

            datasets[group_name] = ds

        # Merge all variable datasets
        if len(datasets) == 0:
            raise ValueError("No valid meteorological data to save")

        merged_ds = xr.merge(datasets.values(), compat="no_conflicts")

        # Add global attributes
        merged_ds.attrs["title"] = "MOBIDIC meteorological forcing data"
        merged_ds.attrs["source"] = "Station observations"
        merged_ds.attrs["Conventions"] = "CF-1.12"
        merged_ds.attrs["history"] = f"Created by MOBIDICpy v{__version__}"

        if add_metadata:
            merged_ds.attrs.update(add_metadata)

        # Set encoding for compression
        encoding = {}
        for var in merged_ds.data_vars:
            encoding[var] = {
                "zlib": True,
                "complevel": compression_level,
                "dtype": "float32",
                "_FillValue": np.nan,
            }

        # Save to NetCDF
        merged_ds.to_netcdf(
            output_path,
            encoding=encoding,
            format="NETCDF4",
        )

        logger.success(f"Saved meteorological data: {len(datasets)} variables, {output_path}")

    @classmethod
    def from_mat(cls, mat_path: str | Path) -> "MeteoData":
        """Load meteorological data from MATLAB .mat file (legacy MOBIDIC format).

        Args:
            mat_path: Path to MATLAB .mat file containing meteodata

        Returns:
            MeteoData object with station data

        Examples:
            >>> meteo = MeteoData.from_mat("examples/Arno/meteodata/meteodata.mat")
            >>> print(meteo)
            >>> meteo.to_netcdf("meteodata.nc")
        """
        reader = MATMeteoReader(mat_path)
        return reader.read()

    @classmethod
    def from_netcdf(cls, nc_path: str | Path) -> "MeteoData":
        """Load meteorological data from NetCDF file.

        Args:
            nc_path: Path to NetCDF file containing meteodata

        Returns:
            MeteoData object with station data

        Examples:
            >>> meteo = MeteoData.from_netcdf("meteodata.nc")
            >>> print(meteo)
        """
        reader = NetCDFMeteoReader(nc_path)
        return reader.read()

    @classmethod
    def from_csv(cls, csv_path: str | Path, config: dict[str, Any]) -> "MeteoData":
        """Load meteorological data from CSV file(s).

        Args:
            csv_path: Path to CSV file or directory containing CSV files
            config: Configuration dictionary specifying CSV structure

        Returns:
            MeteoData object with station data

        Note:
            This method is planned for future implementation.
        """
        raise NotImplementedError("CSV reader not yet implemented")

__init__(stations)

Initialize MeteoData container.

Parameters:

Name Type Description Default
stations dict[str, list[dict[str, Any]]]

Dictionary with variable names as keys (‘precipitation’, ‘temperature_min’, ‘temperature_max’, ‘humidity’, ‘wind_speed’, ‘radiation’) and lists of station dictionaries as values. Each station dict contains: ‘code’, ‘x’, ‘y’, ‘elevation’, ‘time’, ‘data’, ‘name’ (optional).

 required
Source code in mobidic/preprocessing/meteo_preprocessing.py 
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
def __init__(self, stations: dict[str, list[dict[str, Any]]]):
    """Initialize MeteoData container.

    Args:
        stations: Dictionary with variable names as keys ('precipitation',
            'temperature_min', 'temperature_max', 'humidity', 'wind_speed',
            'radiation') and lists of station dictionaries as values.
            Each station dict contains: 'code', 'x', 'y', 'elevation',
            'time', 'data', 'name' (optional).
    """
    self.stations = stations
    self.variables = list(stations.keys())

    # Extract date range from all stations
    all_times = []
    for var_stations in stations.values():
        for station in var_stations:
            if len(station["time"]) > 0:
                all_times.extend(station["time"])

    if all_times:
        self.start_date = pd.Timestamp(min(all_times))
        self.end_date = pd.Timestamp(max(all_times))
    else:
        self.start_date = None
        self.end_date = None

from_csv(csv_path, config) classmethod

Load meteorological data from CSV file(s).

Parameters:

Name Type Description Default
csv_path str | Path

Path to CSV file or directory containing CSV files

 required
config dict[str, Any]

Configuration dictionary specifying CSV structure

 required

Returns:

Type Description
MeteoData

MeteoData object with station data
Note

This method is planned for future implementation.

Source code in mobidic/preprocessing/meteo_preprocessing.py 
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
@classmethod
def from_csv(cls, csv_path: str | Path, config: dict[str, Any]) -> "MeteoData":
    """Load meteorological data from CSV file(s).

    Args:
        csv_path: Path to CSV file or directory containing CSV files
        config: Configuration dictionary specifying CSV structure

    Returns:
        MeteoData object with station data

    Note:
        This method is planned for future implementation.
    """
    raise NotImplementedError("CSV reader not yet implemented")

from_mat(mat_path) classmethod

Load meteorological data from MATLAB .mat file (legacy MOBIDIC format).

Parameters:

Name Type Description Default
mat_path str | Path

Path to MATLAB .mat file containing meteodata

 required

Returns:

Type Description
MeteoData

MeteoData object with station data

Examples:

>>> meteo = MeteoData.from_mat("examples/Arno/meteodata/meteodata.mat")
>>> print(meteo)
>>> meteo.to_netcdf("meteodata.nc")
Source code in mobidic/preprocessing/meteo_preprocessing.py 
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
@classmethod
def from_mat(cls, mat_path: str | Path) -> "MeteoData":
    """Load meteorological data from MATLAB .mat file (legacy MOBIDIC format).

    Args:
        mat_path: Path to MATLAB .mat file containing meteodata

    Returns:
        MeteoData object with station data

    Examples:
        >>> meteo = MeteoData.from_mat("examples/Arno/meteodata/meteodata.mat")
        >>> print(meteo)
        >>> meteo.to_netcdf("meteodata.nc")
    """
    reader = MATMeteoReader(mat_path)
    return reader.read()

from_netcdf(nc_path) classmethod

Load meteorological data from NetCDF file.

Parameters:

Name Type Description Default
nc_path str | Path

Path to NetCDF file containing meteodata

 required

Returns:

Type Description
MeteoData

MeteoData object with station data

Examples:

>>> meteo = MeteoData.from_netcdf("meteodata.nc")
>>> print(meteo)
Source code in mobidic/preprocessing/meteo_preprocessing.py 
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
@classmethod
def from_netcdf(cls, nc_path: str | Path) -> "MeteoData":
    """Load meteorological data from NetCDF file.

    Args:
        nc_path: Path to NetCDF file containing meteodata

    Returns:
        MeteoData object with station data

    Examples:
        >>> meteo = MeteoData.from_netcdf("meteodata.nc")
        >>> print(meteo)
    """
    reader = NetCDFMeteoReader(nc_path)
    return reader.read()

to_netcdf(output_path, compression_level=4, add_metadata=None)

Save meteorological data to CF-compliant NetCDF file.

Parameters:

Name Type Description Default
output_path str | Path

Path to output NetCDF file

 required
compression_level int

Compression level for NetCDF (0-9, default 4)

 4
add_metadata dict[str, Any] | None

Additional global attributes to add to NetCDF file

 None

Examples:

>>> meteo = MeteoData.from_mat("meteodata.mat")
>>> meteo.to_netcdf("meteodata.nc", add_metadata={"basin": "Arno"})
Source code in mobidic/preprocessing/meteo_preprocessing.py 
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
def to_netcdf(
    self,
    output_path: str | Path,
    compression_level: int = 4,
    add_metadata: dict[str, Any] | None = None,
) -> None:
    """Save meteorological data to CF-compliant NetCDF file.

    Args:
        output_path: Path to output NetCDF file
        compression_level: Compression level for NetCDF (0-9, default 4)
        add_metadata: Additional global attributes to add to NetCDF file

    Examples:
        >>> meteo = MeteoData.from_mat("meteodata.mat")
        >>> meteo.to_netcdf("meteodata.nc", add_metadata={"basin": "Arno"})
    """
    output_path = Path(output_path)
    logger.info(f"Saving meteorological data to NetCDF: {output_path}")

    # Create xarray Dataset with station data organized by variable
    # Group related variables that should share station dimensions
    # temperature_min and temperature_max share the same stations
    var_groups = {
        "temperature": ["temperature_min", "temperature_max"],
        "precipitation": ["precipitation"],
        "humidity": ["humidity"],
        "wind_speed": ["wind_speed"],
        "radiation": ["radiation"],
    }

    datasets = {}

    for group_name, var_list in var_groups.items():
        # Check which variables in this group are available
        available_vars = [v for v in var_list if v in self.stations and len(self.stations[v]) > 0]

        if len(available_vars) == 0:
            continue

        # Use the first available variable to define the station dimension
        # For temperature, this ensures min and max share the same stations
        primary_var = available_vars[0]
        var_stations = self.stations[primary_var]

        # Extract station metadata from primary variable
        n_stations = len(var_stations)
        station_codes = np.array([s["code"] for s in var_stations])
        station_x = np.array([s["x"] for s in var_stations])
        station_y = np.array([s["y"] for s in var_stations])
        station_elevation = np.array([s["elevation"] for s in var_stations])
        station_names = np.array([s.get("name", "") for s in var_stations], dtype=str)

        # Collect all unique timestamps across all variables in this group
        all_times = []
        for var_name in available_vars:
            for station in self.stations[var_name]:
                all_times.extend(station["time"])
        unique_times = pd.DatetimeIndex(sorted(set(all_times)))
        n_times = len(unique_times)

        # Station dimension name based on group (e.g., "station_temperature")
        station_dim = f"station_{group_name}"

        # Create data variables for all variables in this group
        data_vars = {}
        for var_name in available_vars:
            var_stations = self.stations[var_name]

            # Create data array (time x station)
            data_array = np.full((n_times, n_stations), np.nan, dtype=np.float32)

            # Fill data array
            for i, station in enumerate(var_stations):
                if len(station["time"]) > 0:
                    # Find indices of this station's times in the unified time array
                    time_indices = unique_times.get_indexer(pd.DatetimeIndex(station["time"]))
                    valid_mask = time_indices >= 0
                    data_array[time_indices[valid_mask], i] = station["data"][valid_mask]

            data_vars[var_name] = (
                ["time", station_dim],
                data_array,
                {
                    "long_name": _get_variable_longname(var_name),
                    "units": _get_variable_units(var_name),
                    "missing_value": np.nan,
                },
            )

        # Create xarray Dataset for this group
        ds = xr.Dataset(
            data_vars,
            coords={
                "time": (["time"], unique_times, {"long_name": "time", "axis": "T"}),
                station_dim: ([station_dim], np.arange(n_stations), {"long_name": "station index"}),
                f"station_code_{group_name}": (
                    [station_dim],
                    station_codes,
                    {"long_name": f"{group_name} station code"},
                ),
                f"x_{group_name}": (
                    [station_dim],
                    station_x,
                    {"long_name": f"{group_name} station x coordinate (easting)", "units": "m"},
                ),
                f"y_{group_name}": (
                    [station_dim],
                    station_y,
                    {"long_name": f"{group_name} station y coordinate (northing)", "units": "m"},
                ),
                f"elevation_{group_name}": (
                    [station_dim],
                    station_elevation,
                    {"long_name": f"{group_name} station elevation", "units": "m"},
                ),
                f"station_name_{group_name}": (
                    [station_dim],
                    station_names,
                    {"long_name": f"{group_name} station name"},
                ),
            },
        )

        datasets[group_name] = ds

    # Merge all variable datasets
    if len(datasets) == 0:
        raise ValueError("No valid meteorological data to save")

    merged_ds = xr.merge(datasets.values(), compat="no_conflicts")

    # Add global attributes
    merged_ds.attrs["title"] = "MOBIDIC meteorological forcing data"
    merged_ds.attrs["source"] = "Station observations"
    merged_ds.attrs["Conventions"] = "CF-1.12"
    merged_ds.attrs["history"] = f"Created by MOBIDICpy v{__version__}"

    if add_metadata:
        merged_ds.attrs.update(add_metadata)

    # Set encoding for compression
    encoding = {}
    for var in merged_ds.data_vars:
        encoding[var] = {
            "zlib": True,
            "complevel": compression_level,
            "dtype": "float32",
            "_FillValue": np.nan,
        }

    # Save to NetCDF
    merged_ds.to_netcdf(
        output_path,
        encoding=encoding,
        format="NETCDF4",
    )

    logger.success(f"Saved meteorological data: {len(datasets)} variables, {output_path}")

MeteoRaster

Container for pre-interpolated raster meteorological forcing.

Container for gridded meteorological data from rasters.

Loads CF-1.12 compliant NetCDF files with dimensions (time, y, x). By default, preloads all data into memory for optimal performance.

Attributes:

Name Type Description
nc_path

Path to NetCDF file
ds

xarray Dataset (preloaded by default, or lazy-loaded if preload=False)
variables

List of available meteorological variables
time_range

Tuple of (start_time, end_time)
start_date

First timestamp in the dataset (property)
end_date

Last timestamp in the dataset (property)
grid_metadata

Dictionary with shape, resolution, crs, origin
Performance

Default preload=True: Loads entire dataset into memory at initialization. This provides fast access during simulation (comparable to station interpolation). Use preload=False for very large datasets to use lazy loading (slower but less memory).

Source code in mobidic/preprocessing/meteo_raster.py 
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
class MeteoRaster:
    """Container for gridded meteorological data from rasters.

    Loads CF-1.12 compliant NetCDF files with dimensions (time, y, x).
    By default, preloads all data into memory for optimal performance.

    Attributes:
        nc_path: Path to NetCDF file
        ds: xarray Dataset (preloaded by default, or lazy-loaded if preload=False)
        variables: List of available meteorological variables
        time_range: Tuple of (start_time, end_time)
        start_date: First timestamp in the dataset (property)
        end_date: Last timestamp in the dataset (property)
        grid_metadata: Dictionary with shape, resolution, crs, origin

    Performance:
        Default preload=True: Loads entire dataset into memory at initialization.
        This provides fast access during simulation (comparable to station interpolation).
        Use preload=False for very large datasets to use lazy loading (slower but less memory).
    """

    def __init__(self, nc_path: str | Path, preload: bool = True):
        """Initialize from NetCDF file.

        Args:
            nc_path: Path to NetCDF file containing raster meteorological data
            preload: If True, load all data into memory immediately (default: True).
                If False, use lazy loading (slower but uses less memory).

        Raises:
            FileNotFoundError: If the specified file does not exist
            ValueError: If the NetCDF file has invalid structure
        """
        self.nc_path = Path(nc_path)
        if not self.nc_path.exists():
            raise FileNotFoundError(f"Meteorological raster file not found: {self.nc_path}")

        logger.info(f"Loading meteorological rasters from: {self.nc_path}")

        # Load NetCDF file with xarray
        if preload:
            # Load all data into memory for fast access
            logger.debug("Preloading all meteorological data into memory...")
            self.ds = xr.open_dataset(self.nc_path).load()
            logger.debug("Preloading complete")
        else:
            # Use lazy loading (slower but uses less memory)
            logger.debug("Using lazy loading (data read on-demand from disk)")
            self.ds = xr.open_dataset(self.nc_path)

        # Validate structure
        self._validate_structure()

        # Extract metadata
        self._extract_metadata()

        # Initialize cache for current timestep (less useful when preloaded, but kept for compatibility)
        self._cache: dict[tuple[str, str], np.ndarray] = {}

        logger.success(
            f"Loaded meteorological rasters: {len(self.variables)} variables, "
            f"{len(self.ds.time)} timesteps, grid={self.grid_metadata['shape']}"
        )

    def _validate_structure(self) -> None:
        """Validate that NetCDF has required dimensions and structure.

        Raises:
            ValueError: If structure is invalid
        """
        # Check required dimensions
        if "time" not in self.ds.dims:
            raise ValueError("NetCDF file must have 'time' dimension")
        if "y" not in self.ds.dims:
            raise ValueError("NetCDF file must have 'y' dimension")
        if "x" not in self.ds.dims:
            raise ValueError("NetCDF file must have 'x' dimension")

        # Check that we have at least one data variable (besides crs and z)
        data_vars = [v for v in self.ds.data_vars if v not in ("crs", "z")]
        if len(data_vars) == 0:
            raise ValueError("NetCDF file must have at least one meteorological variable")

        logger.debug(f"NetCDF structure validated: {self.ds.dims}")

    def _extract_metadata(self) -> None:
        """Extract grid metadata from NetCDF file."""
        # Extract variables (exclude 'crs' and optional 'z' coordinate variables)
        self.variables = [v for v in self.ds.data_vars if v not in ("crs", "z")]

        # Extract time range
        times = self.ds.time.values
        self.time_range = (times[0], times[-1])

        # Extract grid metadata
        nrows = len(self.ds.y)
        ncols = len(self.ds.x)

        # Calculate resolution (assumes regular grid)
        if nrows > 1:
            y_res = abs(float(self.ds.y[1] - self.ds.y[0]))
        else:
            y_res = None

        if ncols > 1:
            x_res = abs(float(self.ds.x[1] - self.ds.x[0]))
        else:
            x_res = None

        # Check that resolution is consistent
        if y_res is not None and x_res is not None and not np.isclose(y_res, x_res, rtol=1e-6):
            logger.warning(f"Grid has different x and y resolution: x={x_res}, y={y_res}")

        resolution = x_res if x_res is not None else y_res

        # Extract origin (lower-left corner)
        # y coordinate is typically in descending order (north to south)
        y_values = self.ds.y.values
        x_values = self.ds.x.values
        yllcorner = float(y_values.min())
        xllcorner = float(x_values.min())

        # Extract CRS if available
        crs = None
        if "crs" in self.ds:
            crs_var = self.ds["crs"]
            if hasattr(crs_var, "spatial_ref"):
                crs = crs_var.attrs.get("spatial_ref")
            elif hasattr(crs_var, "crs_wkt"):
                crs = crs_var.attrs.get("crs_wkt")

        # Convert CRS to EPSG code for cleaner logging output
        epsg = get_epsg_code(crs)
        crs_display = f"EPSG:{epsg}" if epsg else crs

        self.grid_metadata = {
            "shape": (nrows, ncols),
            "resolution": resolution,
            "xllcorner": xllcorner,
            "yllcorner": yllcorner,
            "crs": crs_display,
        }

        logger.debug(f"Grid metadata: {self.grid_metadata}")

    @property
    def start_date(self):
        """First timestamp in the dataset.

        Returns the start date as pandas Timestamp for consistency with MeteoData API.
        """
        return pd.Timestamp(self.time_range[0])

    @property
    def end_date(self):
        """Last timestamp in the dataset.

        Returns the end date as pandas Timestamp for consistency with MeteoData API.
        """
        return pd.Timestamp(self.time_range[1])

    @classmethod
    def from_netcdf(cls, nc_path: str | Path, preload: bool = True) -> "MeteoRaster":
        """Load meteorological raster data from NetCDF file.

        This is an alias for __init__ to match the MeteoData API.

        Args:
            nc_path: Path to NetCDF file containing raster meteorological data
            preload: If True, load all data into memory immediately (default: True).
                If False, use lazy loading (slower but uses less memory).

        Returns:
            MeteoRaster object

        Examples:
            >>> # Fast loading (preload into memory)
            >>> meteo = MeteoRaster.from_netcdf("Arno_meteoraster.nc")
            >>> print(meteo)
            >>>
            >>> # Lazy loading (less memory, slower)
            >>> meteo = MeteoRaster.from_netcdf("Arno_meteoraster.nc", preload=False)
        """
        return cls(nc_path, preload=preload)

    def get_timestep(self, time: datetime, variable: str) -> np.ndarray:
        """Extract 2D grid for a variable at a specific time.

        Uses nearest neighbor time sampling. Results are cached to avoid
        repeated disk reads for the same timestep.

        Args:
            time: Datetime to extract (uses nearest neighbor)
            variable: Variable name (e.g., 'precipitation', 'pet')

        Returns:
            2D numpy array (nrows, ncols) in file units (mm/h for precip/pet)

        Raises:
            KeyError: If variable not found in dataset
        """
        # Check if variable exists
        if variable not in self.variables:
            raise KeyError(f"Variable '{variable}' not found in raster data. Available variables: {self.variables}")

        # Check cache first
        cache_key = (variable, time.isoformat())
        if cache_key in self._cache:
            return self._cache[cache_key]

        # Extract data using xarray's nearest neighbor time selection
        # Note: When data is preloaded, this operation is fast (in-memory)
        # When using lazy loading, this triggers disk I/O
        data = self.ds[variable].sel(time=time, method="nearest").values

        # Cache the result
        self._cache[cache_key] = data

        return data

    def validate_grid_alignment(self, model_metadata: dict) -> None:
        """Validate that raster grid matches model grid.

        Checks shape, resolution, and origin alignment. Raises detailed
        error if grids don't match.

        Args:
            model_metadata: Dictionary with model grid metadata containing:
                - shape: (nrows, ncols)
                - resolution: grid resolution in meters
                - xllcorner: x coordinate of lower-left corner
                - yllcorner: y coordinate of lower-left corner
                - crs: coordinate reference system (optional)

        Raises:
            ValueError: If grids are not aligned with detailed mismatch info
        """
        logger.info("Validating grid alignment between raster and model")

        errors = []

        # Check shape
        if self.grid_metadata["shape"] != model_metadata["shape"]:
            errors.append(f"Shape mismatch: raster={self.grid_metadata['shape']}, model={model_metadata['shape']}")

        # Check resolution
        raster_res = self.grid_metadata["resolution"]
        model_res = model_metadata["resolution"]
        if raster_res is not None and model_res is not None:
            # Handle model_res as tuple (x_res, y_res) or scalar
            if isinstance(model_res, (tuple, list, np.ndarray)):
                model_res_scalar = model_res[0]  # Use x resolution
            else:
                model_res_scalar = model_res

            if not np.isclose(raster_res, model_res_scalar, rtol=1e-6):
                errors.append(f"Resolution mismatch: raster={raster_res:.6f}m, model={model_res_scalar:.6f}m")

        # Check origin (allow 1mm tolerance for floating point precision)
        if not np.isclose(self.grid_metadata["xllcorner"], model_metadata["xllcorner"], atol=1e-3):
            errors.append(
                f"X origin mismatch: raster={self.grid_metadata['xllcorner']:.6f}, "
                f"model={model_metadata['xllcorner']:.6f}"
            )

        if not np.isclose(self.grid_metadata["yllcorner"], model_metadata["yllcorner"], atol=1e-3):
            errors.append(
                f"Y origin mismatch: raster={self.grid_metadata['yllcorner']:.6f}, "
                f"model={model_metadata['yllcorner']:.6f}"
            )

        # Check CRS (warning only, not an error)
        raster_crs = self.grid_metadata["crs"]
        model_crs = model_metadata.get("crs")
        if raster_crs and model_crs:
            if not crs_equals(raster_crs, model_crs):
                logger.warning(
                    "CRS mismatch (proceeding anyway): "
                    "raster CRS != model CRS. "
                    "Ensure grids are in same coordinate system."
                )

        # Raise error if any mismatches found
        if errors:
            error_msg = "Meteorological raster grid does not match model grid:\n  - " + "\n  - ".join(errors)
            logger.error(error_msg)
            raise ValueError(error_msg)

        logger.success("Grid alignment validated: raster matches model grid")

    def close(self) -> None:
        """Close NetCDF file and clear cache."""
        logger.debug("Closing meteorological raster dataset")
        self.ds.close()
        self._cache.clear()

    def __repr__(self):
        return (
            f"MeteoRaster(variables={self.variables}, "
            f"n_times={len(self.ds.time)}, "
            f"shape={self.grid_metadata['shape']}, "
            f"period={self.start_date} to {self.end_date})"
        )

    def __enter__(self):
        """Context manager entry."""
        return self

    def __exit__(self, exc_type, exc_val, exc_tb):
        """Context manager exit."""
        self.close()
        return False

end_date property

Last timestamp in the dataset.

Returns the end date as pandas Timestamp for consistency with MeteoData API.

start_date property

First timestamp in the dataset.

Returns the start date as pandas Timestamp for consistency with MeteoData API.

__enter__()

Context manager entry.

Source code in mobidic/preprocessing/meteo_raster.py 
332
333
334
def __enter__(self):
    """Context manager entry."""
    return self

__exit__(exc_type, exc_val, exc_tb)

Context manager exit.

Source code in mobidic/preprocessing/meteo_raster.py 
336
337
338
339
def __exit__(self, exc_type, exc_val, exc_tb):
    """Context manager exit."""
    self.close()
    return False

__init__(nc_path, preload=True)

Initialize from NetCDF file.

Parameters:

Name Type Description Default
nc_path str | Path

Path to NetCDF file containing raster meteorological data

 required
preload bool

If True, load all data into memory immediately (default: True). If False, use lazy loading (slower but uses less memory).

 True

Raises:

Type Description
FileNotFoundError

If the specified file does not exist
ValueError

If the NetCDF file has invalid structure
Source code in mobidic/preprocessing/meteo_raster.py 
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
def __init__(self, nc_path: str | Path, preload: bool = True):
    """Initialize from NetCDF file.

    Args:
        nc_path: Path to NetCDF file containing raster meteorological data
        preload: If True, load all data into memory immediately (default: True).
            If False, use lazy loading (slower but uses less memory).

    Raises:
        FileNotFoundError: If the specified file does not exist
        ValueError: If the NetCDF file has invalid structure
    """
    self.nc_path = Path(nc_path)
    if not self.nc_path.exists():
        raise FileNotFoundError(f"Meteorological raster file not found: {self.nc_path}")

    logger.info(f"Loading meteorological rasters from: {self.nc_path}")

    # Load NetCDF file with xarray
    if preload:
        # Load all data into memory for fast access
        logger.debug("Preloading all meteorological data into memory...")
        self.ds = xr.open_dataset(self.nc_path).load()
        logger.debug("Preloading complete")
    else:
        # Use lazy loading (slower but uses less memory)
        logger.debug("Using lazy loading (data read on-demand from disk)")
        self.ds = xr.open_dataset(self.nc_path)

    # Validate structure
    self._validate_structure()

    # Extract metadata
    self._extract_metadata()

    # Initialize cache for current timestep (less useful when preloaded, but kept for compatibility)
    self._cache: dict[tuple[str, str], np.ndarray] = {}

    logger.success(
        f"Loaded meteorological rasters: {len(self.variables)} variables, "
        f"{len(self.ds.time)} timesteps, grid={self.grid_metadata['shape']}"
    )

close()

Close NetCDF file and clear cache.

Source code in mobidic/preprocessing/meteo_raster.py 
318
319
320
321
322
def close(self) -> None:
    """Close NetCDF file and clear cache."""
    logger.debug("Closing meteorological raster dataset")
    self.ds.close()
    self._cache.clear()

from_netcdf(nc_path, preload=True) classmethod

Load meteorological raster data from NetCDF file.

This is an alias for init to match the MeteoData API.

Parameters:

Name Type Description Default
nc_path str | Path

Path to NetCDF file containing raster meteorological data

 required
preload bool

If True, load all data into memory immediately (default: True). If False, use lazy loading (slower but uses less memory).

 True

Returns:

Type Description
MeteoRaster

MeteoRaster object

Examples:

>>> # Fast loading (preload into memory)
>>> meteo = MeteoRaster.from_netcdf("Arno_meteoraster.nc")
>>> print(meteo)
>>>
>>> # Lazy loading (less memory, slower)
>>> meteo = MeteoRaster.from_netcdf("Arno_meteoraster.nc", preload=False)
Source code in mobidic/preprocessing/meteo_raster.py 
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
@classmethod
def from_netcdf(cls, nc_path: str | Path, preload: bool = True) -> "MeteoRaster":
    """Load meteorological raster data from NetCDF file.

    This is an alias for __init__ to match the MeteoData API.

    Args:
        nc_path: Path to NetCDF file containing raster meteorological data
        preload: If True, load all data into memory immediately (default: True).
            If False, use lazy loading (slower but uses less memory).

    Returns:
        MeteoRaster object

    Examples:
        >>> # Fast loading (preload into memory)
        >>> meteo = MeteoRaster.from_netcdf("Arno_meteoraster.nc")
        >>> print(meteo)
        >>>
        >>> # Lazy loading (less memory, slower)
        >>> meteo = MeteoRaster.from_netcdf("Arno_meteoraster.nc", preload=False)
    """
    return cls(nc_path, preload=preload)

get_timestep(time, variable)

Extract 2D grid for a variable at a specific time.

Uses nearest neighbor time sampling. Results are cached to avoid repeated disk reads for the same timestep.

Parameters:

Name Type Description Default
time datetime

Datetime to extract (uses nearest neighbor)

 required
variable str

Variable name (e.g., ‘precipitation’, ‘pet’)

 required

Returns:

Type Description
ndarray

2D numpy array (nrows, ncols) in file units (mm/h for precip/pet)

Raises:

Type Description
KeyError

If variable not found in dataset
Source code in mobidic/preprocessing/meteo_raster.py 
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
def get_timestep(self, time: datetime, variable: str) -> np.ndarray:
    """Extract 2D grid for a variable at a specific time.

    Uses nearest neighbor time sampling. Results are cached to avoid
    repeated disk reads for the same timestep.

    Args:
        time: Datetime to extract (uses nearest neighbor)
        variable: Variable name (e.g., 'precipitation', 'pet')

    Returns:
        2D numpy array (nrows, ncols) in file units (mm/h for precip/pet)

    Raises:
        KeyError: If variable not found in dataset
    """
    # Check if variable exists
    if variable not in self.variables:
        raise KeyError(f"Variable '{variable}' not found in raster data. Available variables: {self.variables}")

    # Check cache first
    cache_key = (variable, time.isoformat())
    if cache_key in self._cache:
        return self._cache[cache_key]

    # Extract data using xarray's nearest neighbor time selection
    # Note: When data is preloaded, this operation is fast (in-memory)
    # When using lazy loading, this triggers disk I/O
    data = self.ds[variable].sel(time=time, method="nearest").values

    # Cache the result
    self._cache[cache_key] = data

    return data

validate_grid_alignment(model_metadata)

Validate that raster grid matches model grid.

Checks shape, resolution, and origin alignment. Raises detailed error if grids don’t match.

Parameters:

Name Type Description Default
model_metadata dict

Dictionary with model grid metadata containing: - shape: (nrows, ncols) - resolution: grid resolution in meters - xllcorner: x coordinate of lower-left corner - yllcorner: y coordinate of lower-left corner - crs: coordinate reference system (optional)

 required

Raises:

Type Description
ValueError

If grids are not aligned with detailed mismatch info
Source code in mobidic/preprocessing/meteo_raster.py 
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
def validate_grid_alignment(self, model_metadata: dict) -> None:
    """Validate that raster grid matches model grid.

    Checks shape, resolution, and origin alignment. Raises detailed
    error if grids don't match.

    Args:
        model_metadata: Dictionary with model grid metadata containing:
            - shape: (nrows, ncols)
            - resolution: grid resolution in meters
            - xllcorner: x coordinate of lower-left corner
            - yllcorner: y coordinate of lower-left corner
            - crs: coordinate reference system (optional)

    Raises:
        ValueError: If grids are not aligned with detailed mismatch info
    """
    logger.info("Validating grid alignment between raster and model")

    errors = []

    # Check shape
    if self.grid_metadata["shape"] != model_metadata["shape"]:
        errors.append(f"Shape mismatch: raster={self.grid_metadata['shape']}, model={model_metadata['shape']}")

    # Check resolution
    raster_res = self.grid_metadata["resolution"]
    model_res = model_metadata["resolution"]
    if raster_res is not None and model_res is not None:
        # Handle model_res as tuple (x_res, y_res) or scalar
        if isinstance(model_res, (tuple, list, np.ndarray)):
            model_res_scalar = model_res[0]  # Use x resolution
        else:
            model_res_scalar = model_res

        if not np.isclose(raster_res, model_res_scalar, rtol=1e-6):
            errors.append(f"Resolution mismatch: raster={raster_res:.6f}m, model={model_res_scalar:.6f}m")

    # Check origin (allow 1mm tolerance for floating point precision)
    if not np.isclose(self.grid_metadata["xllcorner"], model_metadata["xllcorner"], atol=1e-3):
        errors.append(
            f"X origin mismatch: raster={self.grid_metadata['xllcorner']:.6f}, "
            f"model={model_metadata['xllcorner']:.6f}"
        )

    if not np.isclose(self.grid_metadata["yllcorner"], model_metadata["yllcorner"], atol=1e-3):
        errors.append(
            f"Y origin mismatch: raster={self.grid_metadata['yllcorner']:.6f}, "
            f"model={model_metadata['yllcorner']:.6f}"
        )

    # Check CRS (warning only, not an error)
    raster_crs = self.grid_metadata["crs"]
    model_crs = model_metadata.get("crs")
    if raster_crs and model_crs:
        if not crs_equals(raster_crs, model_crs):
            logger.warning(
                "CRS mismatch (proceeding anyway): "
                "raster CRS != model CRS. "
                "Ensure grids are in same coordinate system."
            )

    # Raise error if any mismatches found
    if errors:
        error_msg = "Meteorological raster grid does not match model grid:\n  - " + "\n  - ".join(errors)
        logger.error(error_msg)
        raise ValueError(error_msg)

    logger.success("Grid alignment validated: raster matches model grid")

HyetographGenerator

Generator for synthetic hyetographs from IDF parameters.

Generator for synthetic hyetographs from IDF parameters.

This class constructs synthetic rainfall time series (hyetographs) from spatially distributed IDF parameters.

Currently implements the Chicago method (decreasing after-peak curve only).

The total precipitation depth for a given duration is computed as

DDF(t) = ka * k * a * t^n

where
  • ka is the areal reduction factor (ARF) coefficient
  • k is the return period factor (spatially distributed)
  • a is the IDF scale parameter (spatially distributed)
  • n is the IDF exponent parameter (spatially distributed)
  • t is duration in hours

Attributes:

Name Type Description
idf_params

IDFParameters object with a, n, k grids
ka

Areal reduction factor coefficient

Examples:

>>> # Simplest workflow: generate from configuration (recommended)
>>> from mobidic import load_config, load_gisdata, Simulation
>>> config = load_config("basin_hyetograph.yaml")
>>> gisdata = load_gisdata(config.paths.gisdata, config.paths.network)
>>> forcing = HyetographGenerator.from_config(
...     config=config,
...     base_path="basin_dir",
...     start_time=datetime(2000, 1, 1)
... )
>>> sim = Simulation(gisdata, forcing, config)
>>> results = sim.run(forcing.start_date, forcing.end_date)
>>>
>>> # Alternative: generate forcing with manual parameters
>>> generator = HyetographGenerator.from_rasters(
...     a_raster="idf/a.tif",
...     n_raster="idf/n.tif",
...     k_raster="idf/k30.tif",
...     ka=0.8,
...     ref_raster="dem.tif"
... )
>>> forcing = generator.generate_forcing(
...     duration_hours=48,
...     start_time=datetime(2023, 11, 1),
...     output_path="design_storm.nc",
...     add_metadata={"return_period": "30 years"}
... )
>>>
>>> # Advanced workflow: manual control over generation and export
>>> times, precip = generator.generate(
...     duration_hours=48,
...     start_time=datetime(2023, 11, 1),
...     method="chicago_decreasing"
... )
>>> generator.to_netcdf(
...     "hyetograph.nc",
...     times=times,
...     precipitation=precip,
...     add_metadata={"event": "design_storm_30yr"}
... )
Source code in mobidic/preprocessing/hyetograph.py 
 485
 486
 487
 488
 489
 490
 491
 492
 493
 494
 495
 496
 497
 498
 499
 500
 501
 502
 503
 504
 505
 506
 507
 508
 509
 510
 511
 512
 513
 514
 515
 516
 517
 518
 519
 520
 521
 522
 523
 524
 525
 526
 527
 528
 529
 530
 531
 532
 533
 534
 535
 536
 537
 538
 539
 540
 541
 542
 543
 544
 545
 546
 547
 548
 549
 550
 551
 552
 553
 554
 555
 556
 557
 558
 559
 560
 561
 562
 563
 564
 565
 566
 567
 568
 569
 570
 571
 572
 573
 574
 575
 576
 577
 578
 579
 580
 581
 582
 583
 584
 585
 586
 587
 588
 589
 590
 591
 592
 593
 594
 595
 596
 597
 598
 599
 600
 601
 602
 603
 604
 605
 606
 607
 608
 609
 610
 611
 612
 613
 614
 615
 616
 617
 618
 619
 620
 621
 622
 623
 624
 625
 626
 627
 628
 629
 630
 631
 632
 633
 634
 635
 636
 637
 638
 639
 640
 641
 642
 643
 644
 645
 646
 647
 648
 649
 650
 651
 652
 653
 654
 655
 656
 657
 658
 659
 660
 661
 662
 663
 664
 665
 666
 667
 668
 669
 670
 671
 672
 673
 674
 675
 676
 677
 678
 679
 680
 681
 682
 683
 684
 685
 686
 687
 688
 689
 690
 691
 692
 693
 694
 695
 696
 697
 698
 699
 700
 701
 702
 703
 704
 705
 706
 707
 708
 709
 710
 711
 712
 713
 714
 715
 716
 717
 718
 719
 720
 721
 722
 723
 724
 725
 726
 727
 728
 729
 730
 731
 732
 733
 734
 735
 736
 737
 738
 739
 740
 741
 742
 743
 744
 745
 746
 747
 748
 749
 750
 751
 752
 753
 754
 755
 756
 757
 758
 759
 760
 761
 762
 763
 764
 765
 766
 767
 768
 769
 770
 771
 772
 773
 774
 775
 776
 777
 778
 779
 780
 781
 782
 783
 784
 785
 786
 787
 788
 789
 790
 791
 792
 793
 794
 795
 796
 797
 798
 799
 800
 801
 802
 803
 804
 805
 806
 807
 808
 809
 810
 811
 812
 813
 814
 815
 816
 817
 818
 819
 820
 821
 822
 823
 824
 825
 826
 827
 828
 829
 830
 831
 832
 833
 834
 835
 836
 837
 838
 839
 840
 841
 842
 843
 844
 845
 846
 847
 848
 849
 850
 851
 852
 853
 854
 855
 856
 857
 858
 859
 860
 861
 862
 863
 864
 865
 866
 867
 868
 869
 870
 871
 872
 873
 874
 875
 876
 877
 878
 879
 880
 881
 882
 883
 884
 885
 886
 887
 888
 889
 890
 891
 892
 893
 894
 895
 896
 897
 898
 899
 900
 901
 902
 903
 904
 905
 906
 907
 908
 909
 910
 911
 912
 913
 914
 915
 916
 917
 918
 919
 920
 921
 922
 923
 924
 925
 926
 927
 928
 929
 930
 931
 932
 933
 934
 935
 936
 937
 938
 939
 940
 941
 942
 943
 944
 945
 946
 947
 948
 949
 950
 951
 952
 953
 954
 955
 956
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
class HyetographGenerator:
    """Generator for synthetic hyetographs from IDF parameters.

    This class constructs synthetic rainfall time series (hyetographs) from
    spatially distributed IDF parameters.

    Currently implements the Chicago method (decreasing after-peak curve only).

    The total precipitation depth for a given duration is computed as:
        DDF(t) = ka * k * a * t^n

    where:
        - ka is the areal reduction factor (ARF) coefficient
        - k is the return period factor (spatially distributed)
        - a is the IDF scale parameter (spatially distributed)
        - n is the IDF exponent parameter (spatially distributed)
        - t is duration in hours

    Attributes:
        idf_params: IDFParameters object with a, n, k grids
        ka: Areal reduction factor coefficient

    Examples:
        >>> # Simplest workflow: generate from configuration (recommended)
        >>> from mobidic import load_config, load_gisdata, Simulation
        >>> config = load_config("basin_hyetograph.yaml")
        >>> gisdata = load_gisdata(config.paths.gisdata, config.paths.network)
        >>> forcing = HyetographGenerator.from_config(
        ...     config=config,
        ...     base_path="basin_dir",
        ...     start_time=datetime(2000, 1, 1)
        ... )
        >>> sim = Simulation(gisdata, forcing, config)
        >>> results = sim.run(forcing.start_date, forcing.end_date)
        >>>
        >>> # Alternative: generate forcing with manual parameters
        >>> generator = HyetographGenerator.from_rasters(
        ...     a_raster="idf/a.tif",
        ...     n_raster="idf/n.tif",
        ...     k_raster="idf/k30.tif",
        ...     ka=0.8,
        ...     ref_raster="dem.tif"
        ... )
        >>> forcing = generator.generate_forcing(
        ...     duration_hours=48,
        ...     start_time=datetime(2023, 11, 1),
        ...     output_path="design_storm.nc",
        ...     add_metadata={"return_period": "30 years"}
        ... )
        >>>
        >>> # Advanced workflow: manual control over generation and export
        >>> times, precip = generator.generate(
        ...     duration_hours=48,
        ...     start_time=datetime(2023, 11, 1),
        ...     method="chicago_decreasing"
        ... )
        >>> generator.to_netcdf(
        ...     "hyetograph.nc",
        ...     times=times,
        ...     precipitation=precip,
        ...     add_metadata={"event": "design_storm_30yr"}
        ... )
    """

    def __init__(self, idf_params: IDFParameters, ka: float = 1.0):
        """Initialize HyetographGenerator with IDF parameters.

        Args:
            idf_params: IDFParameters object containing a, n, k grids
            ka: Areal reduction factor (ARF) coefficient (default: 1.0)
        """
        self.idf_params = idf_params
        self.ka = ka

        logger.debug(f"HyetographGenerator initialized: ka={ka}, shape={idf_params.shape}")

    @classmethod
    def from_rasters(
        cls,
        a_raster: str | Path,
        n_raster: str | Path,
        k_raster: str | Path,
        ka: float = 1.0,
        ref_raster: str | Path | None = None,
    ) -> "HyetographGenerator":
        """Create HyetographGenerator by loading IDF parameters from raster files.

        If a reference raster is provided, the IDF parameters will be resampled
        to match its extent, resolution, and CRS using nearest neighbor interpolation.
        This is the typical workflow when the IDF rasters have different resolution
        than the model grid (e.g., DEM).

        Args:
            a_raster: Path to raster file containing IDF 'a' parameter
            n_raster: Path to raster file containing IDF 'n' parameter
            k_raster: Path to raster file containing IDF 'k' parameter
            ka: Areal reduction factor (ARF) coefficient (default: 1.0)
            ref_raster: Optional path to reference raster (e.g., DEM) for resampling.
                If provided, IDF parameters will be resampled to match this grid.

        Returns:
            HyetographGenerator instance

        Examples:
            >>> # Without resampling (IDF rasters already aligned)
            >>> generator = HyetographGenerator.from_rasters(
            ...     a_raster="idf/a.tif",
            ...     n_raster="idf/n.tif",
            ...     k_raster="idf/k30.tif",
            ...     ka=0.8
            ... )
            >>>
            >>> # With resampling to DEM grid
            >>> generator = HyetographGenerator.from_rasters(
            ...     a_raster="idf/a.tif",
            ...     n_raster="idf/n.tif",
            ...     k_raster="idf/k30.tif",
            ...     ka=0.8,
            ...     ref_raster="dem.tif"
            ... )
        """
        if ref_raster is not None:
            idf_params = read_idf_parameters_resampled(a_raster, n_raster, k_raster, ref_raster)
        else:
            idf_params = read_idf_parameters(a_raster, n_raster, k_raster)
        return cls(idf_params, ka=ka)

    @classmethod
    def from_config(
        cls,
        config: "MOBIDICConfig",
        base_path: str | Path,
        start_time: datetime,
        preload: bool = True,
    ) -> "MeteoRaster":
        """Create hyetograph forcing from parameters specified in configuration file.

        Convenience method that reads hyetograph parameters from a yaml
        configuration object, generates the hyetograph, saves to NetCDF, and
        returns a MeteoRaster ready for simulation.

        All parameters (duration, timestep, method, IDF rasters, output path) are
        read from the configuration file. Only the start time needs to be specified.
        The start time is a reference datetime for the hyetograph event.

        Args:
            config: MOBIDICConfig object with hyetograph configuration section
            base_path: Base path for resolving relative paths in config (typically
                the directory containing the config file)
            start_time: Start datetime for the hyetograph
            preload: If True, preload all data into memory for fast access
                (default: True, recommended for normal use)

        Returns:
            MeteoRaster object ready for use in Simulation

        Raises:
            AttributeError: If config does not have a hyetograph section
            ValueError: If required configuration parameters are missing

        Examples:
            >>> from mobidic import load_config, load_gisdata, Simulation
            >>> from mobidic.preprocessing.hyetograph import HyetographGenerator
            >>> from datetime import datetime
            >>> from pathlib import Path
            >>>
            >>> # Load configuration
            >>> config_file = Path("basin_hyetograph.yaml")
            >>> config = load_config(config_file)
            >>> gisdata = load_gisdata(config.paths.gisdata, config.paths.network)
            >>>
            >>> # Generate hyetograph
            >>> forcing = HyetographGenerator.from_config(
            ...     config=config,
            ...     base_path=config_file.parent,
            ...     start_time=datetime(2000, 1, 1)
            ... )
            >>>
            >>> # Run simulation
            >>> sim = Simulation(gisdata, forcing, config)
            >>> results = sim.run(forcing.start_date, forcing.end_date)

        Notes:
            - Automatically resamples IDF parameters to DEM grid
            - Uses all hyetograph parameters from config (duration, timestep, method, ka)
            - Output path read from config.paths.hyetograph
            - Creates metadata from config basin and hyetograph sections
            - Returns MeteoRaster ready for simulation with proper date range
        """
        # Import here to avoid circular dependency

        base_path = Path(base_path)

        # Validate config has hyetograph section
        if not hasattr(config, "hyetograph"):
            raise AttributeError("Configuration does not have a 'hyetograph' section")

        # Validate config has paths.hyetograph
        if not hasattr(config.paths, "hyetograph") or config.paths.hyetograph is None:
            raise ValueError(
                "Configuration must specify 'paths.hyetograph' for the output NetCDF file. "
                "Add 'hyetograph: path/to/output.nc' to the 'paths' section in the config file."
            )

        hyeto_config = config.hyetograph

        # Resolve paths relative to base_path
        a_raster_path = base_path / hyeto_config.a_raster
        n_raster_path = base_path / hyeto_config.n_raster
        k_raster_path = base_path / hyeto_config.k_raster
        ref_raster_path = base_path / config.raster_files.dtm
        output_path = base_path / config.paths.hyetograph

        # Ensure output directory exists
        output_path.parent.mkdir(parents=True, exist_ok=True)

        logger.info("Generating hyetograph forcing from configuration")
        logger.debug(f"  Base path: {base_path}")
        logger.debug(f"  Duration: {hyeto_config.duration_hours} hours")
        logger.debug(f"  Timestep: {hyeto_config.timestep_hours} hour(s)")
        logger.debug(f"  Method: {hyeto_config.hyetograph_type}")
        logger.debug(f"  Output: {output_path}")

        # Create generator with IDF parameters resampled to DEM grid
        generator = cls.from_rasters(
            a_raster=a_raster_path,
            n_raster=n_raster_path,
            k_raster=k_raster_path,
            ka=hyeto_config.ka,
            ref_raster=ref_raster_path,
        )

        # Prepare metadata from config
        add_metadata = {
            "hyetograph_method": hyeto_config.hyetograph_type,
            "duration_hours": hyeto_config.duration_hours,
            "timestep_hours": hyeto_config.timestep_hours,
            "areal_reduction_factor": hyeto_config.ka,
            "k_raster": str(hyeto_config.k_raster),
        }

        # Add basin metadata if available
        if hasattr(config, "basin"):
            if hasattr(config.basin, "id") and config.basin.id:
                add_metadata["basin"] = config.basin.id
            if hasattr(config.basin, "paramset_id") and config.basin.paramset_id:
                add_metadata["scenario"] = config.basin.paramset_id

        # Generate forcing using generate_forcing method
        forcing = generator.generate_forcing(
            duration_hours=hyeto_config.duration_hours,
            start_time=start_time,
            output_path=output_path,
            method=hyeto_config.hyetograph_type,
            timestep_hours=hyeto_config.timestep_hours,
            add_metadata=add_metadata,
            preload=preload,
        )

        return forcing

    def generate(
        self,
        duration_hours: int,
        start_time: datetime,
        method: Literal["chicago_decreasing"] = "chicago_decreasing",
        timestep_hours: int = 1,
    ) -> tuple[list[datetime], np.ndarray]:
        """Generate hyetograph precipitation time series.

        Args:
            duration_hours: Total duration of the hyetograph in hours
            start_time: Start datetime for the hyetograph
            method: Hyetograph construction method. Currently only
                "chicago_decreasing" is implemented.
            timestep_hours: Time step in hours (default: 1)

        Returns:
            Tuple of (times, precipitation) where:
                - times: List of datetime objects for each timestep
                - precipitation: 3D array (time, y, x) of precipitation (mm/h)

        Raises:
            ValueError: If method is not supported

        Notes:
            - The Chicago decreasing method generates only the falling part
              of the Chicago hyetograph (after the peak)
            - Precipitation values are in mm/h (intensity)
            - NaN values in IDF parameters propagate to output
            - If NaN values are > 90% in any parameter, a ValueError is raised during
              initialization of the HyetographGenerator
        """
        if method != "chicago_decreasing":
            raise ValueError(f"Unsupported hyetograph method: {method}. Only 'chicago_decreasing' is implemented.")

        logger.info(
            f"Generating {method} hyetograph: duration={duration_hours}h, "
            f"timestep={timestep_hours}h, start={start_time}"
        )

        return self._chicago_decreasing(duration_hours, start_time, timestep_hours)

    def _chicago_decreasing(
        self,
        duration_hours: int,
        start_time: datetime,
        timestep_hours: int,
    ) -> tuple[list[datetime], np.ndarray]:
        """Generate Chicago decreasing hyetograph.

        The Chicago method constructs a hyetograph where precipitation intensity
        decreases monotonically from the peak. This implementation generates only
        the descending part (after peak).

        Args:
            duration_hours: Total duration in hours
            start_time: Start datetime
            timestep_hours: Time step in hours

        Returns:
            Tuple of (times, precipitation) where precipitation is in mm/h
        """
        n_steps = duration_hours // timestep_hours
        nrows, ncols = self.idf_params.shape

        # Extract parameters
        a = self.idf_params.a
        n = self.idf_params.n
        k = self.idf_params.k
        ka = self.ka

        # Initialize arrays
        # DDF: Depth-Duration-Frequency, cumulative depth (mm)
        # P: Incremental precipitation (mm per timestep)
        ddf = np.zeros((n_steps, nrows, ncols))
        precip = np.zeros((n_steps, nrows, ncols))

        # Generate times
        times = [start_time + timedelta(hours=i * timestep_hours) for i in range(n_steps)]

        # Chicago hyetograph calculation
        # Note: t is 1-indexed (t=1 for first hour)
        for i in range(n_steps):
            t = (i + 1) * timestep_hours  # Duration in hours (1, 2, 3, ...)

            # DDF calculation: h = ka * k * a * t^n (mm)
            ddf[i, :, :] = ka * k * (a * np.power(t, n))

            if i > 0:
                # Rainfall increment (Chicago decreasing)
                precip[i, :, :] = ddf[i, :, :] - ddf[i - 1, :, :]
            else:
                # First timestep: P = DDF
                precip[i, :, :] = ddf[i, :, :]

        # Convert from mm/timestep to mm/h (intensity)
        precip_intensity = precip / timestep_hours

        logger.success(
            f"Hyetograph generated: {n_steps} timesteps, "
            f"total depth range=[{np.nanmin(ddf[-1]):.1f}, {np.nanmax(ddf[-1]):.1f}] mm, "
            f"peak intensity range=[{np.nanmin(precip_intensity[0]):.2f}, {np.nanmax(precip_intensity[0]):.2f}] mm/h"
        )

        return times, precip_intensity

    def generate_forcing(
        self,
        duration_hours: int,
        start_time: datetime,
        output_path: str | Path,
        method: Literal["chicago_decreasing"] = "chicago_decreasing",
        timestep_hours: int = 1,
        add_metadata: dict | None = None,
        preload: bool = True,
    ) -> "MeteoRaster":
        """Generate hyetograph and return as MeteoRaster ready for simulation.

        Convenience method that combines generate(), to_netcdf(), and
        MeteoRaster.from_netcdf() into a single call. This simplifies the
        workflow for design storm simulations.

        Args:
            duration_hours: Total duration of the hyetograph in hours
            start_time: Start datetime for the hyetograph
            output_path: Path for output NetCDF file
            method: Hyetograph construction method (default: "chicago_decreasing")
            timestep_hours: Time step in hours (default: 1)
            add_metadata: Optional dictionary of additional global attributes
            preload: If True, preload all data into memory for fast access
                (default: True, recommended for normal use)

        Returns:
            MeteoRaster object ready for use in Simulation

        Examples:
            >>> # Create generator from IDF rasters
            >>> generator = HyetographGenerator.from_rasters(
            ...     a_raster="idf/a.tif",
            ...     n_raster="idf/n.tif",
            ...     k_raster="idf/k30.tif",
            ...     ka=0.8,
            ...     ref_raster="dem.tif"
            ... )
            >>>
            >>> # Generate forcing and get MeteoRaster in one call
            >>> forcing = generator.generate_forcing(
            ...     duration_hours=48,
            ...     start_time=datetime(2023, 11, 1),
            ...     output_path="design_hyetograph.nc",
            ...     add_metadata={"return_period": "30 years"}
            ... )
            >>>
            >>> # Use directly in simulation
            >>> sim = Simulation(gisdata, forcing, config)
            >>> results = sim.run(forcing.start_date, forcing.end_date)

        Notes:
            - This method is equivalent to calling generate(), to_netcdf(),
              and MeteoRaster.from_netcdf() sequentially
            - The NetCDF file is still created at output_path for later use
        """
        # Import here to avoid circular dependency
        from mobidic.preprocessing.meteo_raster import MeteoRaster

        # Generate hyetograph
        times, precipitation = self.generate(
            duration_hours=duration_hours,
            start_time=start_time,
            method=method,
            timestep_hours=timestep_hours,
        )

        # Save to NetCDF
        self.to_netcdf(
            output_path=output_path,
            times=times,
            precipitation=precipitation,
            add_metadata=add_metadata,
        )

        # Load as MeteoRaster
        forcing = MeteoRaster.from_netcdf(output_path, preload=preload)

        logger.success(f"Forcing data ready for simulation: {forcing.start_date} to {forcing.end_date}")

        return forcing

    def to_netcdf(
        self,
        output_path: str | Path,
        times: list[datetime],
        precipitation: np.ndarray,
        add_metadata: dict | None = None,
    ) -> None:
        """Export hyetograph to CF-compliant NetCDF file.

        Creates a NetCDF file compatible with MeteoRaster.from_netcdf() for
        use as meteorological forcing in MOBIDIC simulations.

        Args:
            output_path: Path for output NetCDF file
            times: List of datetime objects for each timestep
            precipitation: 3D array (time, y, x) of precipitation [mm/h]
            add_metadata: Optional dictionary of additional global attributes

        Notes:
            - Output follows CF-1.12 conventions
            - Includes CRS information as grid mapping variable
            - Precipitation units are mm/h (compatible with MeteoRaster)

        Examples:
            >>> times, precip = generator.generate(48, datetime(2023, 11, 1))
            >>> generator.to_netcdf(
            ...     "design_storm.nc",
            ...     times=times,
            ...     precipitation=precip,
            ...     add_metadata={"return_period": "30 years"}
            ... )
        """
        output_path = Path(output_path)
        logger.info(f"Writing hyetograph to NetCDF: {output_path}")

        # Build coordinate arrays
        nrows, ncols = self.idf_params.shape
        cellsize = self.idf_params.cellsize

        # X coordinates (cell centers, west to east)
        x = np.arange(ncols) * cellsize + self.idf_params.xllcorner

        # Y coordinates (cell centers, south to north)
        # Note: Data is stored with y increasing (south to north) after flipud in grid_to_matrix
        y = np.arange(nrows) * cellsize + self.idf_params.yllcorner

        # Create xarray Dataset
        ds = xr.Dataset(
            data_vars={
                "precipitation": (
                    ["time", "y", "x"],
                    precipitation,
                    {
                        "units": "mm h-1",
                        "long_name": "Precipitation rate",
                        "grid_mapping": "crs",
                    },
                ),
            },
            coords={
                "time": times,
                "y": y,
                "x": x,
            },
            attrs={
                "Conventions": "CF-1.12",
                "title": "Synthetic hyetograph from IDF parameters",
                "source": f"MOBIDICpy version {__version__}",
                "history": f"Created {datetime.now().isoformat()}",
                "hyetograph_method": "chicago_decreasing",
                "areal_reduction_factor": self.ka,
            },
        )

        # Add coordinate attributes
        ds.x.attrs = {
            "units": "m",
            "long_name": "x coordinate",
            "standard_name": "projection_x_coordinate",
        }
        ds.y.attrs = {
            "units": "m",
            "long_name": "y coordinate",
            "standard_name": "projection_y_coordinate",
        }
        ds.time.attrs = {
            "long_name": "time",
            "standard_name": "time",
        }

        # Add CRS variable
        if self.idf_params.crs is not None:
            crs_attrs = crs_to_cf_attrs(self.idf_params.crs)
            # Create scalar CRS variable
            ds["crs"] = xr.DataArray(
                0,  # Scalar placeholder
                attrs=crs_attrs,
            )

        # Add custom metadata
        if add_metadata:
            ds.attrs.update(add_metadata)

        # Write to NetCDF with compression
        encoding = {
            "precipitation": {
                "zlib": True,
                "complevel": 4,
                "dtype": "float32",
            },
            "time": {"dtype": "float64"},
        }

        ds.to_netcdf(output_path, encoding=encoding)

        logger.success(f"Hyetograph written to: {output_path} ({precipitation.shape[0]} timesteps)")

__init__(idf_params, ka=1.0)

Initialize HyetographGenerator with IDF parameters.

Parameters:

Name Type Description Default
idf_params IDFParameters

IDFParameters object containing a, n, k grids

 required
ka float

Areal reduction factor (ARF) coefficient (default: 1.0)

 1.0
Source code in mobidic/preprocessing/hyetograph.py 
549
550
551
552
553
554
555
556
557
558
559
def __init__(self, idf_params: IDFParameters, ka: float = 1.0):
    """Initialize HyetographGenerator with IDF parameters.

    Args:
        idf_params: IDFParameters object containing a, n, k grids
        ka: Areal reduction factor (ARF) coefficient (default: 1.0)
    """
    self.idf_params = idf_params
    self.ka = ka

    logger.debug(f"HyetographGenerator initialized: ka={ka}, shape={idf_params.shape}")

from_config(config, base_path, start_time, preload=True) classmethod

Create hyetograph forcing from parameters specified in configuration file.

Convenience method that reads hyetograph parameters from a yaml configuration object, generates the hyetograph, saves to NetCDF, and returns a MeteoRaster ready for simulation.

All parameters (duration, timestep, method, IDF rasters, output path) are read from the configuration file. Only the start time needs to be specified. The start time is a reference datetime for the hyetograph event.

Parameters:

Name Type Description Default
config MOBIDICConfig

MOBIDICConfig object with hyetograph configuration section

 required
base_path str | Path

Base path for resolving relative paths in config (typically the directory containing the config file)

 required
start_time datetime

Start datetime for the hyetograph

 required
preload bool

If True, preload all data into memory for fast access (default: True, recommended for normal use)

 True

Returns:

Type Description
MeteoRaster

MeteoRaster object ready for use in Simulation

Raises:

Type Description
AttributeError

If config does not have a hyetograph section
ValueError

If required configuration parameters are missing

Examples:

>>> from mobidic import load_config, load_gisdata, Simulation
>>> from mobidic.preprocessing.hyetograph import HyetographGenerator
>>> from datetime import datetime
>>> from pathlib import Path
>>>
>>> # Load configuration
>>> config_file = Path("basin_hyetograph.yaml")
>>> config = load_config(config_file)
>>> gisdata = load_gisdata(config.paths.gisdata, config.paths.network)
>>>
>>> # Generate hyetograph
>>> forcing = HyetographGenerator.from_config(
...     config=config,
...     base_path=config_file.parent,
...     start_time=datetime(2000, 1, 1)
... )
>>>
>>> # Run simulation
>>> sim = Simulation(gisdata, forcing, config)
>>> results = sim.run(forcing.start_date, forcing.end_date)
Notes
  • Automatically resamples IDF parameters to DEM grid
  • Uses all hyetograph parameters from config (duration, timestep, method, ka)
  • Output path read from config.paths.hyetograph
  • Creates metadata from config basin and hyetograph sections
  • Returns MeteoRaster ready for simulation with proper date range
Source code in mobidic/preprocessing/hyetograph.py 
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
@classmethod
def from_config(
    cls,
    config: "MOBIDICConfig",
    base_path: str | Path,
    start_time: datetime,
    preload: bool = True,
) -> "MeteoRaster":
    """Create hyetograph forcing from parameters specified in configuration file.

    Convenience method that reads hyetograph parameters from a yaml
    configuration object, generates the hyetograph, saves to NetCDF, and
    returns a MeteoRaster ready for simulation.

    All parameters (duration, timestep, method, IDF rasters, output path) are
    read from the configuration file. Only the start time needs to be specified.
    The start time is a reference datetime for the hyetograph event.

    Args:
        config: MOBIDICConfig object with hyetograph configuration section
        base_path: Base path for resolving relative paths in config (typically
            the directory containing the config file)
        start_time: Start datetime for the hyetograph
        preload: If True, preload all data into memory for fast access
            (default: True, recommended for normal use)

    Returns:
        MeteoRaster object ready for use in Simulation

    Raises:
        AttributeError: If config does not have a hyetograph section
        ValueError: If required configuration parameters are missing

    Examples:
        >>> from mobidic import load_config, load_gisdata, Simulation
        >>> from mobidic.preprocessing.hyetograph import HyetographGenerator
        >>> from datetime import datetime
        >>> from pathlib import Path
        >>>
        >>> # Load configuration
        >>> config_file = Path("basin_hyetograph.yaml")
        >>> config = load_config(config_file)
        >>> gisdata = load_gisdata(config.paths.gisdata, config.paths.network)
        >>>
        >>> # Generate hyetograph
        >>> forcing = HyetographGenerator.from_config(
        ...     config=config,
        ...     base_path=config_file.parent,
        ...     start_time=datetime(2000, 1, 1)
        ... )
        >>>
        >>> # Run simulation
        >>> sim = Simulation(gisdata, forcing, config)
        >>> results = sim.run(forcing.start_date, forcing.end_date)

    Notes:
        - Automatically resamples IDF parameters to DEM grid
        - Uses all hyetograph parameters from config (duration, timestep, method, ka)
        - Output path read from config.paths.hyetograph
        - Creates metadata from config basin and hyetograph sections
        - Returns MeteoRaster ready for simulation with proper date range
    """
    # Import here to avoid circular dependency

    base_path = Path(base_path)

    # Validate config has hyetograph section
    if not hasattr(config, "hyetograph"):
        raise AttributeError("Configuration does not have a 'hyetograph' section")

    # Validate config has paths.hyetograph
    if not hasattr(config.paths, "hyetograph") or config.paths.hyetograph is None:
        raise ValueError(
            "Configuration must specify 'paths.hyetograph' for the output NetCDF file. "
            "Add 'hyetograph: path/to/output.nc' to the 'paths' section in the config file."
        )

    hyeto_config = config.hyetograph

    # Resolve paths relative to base_path
    a_raster_path = base_path / hyeto_config.a_raster
    n_raster_path = base_path / hyeto_config.n_raster
    k_raster_path = base_path / hyeto_config.k_raster
    ref_raster_path = base_path / config.raster_files.dtm
    output_path = base_path / config.paths.hyetograph

    # Ensure output directory exists
    output_path.parent.mkdir(parents=True, exist_ok=True)

    logger.info("Generating hyetograph forcing from configuration")
    logger.debug(f"  Base path: {base_path}")
    logger.debug(f"  Duration: {hyeto_config.duration_hours} hours")
    logger.debug(f"  Timestep: {hyeto_config.timestep_hours} hour(s)")
    logger.debug(f"  Method: {hyeto_config.hyetograph_type}")
    logger.debug(f"  Output: {output_path}")

    # Create generator with IDF parameters resampled to DEM grid
    generator = cls.from_rasters(
        a_raster=a_raster_path,
        n_raster=n_raster_path,
        k_raster=k_raster_path,
        ka=hyeto_config.ka,
        ref_raster=ref_raster_path,
    )

    # Prepare metadata from config
    add_metadata = {
        "hyetograph_method": hyeto_config.hyetograph_type,
        "duration_hours": hyeto_config.duration_hours,
        "timestep_hours": hyeto_config.timestep_hours,
        "areal_reduction_factor": hyeto_config.ka,
        "k_raster": str(hyeto_config.k_raster),
    }

    # Add basin metadata if available
    if hasattr(config, "basin"):
        if hasattr(config.basin, "id") and config.basin.id:
            add_metadata["basin"] = config.basin.id
        if hasattr(config.basin, "paramset_id") and config.basin.paramset_id:
            add_metadata["scenario"] = config.basin.paramset_id

    # Generate forcing using generate_forcing method
    forcing = generator.generate_forcing(
        duration_hours=hyeto_config.duration_hours,
        start_time=start_time,
        output_path=output_path,
        method=hyeto_config.hyetograph_type,
        timestep_hours=hyeto_config.timestep_hours,
        add_metadata=add_metadata,
        preload=preload,
    )

    return forcing

from_rasters(a_raster, n_raster, k_raster, ka=1.0, ref_raster=None) classmethod

Create HyetographGenerator by loading IDF parameters from raster files.

If a reference raster is provided, the IDF parameters will be resampled to match its extent, resolution, and CRS using nearest neighbor interpolation. This is the typical workflow when the IDF rasters have different resolution than the model grid (e.g., DEM).

Parameters:

Name Type Description Default
a_raster str | Path

Path to raster file containing IDF ‘a’ parameter

 required
n_raster str | Path

Path to raster file containing IDF ‘n’ parameter

 required
k_raster str | Path

Path to raster file containing IDF ‘k’ parameter

 required
ka float

Areal reduction factor (ARF) coefficient (default: 1.0)

 1.0
ref_raster str | Path | None

Optional path to reference raster (e.g., DEM) for resampling. If provided, IDF parameters will be resampled to match this grid.

 None

Returns:

Type Description
HyetographGenerator

HyetographGenerator instance

Examples:

>>> # Without resampling (IDF rasters already aligned)
>>> generator = HyetographGenerator.from_rasters(
...     a_raster="idf/a.tif",
...     n_raster="idf/n.tif",
...     k_raster="idf/k30.tif",
...     ka=0.8
... )
>>>
>>> # With resampling to DEM grid
>>> generator = HyetographGenerator.from_rasters(
...     a_raster="idf/a.tif",
...     n_raster="idf/n.tif",
...     k_raster="idf/k30.tif",
...     ka=0.8,
...     ref_raster="dem.tif"
... )
Source code in mobidic/preprocessing/hyetograph.py 
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
@classmethod
def from_rasters(
    cls,
    a_raster: str | Path,
    n_raster: str | Path,
    k_raster: str | Path,
    ka: float = 1.0,
    ref_raster: str | Path | None = None,
) -> "HyetographGenerator":
    """Create HyetographGenerator by loading IDF parameters from raster files.

    If a reference raster is provided, the IDF parameters will be resampled
    to match its extent, resolution, and CRS using nearest neighbor interpolation.
    This is the typical workflow when the IDF rasters have different resolution
    than the model grid (e.g., DEM).

    Args:
        a_raster: Path to raster file containing IDF 'a' parameter
        n_raster: Path to raster file containing IDF 'n' parameter
        k_raster: Path to raster file containing IDF 'k' parameter
        ka: Areal reduction factor (ARF) coefficient (default: 1.0)
        ref_raster: Optional path to reference raster (e.g., DEM) for resampling.
            If provided, IDF parameters will be resampled to match this grid.

    Returns:
        HyetographGenerator instance

    Examples:
        >>> # Without resampling (IDF rasters already aligned)
        >>> generator = HyetographGenerator.from_rasters(
        ...     a_raster="idf/a.tif",
        ...     n_raster="idf/n.tif",
        ...     k_raster="idf/k30.tif",
        ...     ka=0.8
        ... )
        >>>
        >>> # With resampling to DEM grid
        >>> generator = HyetographGenerator.from_rasters(
        ...     a_raster="idf/a.tif",
        ...     n_raster="idf/n.tif",
        ...     k_raster="idf/k30.tif",
        ...     ka=0.8,
        ...     ref_raster="dem.tif"
        ... )
    """
    if ref_raster is not None:
        idf_params = read_idf_parameters_resampled(a_raster, n_raster, k_raster, ref_raster)
    else:
        idf_params = read_idf_parameters(a_raster, n_raster, k_raster)
    return cls(idf_params, ka=ka)

generate(duration_hours, start_time, method='chicago_decreasing', timestep_hours=1)

Generate hyetograph precipitation time series.

Parameters:

Name Type Description Default
duration_hours int

Total duration of the hyetograph in hours

 required
start_time datetime

Start datetime for the hyetograph

 required
method Literal['chicago_decreasing']

Hyetograph construction method. Currently only “chicago_decreasing” is implemented.

 'chicago_decreasing'
timestep_hours int

Time step in hours (default: 1)

 1

Returns:

Type Description
tuple[list[datetime], ndarray]

Tuple of (times, precipitation) where: - times: List of datetime objects for each timestep - precipitation: 3D array (time, y, x) of precipitation (mm/h)

Raises:

Type Description
ValueError

If method is not supported
Notes
  • The Chicago decreasing method generates only the falling part of the Chicago hyetograph (after the peak)
  • Precipitation values are in mm/h (intensity)
  • NaN values in IDF parameters propagate to output
  • If NaN values are > 90% in any parameter, a ValueError is raised during initialization of the HyetographGenerator
Source code in mobidic/preprocessing/hyetograph.py 
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
def generate(
    self,
    duration_hours: int,
    start_time: datetime,
    method: Literal["chicago_decreasing"] = "chicago_decreasing",
    timestep_hours: int = 1,
) -> tuple[list[datetime], np.ndarray]:
    """Generate hyetograph precipitation time series.

    Args:
        duration_hours: Total duration of the hyetograph in hours
        start_time: Start datetime for the hyetograph
        method: Hyetograph construction method. Currently only
            "chicago_decreasing" is implemented.
        timestep_hours: Time step in hours (default: 1)

    Returns:
        Tuple of (times, precipitation) where:
            - times: List of datetime objects for each timestep
            - precipitation: 3D array (time, y, x) of precipitation (mm/h)

    Raises:
        ValueError: If method is not supported

    Notes:
        - The Chicago decreasing method generates only the falling part
          of the Chicago hyetograph (after the peak)
        - Precipitation values are in mm/h (intensity)
        - NaN values in IDF parameters propagate to output
        - If NaN values are > 90% in any parameter, a ValueError is raised during
          initialization of the HyetographGenerator
    """
    if method != "chicago_decreasing":
        raise ValueError(f"Unsupported hyetograph method: {method}. Only 'chicago_decreasing' is implemented.")

    logger.info(
        f"Generating {method} hyetograph: duration={duration_hours}h, "
        f"timestep={timestep_hours}h, start={start_time}"
    )

    return self._chicago_decreasing(duration_hours, start_time, timestep_hours)

generate_forcing(duration_hours, start_time, output_path, method='chicago_decreasing', timestep_hours=1, add_metadata=None, preload=True)

Generate hyetograph and return as MeteoRaster ready for simulation.

Convenience method that combines generate(), to_netcdf(), and MeteoRaster.from_netcdf() into a single call. This simplifies the workflow for design storm simulations.

Parameters:

Name Type Description Default
duration_hours int

Total duration of the hyetograph in hours

 required
start_time datetime

Start datetime for the hyetograph

 required
output_path str | Path

Path for output NetCDF file

 required
method Literal['chicago_decreasing']

Hyetograph construction method (default: “chicago_decreasing”)

 'chicago_decreasing'
timestep_hours int

Time step in hours (default: 1)

 1
add_metadata dict | None

Optional dictionary of additional global attributes

 None
preload bool

If True, preload all data into memory for fast access (default: True, recommended for normal use)

 True

Returns:

Type Description
MeteoRaster

MeteoRaster object ready for use in Simulation

Examples:

>>> # Create generator from IDF rasters
>>> generator = HyetographGenerator.from_rasters(
...     a_raster="idf/a.tif",
...     n_raster="idf/n.tif",
...     k_raster="idf/k30.tif",
...     ka=0.8,
...     ref_raster="dem.tif"
... )
>>>
>>> # Generate forcing and get MeteoRaster in one call
>>> forcing = generator.generate_forcing(
...     duration_hours=48,
...     start_time=datetime(2023, 11, 1),
...     output_path="design_hyetograph.nc",
...     add_metadata={"return_period": "30 years"}
... )
>>>
>>> # Use directly in simulation
>>> sim = Simulation(gisdata, forcing, config)
>>> results = sim.run(forcing.start_date, forcing.end_date)
Notes
  • This method is equivalent to calling generate(), to_netcdf(), and MeteoRaster.from_netcdf() sequentially
  • The NetCDF file is still created at output_path for later use
Source code in mobidic/preprocessing/hyetograph.py 
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
def generate_forcing(
    self,
    duration_hours: int,
    start_time: datetime,
    output_path: str | Path,
    method: Literal["chicago_decreasing"] = "chicago_decreasing",
    timestep_hours: int = 1,
    add_metadata: dict | None = None,
    preload: bool = True,
) -> "MeteoRaster":
    """Generate hyetograph and return as MeteoRaster ready for simulation.

    Convenience method that combines generate(), to_netcdf(), and
    MeteoRaster.from_netcdf() into a single call. This simplifies the
    workflow for design storm simulations.

    Args:
        duration_hours: Total duration of the hyetograph in hours
        start_time: Start datetime for the hyetograph
        output_path: Path for output NetCDF file
        method: Hyetograph construction method (default: "chicago_decreasing")
        timestep_hours: Time step in hours (default: 1)
        add_metadata: Optional dictionary of additional global attributes
        preload: If True, preload all data into memory for fast access
            (default: True, recommended for normal use)

    Returns:
        MeteoRaster object ready for use in Simulation

    Examples:
        >>> # Create generator from IDF rasters
        >>> generator = HyetographGenerator.from_rasters(
        ...     a_raster="idf/a.tif",
        ...     n_raster="idf/n.tif",
        ...     k_raster="idf/k30.tif",
        ...     ka=0.8,
        ...     ref_raster="dem.tif"
        ... )
        >>>
        >>> # Generate forcing and get MeteoRaster in one call
        >>> forcing = generator.generate_forcing(
        ...     duration_hours=48,
        ...     start_time=datetime(2023, 11, 1),
        ...     output_path="design_hyetograph.nc",
        ...     add_metadata={"return_period": "30 years"}
        ... )
        >>>
        >>> # Use directly in simulation
        >>> sim = Simulation(gisdata, forcing, config)
        >>> results = sim.run(forcing.start_date, forcing.end_date)

    Notes:
        - This method is equivalent to calling generate(), to_netcdf(),
          and MeteoRaster.from_netcdf() sequentially
        - The NetCDF file is still created at output_path for later use
    """
    # Import here to avoid circular dependency
    from mobidic.preprocessing.meteo_raster import MeteoRaster

    # Generate hyetograph
    times, precipitation = self.generate(
        duration_hours=duration_hours,
        start_time=start_time,
        method=method,
        timestep_hours=timestep_hours,
    )

    # Save to NetCDF
    self.to_netcdf(
        output_path=output_path,
        times=times,
        precipitation=precipitation,
        add_metadata=add_metadata,
    )

    # Load as MeteoRaster
    forcing = MeteoRaster.from_netcdf(output_path, preload=preload)

    logger.success(f"Forcing data ready for simulation: {forcing.start_date} to {forcing.end_date}")

    return forcing

to_netcdf(output_path, times, precipitation, add_metadata=None)

Export hyetograph to CF-compliant NetCDF file.

Creates a NetCDF file compatible with MeteoRaster.from_netcdf() for use as meteorological forcing in MOBIDIC simulations.

Parameters:

Name Type Description Default
output_path str | Path

Path for output NetCDF file

 required
times list[datetime]

List of datetime objects for each timestep

 required
precipitation ndarray

3D array (time, y, x) of precipitation [mm/h]

 required
add_metadata dict | None

Optional dictionary of additional global attributes

 None
Notes
  • Output follows CF-1.12 conventions
  • Includes CRS information as grid mapping variable
  • Precipitation units are mm/h (compatible with MeteoRaster)

Examples:

>>> times, precip = generator.generate(48, datetime(2023, 11, 1))
>>> generator.to_netcdf(
...     "design_storm.nc",
...     times=times,
...     precipitation=precip,
...     add_metadata={"return_period": "30 years"}
... )
Source code in mobidic/preprocessing/hyetograph.py 
 934
 935
 936
 937
 938
 939
 940
 941
 942
 943
 944
 945
 946
 947
 948
 949
 950
 951
 952
 953
 954
 955
 956
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
def to_netcdf(
    self,
    output_path: str | Path,
    times: list[datetime],
    precipitation: np.ndarray,
    add_metadata: dict | None = None,
) -> None:
    """Export hyetograph to CF-compliant NetCDF file.

    Creates a NetCDF file compatible with MeteoRaster.from_netcdf() for
    use as meteorological forcing in MOBIDIC simulations.

    Args:
        output_path: Path for output NetCDF file
        times: List of datetime objects for each timestep
        precipitation: 3D array (time, y, x) of precipitation [mm/h]
        add_metadata: Optional dictionary of additional global attributes

    Notes:
        - Output follows CF-1.12 conventions
        - Includes CRS information as grid mapping variable
        - Precipitation units are mm/h (compatible with MeteoRaster)

    Examples:
        >>> times, precip = generator.generate(48, datetime(2023, 11, 1))
        >>> generator.to_netcdf(
        ...     "design_storm.nc",
        ...     times=times,
        ...     precipitation=precip,
        ...     add_metadata={"return_period": "30 years"}
        ... )
    """
    output_path = Path(output_path)
    logger.info(f"Writing hyetograph to NetCDF: {output_path}")

    # Build coordinate arrays
    nrows, ncols = self.idf_params.shape
    cellsize = self.idf_params.cellsize

    # X coordinates (cell centers, west to east)
    x = np.arange(ncols) * cellsize + self.idf_params.xllcorner

    # Y coordinates (cell centers, south to north)
    # Note: Data is stored with y increasing (south to north) after flipud in grid_to_matrix
    y = np.arange(nrows) * cellsize + self.idf_params.yllcorner

    # Create xarray Dataset
    ds = xr.Dataset(
        data_vars={
            "precipitation": (
                ["time", "y", "x"],
                precipitation,
                {
                    "units": "mm h-1",
                    "long_name": "Precipitation rate",
                    "grid_mapping": "crs",
                },
            ),
        },
        coords={
            "time": times,
            "y": y,
            "x": x,
        },
        attrs={
            "Conventions": "CF-1.12",
            "title": "Synthetic hyetograph from IDF parameters",
            "source": f"MOBIDICpy version {__version__}",
            "history": f"Created {datetime.now().isoformat()}",
            "hyetograph_method": "chicago_decreasing",
            "areal_reduction_factor": self.ka,
        },
    )

    # Add coordinate attributes
    ds.x.attrs = {
        "units": "m",
        "long_name": "x coordinate",
        "standard_name": "projection_x_coordinate",
    }
    ds.y.attrs = {
        "units": "m",
        "long_name": "y coordinate",
        "standard_name": "projection_y_coordinate",
    }
    ds.time.attrs = {
        "long_name": "time",
        "standard_name": "time",
    }

    # Add CRS variable
    if self.idf_params.crs is not None:
        crs_attrs = crs_to_cf_attrs(self.idf_params.crs)
        # Create scalar CRS variable
        ds["crs"] = xr.DataArray(
            0,  # Scalar placeholder
            attrs=crs_attrs,
        )

    # Add custom metadata
    if add_metadata:
        ds.attrs.update(add_metadata)

    # Write to NetCDF with compression
    encoding = {
        "precipitation": {
            "zlib": True,
            "complevel": 4,
            "dtype": "float32",
        },
        "time": {"dtype": "float64"},
    }

    ds.to_netcdf(output_path, encoding=encoding)

    logger.success(f"Hyetograph written to: {output_path} ({precipitation.shape[0]} timesteps)")

Output classes

MeteoWriter

Writer for saving interpolated meteorological grids from station-based simulations.

Writer for interpolated meteorological data grids.

This class stores interpolated meteorological grids during simulation and writes them to a single CF-compliant NetCDF file at the end.

Unlike StateWriter which uses incremental writing with chunking, this writer buffers all data in memory and writes once. This is more efficient for meteo data which is typically much smaller than state data.

Parameters:

Name Type Description Default
output_path str | Path

Path to output NetCDF file (will be created/overwritten)

 required
grid_metadata dict

Dictionary with grid metadata (shape, resolution, crs, etc.)

 required
variables list[str]

List of meteorological variable names to save (e.g., [‘precipitation’, ‘temperature’])

 required
add_metadata dict | None

Additional global attributes (optional)

 None

Examples:

>>> # Create writer for precipitation and temperature
>>> writer = MeteoWriter(
...     "meteo_forcing.nc",
...     metadata,
...     variables=['precipitation', 'temperature']
... )
>>> for step in range(num_steps):
...     writer.append(current_time, precipitation=precip_grid, temperature=temp_grid)
>>> writer.close()  # Writes all data to file
Source code in mobidic/io/meteo.py 
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
class MeteoWriter:
    """
    Writer for interpolated meteorological data grids.

    This class stores interpolated meteorological grids during simulation
    and writes them to a single CF-compliant NetCDF file at the end.

    Unlike StateWriter which uses incremental writing with chunking, this writer
    buffers all data in memory and writes once. This is more efficient for meteo
    data which is typically much smaller than state data.

    Args:
        output_path: Path to output NetCDF file (will be created/overwritten)
        grid_metadata: Dictionary with grid metadata (shape, resolution, crs, etc.)
        variables: List of meteorological variable names to save (e.g., ['precipitation', 'temperature'])
        add_metadata: Additional global attributes (optional)

    Examples:
        >>> # Create writer for precipitation and temperature
        >>> writer = MeteoWriter(
        ...     "meteo_forcing.nc",
        ...     metadata,
        ...     variables=['precipitation', 'temperature']
        ... )
        >>> for step in range(num_steps):
        ...     writer.append(current_time, precipitation=precip_grid, temperature=temp_grid)
        >>> writer.close()  # Writes all data to file
    """

    def __init__(
        self,
        output_path: str | Path,
        grid_metadata: dict,
        variables: list[str],
        add_metadata: dict | None = None,
    ):
        """Initialize the meteo writer."""
        self.output_path = Path(output_path)
        self.output_path.parent.mkdir(parents=True, exist_ok=True)

        self.grid_metadata = grid_metadata
        self.variables = variables
        self.add_metadata = add_metadata or {}

        # Get grid dimensions
        self.nrows, self.ncols = grid_metadata["shape"]
        resolution = grid_metadata["resolution"]
        xllcorner = grid_metadata["xllcorner"]
        yllcorner = grid_metadata["yllcorner"]

        # Create coordinate arrays
        self.x = xllcorner + np.arange(self.ncols) * resolution[0]
        self.y = yllcorner + np.arange(self.nrows) * resolution[1]

        # Initialize buffers for data
        self.data_buffer = {var: [] for var in variables}
        self.time_buffer = []

        logger.debug(f"MeteoWriter initialized: {self.output_path}, variables={variables}")

    def append(self, time: datetime, **grids: np.ndarray) -> None:
        """
        Append meteorological grids for a timestep.

        Args:
            time: Timestamp for this data
            **grids: Keyword arguments for each variable grid (e.g., precipitation=precip_grid).
                Each grid should be a 2D numpy array with shape (nrows, ncols).

        Raises:
            ValueError: If a required variable is missing or has wrong shape

        Examples:
            >>> writer.append(current_time, precipitation=precip_grid, temperature=temp_grid)
        """
        # Validate that all required variables are provided
        for var in self.variables:
            if var not in grids:
                raise ValueError(f"Missing required variable '{var}' in append call")

            grid = grids[var]
            if grid.shape != (self.nrows, self.ncols):
                raise ValueError(
                    f"Variable '{var}' has incorrect shape {grid.shape}, expected {(self.nrows, self.ncols)}"
                )

        # Append data to buffers
        self.time_buffer.append(time)
        for var in self.variables:
            self.data_buffer[var].append(grids[var].copy())

    def close(self) -> None:
        """
        Write all buffered data to NetCDF file and close writer.

        This creates a CF-1.12 compliant NetCDF file with all meteorological
        variables and metadata.
        """
        if len(self.time_buffer) == 0:
            logger.warning("No data to write (buffer is empty)")
            return

        logger.info(f"Writing interpolated meteorological data to NetCDF: {self.output_path}")

        # Convert buffers to numpy arrays
        n_times = len(self.time_buffer)
        data_vars = {}

        for var in self.variables:
            # Stack list of 2D arrays into 3D array (time, y, x)
            data_array = np.stack(self.data_buffer[var], axis=0)

            # Convert units from m/s to mm/h for precipitation and PET
            if var in ["precipitation", "pet"]:
                data_array = data_array * 1000.0 * 3600.0

            # Create DataArray with metadata
            data_vars[var] = (
                ["time", "y", "x"],
                data_array,
                {
                    "long_name": _get_variable_longname(var),
                    "units": _get_variable_units(var),
                    "grid_mapping": "crs",
                },
            )

        # Create xarray Dataset
        ds = xr.Dataset(
            data_vars,
            coords={
                "time": (["time"], self.time_buffer, {"long_name": "time", "axis": "T"}),
                "x": (
                    ["x"],
                    self.x,
                    {
                        "standard_name": "projection_x_coordinate",
                        "long_name": "x coordinate of projection",
                        "units": "m",
                        "axis": "X",
                    },
                ),
                "y": (
                    ["y"],
                    self.y,
                    {
                        "standard_name": "projection_y_coordinate",
                        "long_name": "y coordinate of projection",
                        "units": "m",
                        "axis": "Y",
                    },
                ),
            },
        )

        # Add grid mapping variable for CRS (CF-1.12 compliance)
        ds["crs"] = ([], 0)
        ds["crs"].attrs = crs_to_cf_attrs(self.grid_metadata.get("crs"))

        # Add global attributes
        ds.attrs["title"] = "MOBIDIC meteorological forcing data"
        ds.attrs["source"] = "Meteo forcing data used in the simulation"
        ds.attrs["Conventions"] = "CF-1.12"
        ds.attrs["history"] = f"Created by MOBIDICpy v{__version__}"
        ds.attrs["date_created"] = datetime.now().isoformat()

        # Add user-provided metadata
        ds.attrs.update(self.add_metadata)

        # Set encoding for compression
        encoding = {}
        for var in self.variables:
            encoding[var] = {
                "zlib": True,
                "complevel": 4,
                "dtype": "float32",
                "_FillValue": np.nan,
            }

        # Write to NetCDF
        ds.to_netcdf(
            self.output_path,
            encoding=encoding,
            format="NETCDF4",
        )

        logger.success(
            f"Saved interpolated meteorological data: {len(self.variables)} variables, "
            f"{n_times} timesteps, {self.output_path}"
        )

        # Clear buffers to free memory
        self.data_buffer.clear()
        self.time_buffer.clear()

    def __enter__(self):
        """Context manager entry."""
        return self

    def __exit__(self, exc_type, exc_val, exc_tb):
        """Context manager exit - automatically close and write data."""
        if exc_type is None:
            # No exception - write data
            self.close()
        else:
            # Exception occurred - log and don't write
            logger.error(f"MeteoWriter exiting due to exception: {exc_val}")
        return False

__enter__()

Context manager entry.

Source code in mobidic/io/meteo.py 
211
212
213
def __enter__(self):
    """Context manager entry."""
    return self

__exit__(exc_type, exc_val, exc_tb)

Context manager exit - automatically close and write data.

Source code in mobidic/io/meteo.py 
215
216
217
218
219
220
221
222
223
def __exit__(self, exc_type, exc_val, exc_tb):
    """Context manager exit - automatically close and write data."""
    if exc_type is None:
        # No exception - write data
        self.close()
    else:
        # Exception occurred - log and don't write
        logger.error(f"MeteoWriter exiting due to exception: {exc_val}")
    return False

__init__(output_path, grid_metadata, variables, add_metadata=None)

Initialize the meteo writer.

Source code in mobidic/io/meteo.py 
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
def __init__(
    self,
    output_path: str | Path,
    grid_metadata: dict,
    variables: list[str],
    add_metadata: dict | None = None,
):
    """Initialize the meteo writer."""
    self.output_path = Path(output_path)
    self.output_path.parent.mkdir(parents=True, exist_ok=True)

    self.grid_metadata = grid_metadata
    self.variables = variables
    self.add_metadata = add_metadata or {}

    # Get grid dimensions
    self.nrows, self.ncols = grid_metadata["shape"]
    resolution = grid_metadata["resolution"]
    xllcorner = grid_metadata["xllcorner"]
    yllcorner = grid_metadata["yllcorner"]

    # Create coordinate arrays
    self.x = xllcorner + np.arange(self.ncols) * resolution[0]
    self.y = yllcorner + np.arange(self.nrows) * resolution[1]

    # Initialize buffers for data
    self.data_buffer = {var: [] for var in variables}
    self.time_buffer = []

    logger.debug(f"MeteoWriter initialized: {self.output_path}, variables={variables}")

append(time, **grids)

Append meteorological grids for a timestep.

Parameters:

Name Type Description Default
time datetime

Timestamp for this data

 required
**grids ndarray

Keyword arguments for each variable grid (e.g., precipitation=precip_grid). Each grid should be a 2D numpy array with shape (nrows, ncols).

 {}

Raises:

Type Description
ValueError

If a required variable is missing or has wrong shape

Examples:

>>> writer.append(current_time, precipitation=precip_grid, temperature=temp_grid)
Source code in mobidic/io/meteo.py 
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
def append(self, time: datetime, **grids: np.ndarray) -> None:
    """
    Append meteorological grids for a timestep.

    Args:
        time: Timestamp for this data
        **grids: Keyword arguments for each variable grid (e.g., precipitation=precip_grid).
            Each grid should be a 2D numpy array with shape (nrows, ncols).

    Raises:
        ValueError: If a required variable is missing or has wrong shape

    Examples:
        >>> writer.append(current_time, precipitation=precip_grid, temperature=temp_grid)
    """
    # Validate that all required variables are provided
    for var in self.variables:
        if var not in grids:
            raise ValueError(f"Missing required variable '{var}' in append call")

        grid = grids[var]
        if grid.shape != (self.nrows, self.ncols):
            raise ValueError(
                f"Variable '{var}' has incorrect shape {grid.shape}, expected {(self.nrows, self.ncols)}"
            )

    # Append data to buffers
    self.time_buffer.append(time)
    for var in self.variables:
        self.data_buffer[var].append(grids[var].copy())

close()

Write all buffered data to NetCDF file and close writer.

This creates a CF-1.12 compliant NetCDF file with all meteorological variables and metadata.

Source code in mobidic/io/meteo.py 
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
def close(self) -> None:
    """
    Write all buffered data to NetCDF file and close writer.

    This creates a CF-1.12 compliant NetCDF file with all meteorological
    variables and metadata.
    """
    if len(self.time_buffer) == 0:
        logger.warning("No data to write (buffer is empty)")
        return

    logger.info(f"Writing interpolated meteorological data to NetCDF: {self.output_path}")

    # Convert buffers to numpy arrays
    n_times = len(self.time_buffer)
    data_vars = {}

    for var in self.variables:
        # Stack list of 2D arrays into 3D array (time, y, x)
        data_array = np.stack(self.data_buffer[var], axis=0)

        # Convert units from m/s to mm/h for precipitation and PET
        if var in ["precipitation", "pet"]:
            data_array = data_array * 1000.0 * 3600.0

        # Create DataArray with metadata
        data_vars[var] = (
            ["time", "y", "x"],
            data_array,
            {
                "long_name": _get_variable_longname(var),
                "units": _get_variable_units(var),
                "grid_mapping": "crs",
            },
        )

    # Create xarray Dataset
    ds = xr.Dataset(
        data_vars,
        coords={
            "time": (["time"], self.time_buffer, {"long_name": "time", "axis": "T"}),
            "x": (
                ["x"],
                self.x,
                {
                    "standard_name": "projection_x_coordinate",
                    "long_name": "x coordinate of projection",
                    "units": "m",
                    "axis": "X",
                },
            ),
            "y": (
                ["y"],
                self.y,
                {
                    "standard_name": "projection_y_coordinate",
                    "long_name": "y coordinate of projection",
                    "units": "m",
                    "axis": "Y",
                },
            ),
        },
    )

    # Add grid mapping variable for CRS (CF-1.12 compliance)
    ds["crs"] = ([], 0)
    ds["crs"].attrs = crs_to_cf_attrs(self.grid_metadata.get("crs"))

    # Add global attributes
    ds.attrs["title"] = "MOBIDIC meteorological forcing data"
    ds.attrs["source"] = "Meteo forcing data used in the simulation"
    ds.attrs["Conventions"] = "CF-1.12"
    ds.attrs["history"] = f"Created by MOBIDICpy v{__version__}"
    ds.attrs["date_created"] = datetime.now().isoformat()

    # Add user-provided metadata
    ds.attrs.update(self.add_metadata)

    # Set encoding for compression
    encoding = {}
    for var in self.variables:
        encoding[var] = {
            "zlib": True,
            "complevel": 4,
            "dtype": "float32",
            "_FillValue": np.nan,
        }

    # Write to NetCDF
    ds.to_netcdf(
        self.output_path,
        encoding=encoding,
        format="NETCDF4",
    )

    logger.success(
        f"Saved interpolated meteorological data: {len(self.variables)} variables, "
        f"{n_times} timesteps, {self.output_path}"
    )

    # Clear buffers to free memory
    self.data_buffer.clear()
    self.time_buffer.clear()

Utility functions

convert_mat_to_netcdf

Direct conversion from MATLAB .mat format to CF-compliant NetCDF.

Convert MATLAB .mat meteorological data to NetCDF format.

This is a convenience function that combines loading from .mat and saving to NetCDF in a single call.

Parameters:

Name Type Description Default
mat_path str | Path

Path to input MATLAB .mat file

 required
output_path str | Path

Path to output NetCDF file

 required
compression_level int

Compression level (0-9, default 4)

 4
add_metadata dict[str, Any] | None

Additional global attributes for NetCDF file

 None

Examples:

>>> convert_mat_to_netcdf(
...     "examples/Arno/meteodata/meteodata.mat",
...     "examples/Arno/meteodata/meteodata.nc",
...     add_metadata={"basin": "Arno", "basin_id": "1234"}
... )
Source code in mobidic/preprocessing/meteo_preprocessing.py 
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
def convert_mat_to_netcdf(
    mat_path: str | Path,
    output_path: str | Path,
    compression_level: int = 4,
    add_metadata: dict[str, Any] | None = None,
) -> None:
    """Convert MATLAB .mat meteorological data to NetCDF format.

    This is a convenience function that combines loading from .mat and
    saving to NetCDF in a single call.

    Args:
        mat_path: Path to input MATLAB .mat file
        output_path: Path to output NetCDF file
        compression_level: Compression level (0-9, default 4)
        add_metadata: Additional global attributes for NetCDF file

    Examples:
        >>> convert_mat_to_netcdf(
        ...     "examples/Arno/meteodata/meteodata.mat",
        ...     "examples/Arno/meteodata/meteodata.nc",
        ...     add_metadata={"basin": "Arno", "basin_id": "1234"}
        ... )
    """
    logger.info("=" * 80)
    logger.info("MOBIDIC METEO PREPROCESSING")
    logger.info("=" * 80)

    # Load from MAT
    meteo_data = MeteoData.from_mat(mat_path)
    logger.info(f"Loaded {len(meteo_data.variables)} meteorological variables")
    logger.info(f"Date range: {meteo_data.start_date} to {meteo_data.end_date}")

    # Save to NetCDF
    meteo_data.to_netcdf(output_path, compression_level, add_metadata)

    logger.info("=" * 80)
    logger.info("METEO PREPROCESSING COMPLETE")
    logger.info("=" * 80)

Supported variables

The meteorological preprocessing module handles six standard variables:

Variable Description Units (NetCDF) MATLAB field MATLAB units
precipitation Rainfall and snowfall mm/h sp mm (cumulated)
temperature_min Minimum air temperature °C s_ta_min °C
temperature_max Maximum air temperature °C s_ta_max °C
humidity Relative humidity % s_ua %
wind_speed Wind speed m/s s_vv m/s
radiation Solar radiation W/m² s_ra W/m²

Automatic precipitation unit conversion

When loading from MATLAB .mat files, precipitation is automatically converted from mm (cumulated over the sampling interval) to mm/h. The conversion infers the timestep from the time array (dt = time[1] - time[0]) and divides precipitation values by the timestep in hours. This ensures consistency with raster forcing, which also uses mm/h.

Usage examples

Example 1: Quick conversion from MATLAB format

The fastest way to convert legacy MATLAB meteorological data:

from mobidic import convert_mat_to_netcdf

convert_mat_to_netcdf(
    "input/meteodata.mat",
    "output/meteodata.nc",
    compression_level=4,
    add_metadata={
        "basin": "Arno",
        "description": "Meteorological forcing for flood event",
    }
)

Example 2: Load, inspect, and save station data

When you need to examine data before using it in simulations:

from mobidic import MeteoData

# Load from MATLAB file
meteo = MeteoData.from_mat("input/meteodata.mat")

# Inspect loaded data
print(meteo)
print(f"Date range: {meteo.start_date} to {meteo.end_date}")
print(f"Variables: {meteo.variables}")

# Check station counts
for var_name, stations in meteo.stations.items():
    print(f"{var_name}: {len(stations)} stations")

# Save to NetCDF with metadata
meteo.to_netcdf(
    "output/meteodata.nc",
    compression_level=4,
    add_metadata={"basin": "Arno"}
)

Example 3: Work with existing NetCDF station data

Load and access station data from CF-compliant NetCDF files:

from mobidic import MeteoData

# Load from NetCDF
meteo = MeteoData.from_netcdf("output/meteodata.nc")

# Access station data for a specific variable
precip_stations = meteo.stations["precipitation"]
first_station = precip_stations[0]

print(f"Station code: {first_station['code']}")
print(f"Location: ({first_station['x']:.2f}, {first_station['y']:.2f})")
print(f"Elevation: {first_station['elevation']:.1f} m")
print(f"Data shape: {first_station['data'].shape}")

Example 4: Use pre-interpolated raster forcing

Load and run simulations with gridded meteorological data:

from mobidic import MeteoRaster, Simulation, load_gisdata, load_config

# Load raster forcing (default: preload into memory for best performance)
forcing = MeteoRaster.from_netcdf("meteo_raster.nc")

# Inspect raster data
print(f"Variables: {forcing.variables}")
print(f"Time range: {forcing.start_date} to {forcing.end_date}")
print(f"Grid shape: {forcing.grid_metadata['shape']}")
print(f"Resolution: {forcing.grid_metadata['resolution']:.1f} m")

# For very large datasets (>several GB), use lazy loading to save memory
forcing_lazy = MeteoRaster.from_netcdf("meteo_raster.nc", preload=False)

# Use in simulation (automatically detects raster mode)
config = load_config("basin.yaml")
gisdata = load_gisdata("gisdata.nc", "network.parquet")
sim = Simulation(gisdata, forcing, config)
results = sim.run("2023-01-01", "2023-12-31")

Example 5: Export interpolated data for reuse

Create reusable raster forcing from station-based simulations:

from mobidic import MeteoData, Simulation, load_gisdata, load_config

# Load station-based forcing
forcing = MeteoData.from_netcdf("meteo_stations.nc")

# Enable meteo forcing output in configuration
config = load_config("basin.yaml")
config.output_forcing_data.meteo_data = True

# Run simulation (exports forcing grids automatically)
gisdata = load_gisdata("gisdata.nc", "network.parquet")
sim = Simulation(gisdata, forcing, config)
results = sim.run("2023-01-01", "2023-12-31")

# Forcing data saved to: {output_dir}/meteo_forcing.nc
# Use this file as raster forcing in subsequent runs for faster performance

Example 6: Generate design storm hyetograph from configuration

The simplest workflow for design storm simulations reads all parameters from the YAML configuration:

from datetime import datetime
from pathlib import Path
from mobidic import load_config, load_gisdata, Simulation
from mobidic.preprocessing.hyetograph import HyetographGenerator

# Load configuration with hyetograph section
config_file = Path("basin_hyetograph.yaml")
config = load_config(config_file)
gisdata = load_gisdata(config.paths.gisdata, config.paths.network)

# Generate forcing from config - only start_time needed!
# All parameters (IDF rasters, duration, timestep, output path) read from config
forcing = HyetographGenerator.from_config(
    config=config,
    base_path=config_file.parent,
    start_time=datetime(2000, 1, 1)
)

# Run simulation
sim = Simulation(gisdata, forcing, config)
results = sim.run(forcing.start_date, forcing.end_date)

Required configuration (basin_hyetograph.yaml):

paths:
  hyetograph: output/design_storm.nc  # Output NetCDF path

hyetograph:
  a_raster: idf/a.tif        # IDF 'a' parameter raster
  n_raster: idf/n.tif        # IDF 'n' parameter raster
  k_raster: idf/k30.tif      # Return period factor raster (e.g., 30-year)
  duration_hours: 48         # Storm duration
  timestep_hours: 1          # Time step
  hyetograph_type: chicago_decreasing  # Hyetograph method
  ka: 0.8                    # Areal reduction factor

Example 7: Generate hyetograph with manual parameters

For more control over the hyetograph generation process:

from datetime import datetime
from mobidic import Simulation, load_gisdata, load_config
from mobidic.preprocessing.hyetograph import HyetographGenerator

# Create generator from IDF rasters
generator = HyetographGenerator.from_rasters(
    a_raster="idf/a.tif",
    n_raster="idf/n.tif",
    k_raster="idf/k30.tif",
    ka=0.8,                    # Areal reduction factor
    ref_raster="dem.tif"       # Reference grid for resampling
)

# Generate forcing and get MeteoRaster in one call
forcing = generator.generate_forcing(
    duration_hours=48,
    start_time=datetime(2023, 11, 1),
    output_path="design_storm.nc",
    method="chicago_decreasing",
    timestep_hours=1,
    add_metadata={"return_period": "30 years"}
)

# Use directly in simulation
config = load_config("basin.yaml")
gisdata = load_gisdata("gisdata.nc", "network.parquet")
sim = Simulation(gisdata, forcing, config)
results = sim.run(forcing.start_date, forcing.end_date)

Example 8: Advanced hyetograph workflow (full control)

For complete control over generation and export:

from datetime import datetime
from mobidic import MeteoRaster
from mobidic.preprocessing.hyetograph import HyetographGenerator

# Create generator
generator = HyetographGenerator.from_rasters(
    a_raster="idf/a.tif",
    n_raster="idf/n.tif",
    k_raster="idf/k30.tif",
    ka=0.8,
    ref_raster="dem.tif"
)

# Generate time series (returns times and precipitation arrays)
times, precipitation = generator.generate(
    duration_hours=48,
    start_time=datetime(2023, 11, 1),
    method="chicago_decreasing",
    timestep_hours=1,
)

# Inspect generated data
print(f"Time steps: {len(times)}")
print(f"Precipitation shape: {precipitation.shape}")  # (time, y, x)
print(f"Peak intensity: {precipitation[0].max():.2f} mm/h")
print(f"Total depth: {precipitation.sum(axis=0).max():.1f} mm")

# Save to NetCDF with custom metadata
generator.to_netcdf(
    "design_storm.nc",
    times=times,
    precipitation=precipitation,
    add_metadata={
        "return_period": "30 years",
        "event_type": "design_storm",
        "basin": "Arno"
    }
)

# Load as MeteoRaster for simulation
forcing = MeteoRaster.from_netcdf("design_storm.nc")

Data formats

Station data (MeteoData)

MATLAB format

Legacy MATLAB .mat files should contain struct arrays with these fields:

% Example structure for precipitation (variable name: sp)
sp(i).code      % Station identifier (string or number)
sp(i).est       % X coordinate (easting)
sp(i).nord      % Y coordinate (northing)
sp(i).quota     % Elevation (m)
sp(i).name      % Station name (optional)
sp(i).time      % Time vector (MATLAB datenum)
sp(i).dati      % Data vector (values)

Variable names in MATLAB: sp (precipitation), s_ta_min (temperature min), s_ta_max (temperature max), s_ua (humidity), s_vv (wind speed), s_ra (radiation).

Unit conversions during import

  • MATLAB datenum: Automatically converted to pandas datetime, accounting for the epoch difference (MATLAB: January 1, 0000; Unix: January 1, 1970).
  • Precipitation: Automatically converted from mm (cumulated over sampling interval) to mm/h by inferring the timestep from the time array.

Internal Python representation

MeteoData stores station data as a nested dictionary:

meteo.stations = {
    "precipitation": [
        {
            "code": "STATION_001",
            "x": 671234.5,              # Easting
            "y": 4821098.3,             # Northing
            "elevation": 342.0,         # meters
            "name": "Firenze",
            "time": pd.DatetimeIndex([...]),
            "data": np.array([...])     # Time series values
        },
        # ... additional stations
    ],
    "temperature_min": [...],
    "temperature_max": [...],
    # ... other variables
}

NetCDF format (CF-1.12 compliant)

Station data exported to NetCDF follows Climate and Forecast conventions:

Structure: 

Dimensions:
  - time: N timesteps
  - station_precipitation: M stations for precipitation
  - station_temperature_min: K stations for temperature_min
  - ... (one dimension per variable)

Coordinates:
  - time(time): datetime64[ns]
  - x_precipitation(station_precipitation): float64
  - y_precipitation(station_precipitation): float64
  - elevation_precipitation(station_precipitation): float32
  - station_code_precipitation(station_precipitation): str

Data variables:
  - precipitation(time, station_precipitation): float32 [mm/h]
  - temperature_min(time, station_temperature_min): float32 [°C]
  - ... (one variable per meteorological parameter)

Features:

  • CF-1.12 compliant metadata (standard_name, long_name, units)
  • Zlib compression (default level 4) for efficient storage
  • Custom global attributes supported via add_metadata
  • Precipitation stored in mm/h (consistent with raster forcing)

Raster data (MeteoRaster)

NetCDF format (CF-1.12 compliant)

Gridded forcing data follows standard CF conventions:

Required structure: 

Dimensions:
  - time: unlimited (temporal)
  - y: grid rows (north-south)
  - x: grid columns (east-west)

Coordinates:
  - time(time): datetime64[ns] with CF-compliant calendar
  - x(x): float64 (easting or longitude)
  - y(y): float64 (northing or latitude)

Data variables:
  - precipitation(time, y, x): float32 [mm/h]
  - pet(time, y, x): float32 [mm/h]
  - temperature(time, y, x): float32 [°C]
  - ... (additional variables as needed)

Grid mapping:
  - crs(): int32 with spatial_ref or crs_wkt attribute

Example inspection: 

import xarray as xr

ds = xr.open_dataset("meteo_raster.nc")
print(ds)
# <xarray.Dataset>
# Dimensions:         (time: 8761, y: 200, x: 300)
# Coordinates:
#   * time            (time) datetime64[ns] 2023-01-01 ... 2023-12-31
#   * x               (x) float64 600000.0 ... 629900.0
#   * y               (y) float64 4800000.0 ... 4820000.0
# Data variables:
#     precipitation   (time, y, x) float32 ...
#     pet             (time, y, x) float32 ...
#     temperature     (time, y, x) float32 ...
#     crs             () int32 ...

Requirements:

  • Grid must align with model grid (validated automatically)
  • Units: mm/h for precipitation, °C for temperature
  • Missing values: represented as NaN