Merge pull request #13 from EcoExtreML/read_forcing

Added read routines for the matlab forcing data

Author: BSchilperoort

PyStemmusScope/__init__.py

@@ -1,5 +1,7 @@
 """Documentation about PyStemmusScope"""
 import logging
+from . import forcing_io
+from . import iostreamer
 from .iostreamer import create_io_dir
 from .iostreamer import read_config
 

PyStemmusScope/forcing_io.py

@@ -0,0 +1,185 @@
+from pathlib import Path
+import hdf5storage
+import numpy as np
+import xarray as xr
+from . import variable_conversion as vc
+
+
+def _write_matlab_ascii(fname, data, ncols):
+    """Internal function to handle writing data in the Matlab ascii format. Equivalent
+    to `save([-], '-ascii')` in Matlab.
+
+    Args:
+        fname (path or str): Path to the file that should be written
+        data (np.array): Array with data to write to file
+        ncols (int, optional): The number of data columns, required to correctly format
+            the ascii file when writing multiple variables.
+    """
+    matlab_fmt = ' %14.7e'
+    multi_fmt = [matlab_fmt]*ncols
+    multi_fmt[0] = f' {multi_fmt[0]}'
+    np.savetxt(fname, data, multi_fmt)
+
+
+def read_forcing_data(forcing_file):
+    """Reads the forcing data from the provided netCDF file, and applies the required
+    unit conversions before returning the read data.
+
+    Args:
+        forcing_file (Path): Path to the netCDF file containing the forcing data
+
+    Returns:
+        dict: Dictionary containing the different variables required by STEMMUS_SCOPE
+            for the different forcing files.
+    """
+    ds_forcing = xr.open_dataset(forcing_file)
+
+    # remove the x and y coordinates from the data variables to make the numpy arrays 1D
+    ds_forcing = ds_forcing.squeeze(['x', 'y'])
+
+    data = {}
+    
+    # Expected time format is days (as floating point) since Jan 1st 00:00.
+    data['doy_float'] = (
+        ds_forcing['time'].dt.dayofyear - 1 +
+        ds_forcing['time'].dt.hour/24 +
+        ds_forcing['time'].dt.minute/60/24
+    )
+    data['year'] = ds_forcing['time'].dt.year.astype(float)
+
+    data['t_air_celcius'] = vc.kelvin_to_celcius(ds_forcing['Tair'])
+    data['psurf_hpa'] = vc.pa_to_hpa(ds_forcing['Psurf'])
+    data['co2_conv'] = vc.co2_molar_fraction_to_kg_per_m3(ds_forcing['CO2air'])
+    data['precip_conv'] = vc.precipitation_mm_s_to_cm_s(ds_forcing['Precip'])
+    data['lw_down'] = ds_forcing['LWdown']
+    data['sw_down'] = ds_forcing['SWdown']
+    data['wind_speed'] = ds_forcing['Wind']
+    data['rh'] = ds_forcing['RH']
+    data['vpd'] = ds_forcing['VPD']
+    data['lai'] = vc.mask_data(ds_forcing['LAI'], min_value=0.01)
+    data['ea'] = vc.kpa_to_hpa(vc.calculate_ea(data['t_air_celcius'], data['rh']))
+
+    # Load in non-timedependent variables
+    data['sitename'] = forcing_file.name.split('_')[0]
+    # Forcing data timestep size in seconds
+    time_delta = (ds_forcing.time.values[1] -
+                  ds_forcing.time.values[0]) / np.timedelta64(1, 's')
+    data['DELT'] = time_delta.astype(float)
+    data['total_timesteps'] = ds_forcing.time.size
+
+    data['latitude'] = ds_forcing['latitude'].values
+    data['longitude'] = ds_forcing['longitude'].values
+    data['elevation'] = ds_forcing['elevation'].values
+    data['IGBP_veg_long'] = ds_forcing['IGBP_veg_long'].values
+    data['reference_height'] = ds_forcing['reference_height'].values
+    data['canopy_height'] = ds_forcing['canopy_height'].values
+
+    return data
+
+
+def write_dat_files(data, input_dir):
+    """Fuction to write the single-data .dat files for the STEMMUS_SCOPE matlab model.
+
+    Args:
+        data (dict): Data dictionary generated by read_forcing
+        input_dir (path or str): Directory to which the different single-column .dat
+            files should be written to.
+    """
+    write_info = {
+        'doy_float': 't_.dat',
+        't_air_celcius': 'Ta_.dat',
+        'sw_down': 'Rin_.dat',
+        'lw_down': 'Rli_.dat',
+        'psurf_hpa': 'p_.dat',
+        'wind_speed': 'u_.dat',
+        'co2_conv': 'CO2_.dat',
+        'ea': 'ea_.dat',
+        'year': 'year_.dat'
+    }
+    for var, fname in write_info.items():
+        fpath = Path(input_dir) / fname
+        _write_matlab_ascii(fpath, data[var], ncols=1)
+
+
+def write_lai_file(data, fname):
+    """Function to write the ascii LAI_.dat file for STEMMUS_SCOPE.
+
+    Args:
+        data (dict): Dictionary containing the required variables. Generated by the
+            function `read_forcing_data`.
+        fname (Path): Full path, including filename, to which the file should be
+            written to.
+    """
+    lai_file_data = np.vstack([data['doy_float'], data['lai']]).T
+    _write_matlab_ascii(fname, lai_file_data, ncols=2)
+
+
+def write_meteo_file(data, fname):
+    """Function to write the ascii Mdata.txt meteo file for STEMMUS_SCOPE.
+
+    Args:
+        data (dict): Dictionary containing the required variables. Generated by the
+            function `read_forcing_data`.
+        fname (Path): Full path, including filename, to which the file should be
+            written to.
+    """
+    meteo_data_vars = ['doy_float', 't_air_celcius', 'rh',
+        'wind_speed', 'psurf_hpa', 'precip_conv', 'sw_down',
+        'lw_down', 'vpd', 'lai']
+    meteo_file_data = np.vstack([data[var] for var in meteo_data_vars]).T
+    _write_matlab_ascii(fname, meteo_file_data, ncols=len(meteo_data_vars))
+
+
+def prepare_global_variables(data, input_path, config):
+    """Function to read and calculate global variables for STEMMUS_SCOPE from the
+    forcing data. Data will be written to a Matlab binary file (v7.3), under the name
+    'forcing_globals.mat' in the specified input directory.
+
+    Args:
+        data (dict): Dictionary containing the required variables. Generated by the
+            function `read_forcing_data`.
+        input_path (Path): Path to which the file should be written to.
+        config (dict): The PyStemmusScope configuration dictionary.
+    """
+    if int(config['NumberOfTimeSteps']) > data['total_timesteps']:
+        total_duration = data['total_timesteps']
+    else:
+        total_duration = int(config['NumberOfTimeSteps'])
+
+    matfile_vars = ['latitude', 'longitude', 'elevation', 'IGBP_veg_long',
+                    'reference_height', 'canopy_height', 'DELT', 'sitename']
+    matfiledata = {key: data[key] for key in matfile_vars}
+
+    matfiledata['Dur_tot'] = float(total_duration) # Matlab expects a 'double'
+
+    hdf5storage.savemat(input_path / 'forcing_globals.mat', matfiledata, appendmat=False)
+
+
+def prepare_forcing(input_dir, forcing_file, config):
+    """Function to prepare the forcing files required by STEMMUS_SCOPE. The input
+        directory should be taken from the model configuration file.
+
+    Args:
+        input_dir (path or str): Path to the input directory that will be read by
+            STEMMUS_SCOPE.
+        forcing_file (path or str): Path to the netCDF forcing file that will be used
+            to generate the input data.
+        config (dict): The PyStemmusScope configuration dictionary.
+    """
+    input_path = Path(input_dir)
+    
+    # Read the required data from the forcing file into a dictionary
+    data = read_forcing_data(forcing_file)
+
+    # Write the single-column ascii '.dat' files to the input directory
+    write_dat_files(data, input_path)
+
+    # Write the two-column LAI_.dat file to the input directory.
+    write_lai_file(data, input_path / 'LAI_.dat')
+
+    # Write the multi-column Mdata.txt ascii file to the input directory
+    write_meteo_file(data, input_path / 'Mdata.txt')
+
+    # Write the remaining variables (without time dependency) to the matlab v7.3
+    #  file 'forcing_globals.mat'
+    prepare_global_variables(data, input_path, config)

PyStemmusScope/variable_conversion.py

@@ -0,0 +1,117 @@
+import numpy as np
+
+
+def calculate_ea(t_air_celcius, rh):
+    """Function that calculates the actual vapour pressure (kPa) from the
+    air temperature (degree Celcius) and relative humidity (%)
+
+    Args:
+        t_air_celcius: numpy array containing the air temperature in degrees C.
+
+        rh: numpy array of same shape as t_air_celcius, containing the relative humidity
+            as a percentage (e.g. ranging from 0 - 100).
+
+    Returns:
+        numpy array with the calculated actual vapor pressure
+    """
+    # Teten O. Über einige meteorologische Begriffe. Z. Geophys., 1930; 6. 297-309.
+    # Murray FW. On the computation of saturation vapor pressure,
+    #   J. Appl. Meteorol., 1967; 6, 203-204
+    es = 0.61078 * 10**(t_air_celcius*7.5 / (237.3+t_air_celcius))
+    return es * rh/100
+
+
+def kpa_to_hpa(pressure):
+    """Function to convert pressure data in kPa to hPa.
+
+    Args:
+        pressure (xr.DataArray): DataArray with pressure values in [kPa].
+
+    Returns:
+        xr.DataArray: Pressure data [hPa].
+
+    Raises:
+        ValueError if the input pressure is not in [kPa]
+    """
+    return pressure * 10
+
+
+def pa_to_hpa(pressure):
+    """Function to convert pressure data in Pa to hPa.
+
+    Args:
+        pressure (xr.DataArray): DataArray with pressure values in Pascal.
+
+    Returns:
+        xr.DataArray: Pressure data [hPa].
+
+    Raises:
+        ValueError if the input pressure is not in Pascal
+    """
+    if pressure.units == 'Pa':
+        return pressure / 100
+    raise ValueError('Input pressure is not in [Pa]')
+
+
+def kelvin_to_celcius(temperature):
+    """Function to convert temperature data in Kelvin to deg Celcius
+
+    Args:
+        temperature (xr.DataArray): DataArray with temperature values in Kelvin.
+
+    Returns:
+        xr.DataArray: Temperature in degrees Celcius.
+
+    Raises:
+        ValueError if the input temperature is not in Kelvin
+    """
+    if temperature.units=='K':
+        return temperature - 273.15
+    raise ValueError('Input temperature is not in [K]')
+
+
+def co2_molar_fraction_to_kg_per_m3(molar_fraction):
+    """Function to convert CO2 molar fraction [mol cO2/mol air] to CO2
+    concentration in [kg CO2 / m3 air]
+
+    Note: the density of air [kg/m3] used for the calculation is assumed to be constant
+    here, but will vary depending on the air pressure and air temperature.
+
+    Args:
+        molar_fraction (float, np.array): CO2 concentration as molar fraction
+
+    Returns:
+        Same as input: CO2 concentration in [kg CO2 / m3 air]
+    """
+    molecular_weight_co2 = 44.01 # [kg/mol]
+    avg_density_air = 1.292 # [kg/m3]
+    avg_molar_mass_air = 28.9647 # [kg/mol]
+    molar_density_air = avg_molar_mass_air / avg_density_air # [m3/mol]
+    return molar_fraction * molecular_weight_co2 / molar_density_air
+
+
+def precipitation_mm_s_to_cm_s(precipitation_rate):
+    """Function to convert precipitation rate in [kg/m2/s] or [mm/s] to [cm/s]
+
+    Args:
+        precipitation_rate (float, np.array): Precipitation rate in [kg/m2/s] or [mm/s].
+
+    Returns:
+        Same as input: Precipitation rate in [cm/s]
+    """
+    return precipitation_rate/10
+
+
+def mask_data(data, min_value=None, max_value=None):
+    """Function to apply a mask to data.
+
+    Args:
+        data (np.array): Array containing the original data.
+        condition (np.array): Boolean array to mask the data with.
+        value (float): Value to replaced the masked data with.
+
+    Returns:
+        np.array: Array where the values under the mask have been replaced with the
+            given value.
+    """
+    return np.clip(data, min_value, max_value)

notebooks/run_model_in_notebook.ipynb

@@ -17,7 +17,8 @@
     "import os\n",
     "from pathlib import Path\n",
     "import subprocess\n",
-    "import PyStemmusScope"
+    "from PyStemmusScope import forcing_io\n",
+    "from PyStemmusScope import iostreamer"
    ]
   },
   {
@@ -46,7 +47,7 @@
     "# path to config file\n",
     "path_to_config_file = \"path_to_config_file\"\n",
     "# Instantiate working directories handler from PyStemmusScope\n",
-    "config = PyStemmusScope.read_config(path_to_config_file)"
+    "config_template = iostreamer.read_config(path_to_config_file)"
    ]
   },
   {
@@ -79,7 +80,7 @@
     }
    ],
    "source": [
-    "config"
+    "config_template"
    ]
   },
   {
@@ -113,8 +114,8 @@
    ],
    "source": [
     "# edit config\n",
-    "config[\"NumberOfTimeSteps\"] = \"5\"\n",
-    "config"
+    "config_template[\"NumberOfTimeSteps\"] = \"5\"\n",
+    "config_template"
    ]
   },
   {
@@ -145,15 +146,20 @@
     "\n",
     "full_run = False\n",
     "if full_run:\n",
-    "    forcing_filenames_list = [file.name for file in Path(config[\"ForcingPath\"]).iterdir()]\n",
+    "    forcing_filenames_list = [file.name for file in Path(config_template[\"ForcingPath\"]).iterdir()]\n",
     "\n",
     "# empty dict to store path to config file and work directory for each station\n",
     "config_path_dict = {}\n",
     "input_dir_dict = {}\n",
     "for nc_file in forcing_filenames_list:\n",
-    "    input_dir, output_dir, config_path = PyStemmusScope.create_io_dir(nc_file, config)\n",
+    "    input_dir, output_dir, config_path = iostreamer.create_io_dir(nc_file, config_template)\n",
     "    config_path_dict[nc_file] = config_path\n",
     "    input_dir_dict[nc_file] = input_dir\n",
+    "    \n",
+    "    config_run = iostreamer.read_config(config_path)\n",
+    "\n",
+    "    forcing_io.prepare_forcing(input_dir, Path(config_run[\"ForcingPath\"]) / nc_file, config_run)\n",
+    "    \n",
     "    print(input_dir, output_dir, config_path)"
    ]
   },
@@ -277,7 +283,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.12"
+   "version": "3.9.13"
   }
  },
  "nbformat": 4,