Users intending to complete the vector-based workflow on a personal or workstation machine are encouraged to consult the instructions for installing Python and the Windows Subsystem for Linux page before proceeding.

The main structure of the vector-based workflow has been written in the Python language. Certain parts of the workflow require some external libraries and packages that must be installed in python.

1. Introduction

This page explains the workflow of the model setup for the MESH vector-based application. The model setup has three main parts which are explained in detail. The general structure of this workflow is shown in Figure 1. Note that the codes for generating the model setup are available from this Github page: https://github.com/MESH-Model/MESH-Scripts/tree/master/Model_Workflow/vector_based_workflow

Figure 1. MESH vector-based workflow

If users intend to setup the domain of interest from scratch based on any high-resolution Digital Elevation Model (DEM), rather than MeritHydro DEM, catchment, and river files, they can follow a workflow using GRASS GIS.

If a user intends to install and run easymore on a workstation machine, It is very important the required packages, before running easymore, be installed and satisfied appropriately. Otherwise, a user will receive a failure error. The list of required packages are as follow:

requirement

install_requires=[
        'numpy',
        'xarray',
        'pandas',
        'netCDF4',
        'datetime',
    ]
extras_require={
    ":python_version>'3.6'":['geopandas >= 0.8.1',
    'shapely',
    'pyshp',
    'pysheds',
    'gdal',
    'geovoronoi',
    'json5',
    'rasterio']
    }

2. Downloading Scripts & Input Data

2.1. Downloading the MESH Vector-Based Workflow Scripts

The scripts to produce a vector-based setup of MESH are available at: https://github.com/MESH-Model/MESH-Scripts. The best way to retrieve the necessary scripts is to use git. On Linux, run the following command to install git:

sudo apt-get install git-all

once git is installed, the user can use git clone to retrieve the workflow scripts:

git clone https://github.com/MESH-Model/MESH-Scripts.git

this will download the latest version of the workflow scripts to the directory from which the command is entered. To begin the setup workflow, navigate to MESH-Scripts/Model_Workflow/vector_based_workflow/0_control_files/ and create a duplicate of the control_BowAtBanff.txt file, renaming it as necessary to represent the domain of interest. Modify the settings within the control file as needed to represent the domain of interest. Once the settings have been changed, save the control file as control_active.txt – replacing the existing version of control_active.txt as needed.

2.2. Downloading the MERIT-Hydro Shapefiles

The workflow uses vectorized river flowlines and unit catchments based on the MERIT-Hydro raster hydrography dataset to represent the domain of interest in MESH.

The vectorized river networks and catchments are divided into large continental basins, each represented by a unique 'PFAF number'. The individual basins can be downloaded from the following links:

File	Link
Public Access	https://drive.google.com/drive/folders/1owkvZQBMZbvRv3V4Ff3xQPEgmAC48vJo
MERIT_Hydro_catchment Local Access	\\globalwater\MDR\PhysiographicData\MERIT_Hydro_gobal_basin\ MERIT_Hydro_basin_bugfixed
MERIT_Hydro_river Local Access	\\globalwater\MDR\PhysiographicData\MERIT_Hydro_gobal_basin\MERIT_Hydro_river

The map below is provided as a quick reference to help identify the location of different MERIT-Basins. For example, a user wishing to select the Bow River with the outlet at Banff as a subset should download the shapefile representing ‘PFAF 71’.

Figure 2. PFAF reference map

The following section will describe how to provide settings to the control file. Most of the user input is provided in this section. Once the control file settings are provided, the user can execute the processing scripts to create the drainage database.

The following sections do not currently accept input from the control file:

Producing zonal statistics
Reordering RDRSv2.1 climate forcing data
Model configuration

3. Providing Settings to the Control File

The purpose of the control file is to provide all inputs to the scripts in the vector-based workflow to eliminate the need to alter the workflow scripts themselves. In the current workflow version, the control file provides settings to the subbasin selection, network topology, and drainage database scripts.

Example Control File:

control_BowAtBanff.txt

# MESH vector-based workflow setting file.
# Characters '|' and '#' are used as separators to find the actual setting values. Any text behind '|' is assumed to be part of the setting value, unless preceded by '#'.
 
# Note on path specification
# If deviating from default paths, a full path must be specified. E.g. '/home/user/non-default/path'
 

# Modelling domain settings
root_path                      | /path/to/MESH-Scripts/Model_Workflow/                       # Root folder where data will be stored.
domain_name                    | BowAtBanff                                                  # Used as part of the root folder name for the prepared data.

## Subbasin Selection Settings
input_basin_path               | default                                                      # If 'default', uses 'root_path/shape_file/catchment'
input_basin_name               | cat_pfaf_71_MERIT_Hydro_v07_Basins_v01_bugfix1_WGS84.shp     # Name of the catchment shapefile from which to subset. (the basin containing your desired domain)
input_river_path               | default                                                      # If 'default', uses 'root_path/shape_file/river_network'
input_river_name               | rivEndoMERITpfaf_71_WGS84.shp                                # Name of river network shapefile from which to subset
subset_basin_outdir            | default                                                      # If 'default', uses 'root_path/shape_file/catchment (the river network containing your desired domain)
subset_river_outdir            | default                                                      # If 'default', uses 'root_path/shape_file/river_network'
case_name                      | BowAtBanff                                                   # This will define the name of output files  
target_segment                 | 71028585                                                     # The location of the outlet for the basin of interest. You may need to open the shapefile in QGIS or ArcGIS to identify this value from the attribute table.
                                                                                          # See: https://wiki.usask.ca/display/MESH/MESH+vector-based+workflow#MESHvectorbasedworkflow-2.1.MERITHydrosubset

## Network Topology Settings
# Shapefile settings - river network file
river_network_shp_path         | default                                                     # If 'default', uses 'root_path/shape_file/river_network'.
river_network_shp_name         | BowAtBanff_riv.shp                                          # Name of the river network shapefile. Requires extension '.shp'.
river_network_shp_segid        | COMID                                                       # Name of the segment ID column.
river_network_shp_downsegid    | NextDownID                                                  # Name of the downstream segment ID column.
river_network_shp_slope        | slope                                                       # Name of the slope column. Slope must be in in units [length/length].
river_network_shp_length       | lengthkm                                                    # Name of the segment length column. Length must be in units [m].
river_network_shp_outlet_id    | 71028585                                                    # ID of the most downstream segment in the river network (river network shapefile) ... 
                                                                                          # ... Note that this specifies a single segment ID value. It is not a column in the shapefile.
  
# Shapefile settings - catchment file
river_basin_shp_path           | default                                                     # If 'default', uses 'root_path/shape_file/catchment'.
river_basin_shp_name           | BowAtBanff_cat.shp                                          # Name of the routing subbasins shapefile needed for remapping. Requires extension '.shp'.
river_basin_shp_rm_hruid       | COMID                                                       # Name of the routing basin ID column.
river_basin_shp_area           | unitarea                                                    # Name of the catchment area column. Area must be in units [m^2]
river_basin_shp_hru_to_seg     | COMID                                                       # Name of the column that shows which river segment each HRU connects to.

# Experiment settings - routing
settings_routing_path          | default                                                  # If 'default', uses 'root_path/domain_[name]/topology'.
settings_routing_topology      | network_topology_BowBanff.nc                             # Name of the river network topology file.


## Drainage Database Settings
input_lc_zh_path		       | default                                     		      # If 'default', uses 'root_path/domain_[name]/zonalhist/
input_lc_zh_name		       | landsat_bow_stats_NA_NALCMS_2010_v2_land_cover_30m.csv    # Name of the input zonal staistics file. Can be either a .csv file (produced by GIS Tool) or a .shp file (produced in QGIS)
input_topo_path		           | default                                     		      # If 'default', uses 'root_path/domain_[name]/topology/
input_topo_name                | network_topology_BowBanff.nc                              # Name of the network topology file produced previously
DDB_output_dir			       | default                                     		      # If 'default', uses 'root_path/domain_[name]/drainagedatabase/
lc_type_prefix			       | frac_ 										  			  # Typically, use frac_ when using zonal statistics generated with GIS Tool (.csv). use NALCMS_ for QGIS zonal histogram (.shp) 
merit_basin_path	           | default                                     		      # If 'default', uses 'root_path/domain_[name]/subbasin/catchment/
merit_basin_name	           | BowAtBanff_cat.shp 	                                      # The parent shapefile used in subbasin selection (irrelevant if subbasin selection is not used)

# Climate forcing settings
remapping_out			   | default
source_nc_path             | ../forcing/
source_nc_name             | rdrsv2.1_1980.nc
var_lon                    | lon
var_lat                    | lat
var_time                   | time
format_list                | f4
fill_value_list            | -9999.00
var_names                  | RDRS_v2.1_A_PR0_SFC, RDRS_v2.1_P_FI_SFC, RDRS_v2.1_P_FB_SFC, RDRS_v2.1_P_TT_09944, RDRS_v2.1_P_UVC_09944, RDRS_v2.1_P_P0_SFC, RDRS_v2.1_P_HU_09944 # variables must be on 1 line and separated by a comma and a space.  # Default folder structure
# Example of the resulting folder structure in "root_path".
# New domains will go into their own folder.
 
- shape_file
   |
   |_ catchment
   |   |
   |   |_ cat_pfaf_71_MERIT_Hydro_v07_Basins_v01_bugfix1_WGS84.shp
   |   |
   |_ river_network
   |   |
   |   |_ rivEndoMERITpfaf_71_WGS84.shp
   |   |
- vector_based_workflow
   |
   |_ 0_control_files
   |
   |_ workflow_data
      |
      |_ domain_[name]
      |   |
      |   |_ drainagedatabase
      |   |   |
      |   |   |_ BowAtBanff_MESH_drainage_database.nc
      |   |   |
      |   |_ topology
      |   |   |
      |   |   |_ network_topology_BowBanff.nc
      |   |   |
      |   |_ zonalhist
      |   |   |
      |   |   |_ landsat_bow_stats_NA_NALCMS_2010_v2_land_cover_30m.csv

3.1. Root Path and Domain Name

The root path represents the folder within which all workflow scripts and data are stored. The domain name should represent the basin for which the user is setting up MESH. This setting will be used to name many of the workflow outputs.

control_BowAtBanff.txt

# Modeling domain settings
root_path                   | /path/to/MESH-Scripts/Model_Workflow/                     # Root folder where data will be stored.
domain_name                 | BowAtBanff                                                # Used as part of the root folder name for the prepared data.

By default, the root path should be the path to ‘Model_Workflow’ within the previously downloaded ‘MESH-Scripts’ folder. For example, this path may appear as follows:

on a local WSL installation:	/mnt/c/Users/<user>/Desktop/USASK/MESH-Scripts/Model_Workflow/
on Graham:	/project/6008034/<user>/MESH-Scripts/Model_Workflow/

3.2 Subbasin Selection Settings

The subbasin selection settings modify an optional script used to extract a smaller domain of interest from within one of the MERIT-Hydro Basins shapefiles. The script is optional, as some users may need to produce a setup for an entire MERIT-Hydro Basin. If this is the case, there is no need to modify the subbasin selection settings, and the subbasin selection script can be skipped entirely. Otherwise, proceed as follows:

control_BowAtBanff.txt

## Subbasin Selection Settings
input_basin_path           | default                                                      # If 'default', uses 'root_path/shape_file/catchment'
input_basin_name           | cat_pfaf_71_MERIT_Hydro_v07_Basins_v01_bugfix1_WGS84.shp     # Name of the catchment shapefile from which to subset. (the basin containing your desired domain)
input_river_path           | default                                                      # If 'default', uses 'root_path/shape_file/river_network'
input_river_name           | rivEndoMERITpfaf_71_WGS84.shp                                # Name of river network shapefile from which to subset
subset_basin_outdir        | default                                                      # If 'default', uses 'root_path/shape_file/catchment (the river network containing your desired domain)
subset_river_outdir        | default                                                      # If 'default', uses 'root_path/shape_file/river_network'
case_name                  | BowAtBanff                                                   # This will define the name of output files  
target_segment             | 71028585                                                     # The location of the outlet for the basin of interest. You may need to open the shapefile in QGIS or ArcGIS to identify this value from the attribute table.
                                                                                          # See: https://wiki.usask.ca/display/MESH/MESH+vector-based+workflow#MESHvectorbasedworkflow-2.1.MERITHydrosubset

Setting	Description
input_basin_path	Set the input basin path to the location of the MERIT-Hydro Basin catchment within which the domain of interest is contained. By default, the subbasin selection script will use ‘root_path/shape_file/catchment’. It is recommended to use the default path and save the downloaded shapefile here. However, it may be necessary to deviate from the default path. In that case, simply replace ‘default’ with the full path to the MERIT-Hydro Basins shapefile. For Graham users, MERIT-Hydro shapefiles are available at /project/6008034/Model_Output/MESH/NA_workflow/shape_file/catchment. Users may replace ‘default’ with this path to avoid waiting for a new download of the same file(s).
input_basin_name	The input_basin_name setting determines the name of the catchment shapefile used above. Typically, the MERIT-Hydro Basins catchment shapefiles are named “cat_pfaf_{}_MERIT_Hydro_v07_Basins_v01_bugfix1.shp” with ‘{}’ representing the number associated with each individual MERIT Basin. To select the Bow River at Banff, a user would supply cat_pfaf_71_MERIT_Hydro_v07_Basins_v01_bugfix1.shp to the input_basin_name setting. The MERIT-Hydro Basin files available on Graham have been named in a slightly different fashion, and users should review the naming convention there before setting the input_basin_name.
input_river_path	The input river path setting is analogous to the input_basin_path setting, except that it represents the MERIT-Hydro Basin river network shapefile. This setting can be specified in the same way, and files are similarly available on Graham in the directory /project/6008034/Model_Output/MESH/NA_workflow/shape_file/river_network
input_river_name	Analogous to input_basin_name, with the naming convention: “riv_pfaf_{}_MERIT_Hydro_v07_Basins_v01_bugfix1.shp”
subset_basin_outdir & subset_river_outdir	The ‘subset_basin_outdir’ and ‘subset_river_outdir’ settings define the output directories for the subbasin selection script. In other words, the output catchment and river network shapefiles from the subbasin selection script will be placed in the directories defined by these settings. By default, the outputs will be saved to root_path/shape_file/catchment and root_path/shape_file/river_network The output files will be named in the following convention: [domain_name]_cat.shp (e.g., BowAtBanff_cat.shp) [domain_name]_riv.shp (e.g., BowAtBanff_riv.shp)
case_name	case_name defines the name of the output files. Set this to be the same as the domain name. Setting is redundant and will be removed in the future.
target_segment	The target_segment setting represents the ID number of the outlet for the basin of interest. It may be necessary to open the catchment shapefile in GIS software (ArcGIS or QGIS) to identify this value. Inside GIS, locate and select the polygon containing the basin outlet (for the Bow River at Banff basin, this will be the polygon containing the Bow River at Banff – 05BB001 WSC gauging station). Open the attribute table and display only selected features. Use the ‘COMID’ of the selected feature as the setting for ‘target_segment’. The following figure shows the `COMID` of the outlet "71028585" for a basin of interest (e.g., Bow River at Banff). (Click image to enlarge) Figure 3. Finding the ID of a catchment near the outlet

3.3 Network Topology Settings

Creation of the network topology file is a requirement of the vector-based workflow. The network topology file will contain geophysical parameters such as unique indices of stream segments, routing parameters – such as channel length, slope, and basin area – and other required information from the catchment and river network shapefiles for the domain of interest. The network topology file is created with the EASYMORE package and is later converted into a drainage database to be used with MESH.

control_BowAtBanff.txt

## Network Topology Settings
# Shapefile settings - river network file
river_network_shp_path      | default                                                     # If 'default', uses 'root_path/shape_file/river_network'.
river_network_shp_name      | BowAtBanff_riv.shp                                          # Name of the river network shapefile. Requires extension '.shp'.
river_network_shp_segid     | COMID                                                       # Name of the segment ID column.
river_network_shp_downsegid | NextDownID                                                  # Name of the downstream segment ID column.
river_network_shp_slope     | slope                                                       # Name of the slope column. Slope must be in in units [length/length].
river_network_shp_length    | lengthkm                                                    # Name of the segment length column. Length must be in units [m].
river_network_shp_outlet_id | 71028585                                                    # ID of the most downstream segment in the river network (river network shapefile) ... 
                                                                                          # ... Note that this specifies a single segment ID value. It is not a column in the shapefile.
  
# Shapefile settings - catchment file
river_basin_shp_path        | default                                                     # If 'default', uses 'root_path/shape_file/catchment'.
river_basin_shp_name        | BowAtBanff_cat.shp                                          # Name of the routing subbasins shapefile needed for remapping. Requires extension '.shp'.
river_basin_shp_rm_hruid    | COMID                                                       # Name of the routing basin ID column.
river_basin_shp_area        | unitarea                                                    # Name of the catchment area column. Area must be in units [m^2]
river_basin_shp_hru_to_seg  | COMID                                                       # Name of the column that shows which river segment each HRU connects to.

# Experiment settings - routing
settings_routing_path          | default                                                  # If 'default', uses 'root_path/domain_[name]/topology'.
settings_routing_topology      | network_topology_BowBanff.nc                             # Name of the river network topology file.

The network topology settings are split into three sections: river network file, catchment file, and routing.

3.3.1. River Network File Settings

The settings in the river network file section define the name and location of the MERIT-Hydro river network shapefile to be used. This is either the file created with the subbasin selection script, or the full MERIT-Hydro river network file (if setting up for an entire PFAF).

Setting	Description
river_network_shp_segid	refers to the name of the column in the river network shapefile that represents the ID of each stream segment (typically, ‘COMID’).
river_network__shp_downsegid	refers to the name of the column in the river network shapefile that represents the ID of the river segment that is downstream from the current segment.
river_network_shp_slope	refers to the name of the column in the river network shapefile that represents the slope of the current river segment
river_network_shp_length	refers to the name of the column in the river network shapefile that represents the length in meters of the current river segment.
river_network_shp_outlet_id	refers to the most downstream river segment in the river network. Note that this is a single segment and not a column in the shapefile. If the basin of interest has multiple outlets (i.e., the user is setting up for an entire MERIT-Hydro Basin), the user must disable this setting by placing a ‘#’ before it. If the user fails to do so, the network topology file may be inaccurate or may not work at all.
river_order	refers to the name of the column in the river network shapefile that represents the order of each stream segment. Not used if include_IAK is set to False.
include_IAK	Use this setting to specify whether to include the river order column from the MERIT Hydro shapefile. If True, IAK will be included in the network topology file and therefore in the drainage database. Otherwise, set to False.

3.3.2. Catchment File Settings

The settings in the catchment file section define the name and location of the MERIT-Hydro catchment shapefile to be used. This is either the file created with the subbasin selection script, or the full MERIT-Hydro catchment file (if setting up for an entire PFAF).

Setting	Description
river_basin_shp_rm_hruid	refers to the name of the column in the catchment shapefile that represents the ID of each routing basi-n (typically, ‘COMID’).
river_basin_shp_area	refers to the name of the column in the catchment shapefile that represents the area in m² of each subbasin in the shapefile
river_basin_shp_hru_to_seg	refers to the name of the column in the catchment shapefile that indicates which river segment (in the river network shapefile) each subbasin connects to.

3.3.3. Routing Settings

The routing settings define the output of the network topology script. Leaving the path as default will place the resulting netcdf (.nc) file in the directory /root_path/vector_based_workflow/workflow_data/domain_name/topology/. The settings_routing_topology setting defines the name of the output file.

3.4 Drainage Database Settings

Though the settings can be applied now, a zonal statistics file must be created before proceeding to produce the drainage database.

By default, GIS Tool (used to generate zonal statistics) will place the required file in the directory /root_path/vector_based_workflow/workflow_data/[domain_name]/zonalhist/ and gives it the name landsat_[domain_name]_stats_NA_NALCMS_2010_v2_land_cover_30m.csv.

As zonal statistics settings are not currently included in the control file, it may be necessary to produce the zonal statistics before specifying its name and location if the default name and path are not used.

control_BowAtBanff.txt

## Drainage Database Settings
input_lc_zh_path		       | default                                     		      # If 'default', uses 'root_path/domain_[name]/zonalhist/
input_lc_zh_name		       | landsat_bow_stats_NA_NALCMS_2010_v2_land_cover_30m.csv    # Name of the input zonal staistics file. Can be either a .csv file (produced by GIS Tool) or a .shp file (produced in QGIS)
input_topo_path		           | default                                     		      # If 'default', uses 'root_path/domain_[name]/topology/
input_topo_name                | network_topology_BowBanff.nc                              # Name of the network topology file produced previously
DDB_output_dir			       | default                                     		      # If 'default', uses 'root_path/domain_[name]/drainagedatabase/
lc_type_prefix			       | frac_ 										  			  # Typically, use frac_ when using zonal statistics generated with GIS Tool (.csv). use NALCMS_ for QGIS zonal histogram (.shp) 
merit_basin_path	           | default                                     		      # If 'default', uses 'root_path/domain_[name]/subbasin/catchment/
merit_basin_name	           | BowAtBanff_cat.shp 	                                      # The parent shapefile used in subbasin selection (irrelevant if subbasin selection is not used)

To create the drainage database file, two main processes are carried out. First, the required variables for creating the drainage database (a.k.a. the basin information file) are read from the network topology file. Several processes must be carried out after reading the network topology file to conform to the requirements of MESH – such as reindexing, renaming, editing, and adding missing variables and coordinates. Second, the land cover fractions for each subbasin are calculated and appended to the working drainage database file. Creation of the drainage database file requires two inputs: the network topology file created above, and a land cover zonal histogram.

Setting	Description
input_lc_zh_path	refers to the location of the zonal statistics file produced with GIS Tool or QGIS.
input_lc_zh_name	refers to the name of the zonal statistics file produced with GIS Tool or QGIS. Using a QGIS-Generated zonal Histogram If you intend to create the zonal histogram with QGIS (following the steps in section 4.3.2), specify the name of the shapefile (including the extension) in this setting. For example: input_lc_zh_name \| BowAtBanff_lczh.shp
input_topo_path	refers to the location of the network topology file
input_topo_name	refers to the name of the network topology file
DDB_output_dir	refers to the output directory for the drainage database. By default, the output file is saved in /root_path/vector_based_workflow/workflow_data/[domain_name]/drainagedatabase/
lc_type_prefix	refers to the prefix given to each column representing a unique landcover type in the zonal statistics file. If the zonal statistics file was produced by GIS Tool, the output column prefix by default is ‘frac_’. If users utilize QGIS in calculating the zonal histogram, the output column prefix 'NALCMS_' or QGIS default value of 'HISTO_' can be used.
merit_basin_path & merit_basin_name	refer to the location of the catchment shapefile used to produce the network topology file. These settings are redundant and will be removed in a later version of the workflow. They are used in case the user wishes to complete a sanity check of the subbasin selection process (which must be enabled by altering the script – not recommended.)
land_cover_types	The 'land_cover_types' setting defines the names for land cover types (GRUs) in the drainage database. The default setting uses the 19 NALCMS land cover classes. If you have created custom GRUs or intend to use a different land cover product, you can specify your list of classes in 'land_cover_types'. The list of land cover types must be provided on a single line and each class name must be separated by a comma and a single space. The order of classes in 'land_cover_types' must correspond to the order of classes in the specified zonal histogram. EXAMPLE: You have a zonal histogram with columns "HISTO_1", "HISTO_2", and "HISTO_4" which correspond to "Forest-Pine", "Forest-Fir", and "Water" respectively. The 'land_cover_types' setting should read: Forest-Pine, Forest-Fir, Cropland, Water. "Forest-Pine" will represent HISTO_1, "Forest-Fir" will represent HISTO_2, "Cropland" will be exluded from the drainage database, and "Water" will represent HISTO_4. If 'land_cover_types' were to read: Forest, Grassland, Water, Cropland then "Water" would mistakenly be excluded and "Cropland" would be assigned to represent HISTO_4.

3.5 Climate Forcing Settings

Currently, the climate forcing settings only apply to the easymore_remapping.py script. The 'change unity', 'datatool', and 'reordering' scripts must have their input and output locations modified manually.

# Climate forcing settings
remapping_out			   | default
source_nc_path             | ../forcing/
source_nc_name             | rdrsv2.1_1980.nc
var_lon                    | lon
var_lat                    | lat
var_time                   | time
format_list                | f4
fill_value_list            | -9999.00
var_names                  | RDRS_v2.1_A_PR0_SFC, RDRS_v2.1_P_FI_SFC, RDRS_v2.1_P_FB_SFC, RDRS_v2.1_P_TT_09944, RDRS_v2.1_P_UVC_09944, RDRS_v2.1_P_P0_SFC, RDRS_v2.1_P_HU_09944 # variables must be on 1 line and separated by a comma and a space.

Setting	Description
remapping_out	The output location for the remapped forcing file. By default, uses root_path/forcing/domainname
source_nc_path	The path to the netcdf forcing file from which the output file will be remapped. This should be the location of the file obtained with the datatool.sh script.
source_nc_name	The name of the netcdf forcing file from which the output file will be remapped. This should be the name of the file obtained with the datatool script.
var_lon	The name of the longitude variable in the source netcdf file. (Usually leave this as 'lon')
var_lat	The name of the latitude variable in the source netcdf file. (Usually leave this as 'lat')
var_time	The name of the time variable in the source netcdf file. (Usually leave this as 'time'.)
format list	Leave this as 'f4'
fill_value_list	Leave this as '-9999.00'
var_names	The names of the forcing variables from the source file to be remapped. Typically leave these as they are. This setting has strict formatting. Place all variables on a single line and separate each variable with a comma and a single space.

4. Generating the Drainage Database (Basin Information File)

If running on a workstation machine or without submitting to the scheduler on Graham, users should first ensure that a virtual environment with the required Python packages is activated. For information on how to achieve this, consult Installing Python and the Windows Subsystem for Linux.

Scheduling Jobs on Graham

Graham users may wish to create individual shell scripts to submit each of the workflow scripts to the scheduler. If so, ensure that a virtual environment with the requisite packages is created and activated and that the necessary modules are loaded. Graham users are strongly encouraged to consult: running python jobs on HPC cluster and the Digital Alliance of Canada documentation on Graham

4.1. Subbasin Selection

Caution: PFAF 74 (Mississippi Basin)

If downloading the MERIT shapefiles directly from the link in section 2.2, the file riv_pfaf_74_MERIT_Hydro_v07_Basins_v01_bugfix1.shp contains a point feature that will cause the subbasin selection script to fail.

To fix this error, open the shapefile in GIS, find the row in the attribute table where 'COMID' equals 74030207 and change the 'NextDownID' for that row from 74030207 to 0.

To run the script on a Linux-alike system, navigate to the directory 0_subbasin_selection/ and run the following command:

python3 basin_subset.py

Shapefiles representing the subset domain of interest will be output to the folder specified in the control file (/root_path/vector_based_workflow/workflow_data/[domain_name]/shapefiles/catchment or shapefiles/river_network by default).

Example subbasin selection script:

#%% load modules
from easymore.easymore import easymore
import matplotlib.pyplot as plt
import geopandas as gpd
import numpy as np
from pathlib import Path
from shutil import copyfile
from datetime import datetime
import networkx as nx
from typing import (
    Union,
)

######################################### Control file handling #########################################
# Easy access to control file folder
print("Reading control file...")
controlFolder = Path('../0_control_files')
 
# Store the name of the 'active' file in a variable
controlFile = 'control_active.txt'
 
#Function to extract a given setting from the control file
def read_from_control( file, setting ):
     
    # Open 'control_active.txt' and ...
    with open(file) as contents:
        for line in contents:
             
            # ... find the line with the requested setting
            if setting in line and not line.startswith('#'):
                break
     
    # Extract the setting's value
    substring = line.split('|',1)[1]      # Remove the setting's name (split into 2 based on '|', keep only 2nd part)
    substring = substring.split('#',1)[0] # Remove comments, does nothing if no '#' is found
    substring = substring.strip()         # Remove leading and trailing whitespace, tabs, newlines
        
    # Return this value   
    return substring
 
#Function to specify a default path
def make_default_path(suffix):
     
    # Get the root path
    rootPath = Path( read_from_control(controlFolder/controlFile,'root_path') )
     
    # Get the domain folder
    #domain_name = read_from_control(controlFolder/controlFile,'domain_name')
    #domainFolder = 'domain_' + domain_name
     
    # Specify the forcing path
    defaultPath = rootPath / suffix
     
    return defaultPath

# Get the domain folder
domain_name = read_from_control(controlFolder/controlFile,'domain_name')
domainFolder = 'domain_' + domain_name

#%% assigning input files
# Note : this is an example for the Fraser setup
input_basin_path = read_from_control(controlFolder/controlFile,'input_basin_path')
input_basin_name = read_from_control(controlFolder/controlFile,'input_basin_name')
#input_basin = '../cat_pfaf_78_MERIT_Hydro_v07_Basins_v01_bugfix1.shp'
input_river_path = read_from_control(controlFolder/controlFile,'input_river_path')
input_river_name = read_from_control(controlFolder/controlFile,'input_river_name')
#input_river = '../rivEndoMERITpfaf_78.shp'
 
outdir_basin = read_from_control(controlFolder/controlFile,'subset_basin_outdir')
outdir_river = read_from_control(controlFolder/controlFile,'subset_river_outdir')

# Specify default path if needed
if input_basin_path == 'default':
    input_basin_path = make_default_path('shapefiles/catchment/') # outputs a Path()
else:
    input_basin_path = Path(input_basin_path) # make sure a user-specified path is a Path()
    
if input_river_path == 'default':
    input_river_path = make_default_path('shapefiles/river_network/') # outputs a Path()
else:
    input_river_path = Path(input_river_path) # make sure a user-specified path is a Path()
    
if outdir_basin == 'default':
    outdir_basin = make_default_path('shapefiles/catchment/') # outputs a Path()
else:
    outdir_basin = Path(outdir_basin) # make sure a user-specified path is a Path()
 
if outdir_river == 'default':
    outdir_river = make_default_path('shapefiles/river_network/') # outputs a Path()
else:
    outdir_river = Path(outdir_river) # make sure a user-specified path is a Path()

# Get additional settings
hruid = read_from_control(controlFolder/controlFile,'river_basin_shp_rm_hruid')
segid = read_from_control(controlFolder/controlFile,'river_network_shp_segid')
downid = read_from_control(controlFolder/controlFile,'river_network_shp_downsegid')
target_segment = int(read_from_control(controlFolder/controlFile,'target_segment')) # e.g. 78011863
case_name = read_from_control(controlFolder/controlFile,'case_name') #e.g. 05NG001


# Kasra's Fix
def find_upstream(
    gdf: gpd.GeoDataFrame,
    target_id: Union[str, int],
    main_id: str,
    ds_main_id: str,
):
    '''Find "ancestors" or upstream segments in a river network given
    in the from of a geopandas.GeoDataFrame `gdf`
    
    Parameters:
    -----------
    gdf: geopandas.GeoDataFrame
        GeoDataFrame of river segments including at least three pieces
        of information: 1) geometries of segments, 2) segment IDs, and
        3) downstream segment IDs
    target_id: str, int, or any other data type as included in `gdf`
        Indicating the target ID anscestor or upstream of which is
        desired
    main_id: str
        String defining the column of element IDs in the input geopandas
        dataframe
    ds_main_id: str
        String defining the column of downstream element IDs in the
        input geopandas dataframe
    
    Returns:
    --------
    nodes: set
        IDs of nodes being upstream or anscestor of the `target_id`
    
    '''
    # creating a DiGraph out of `gdf` object
    riv_graph = nx.from_pandas_edgelist(gdf, source=main_id, target=ds_main_id, create_using=nx.DiGraph)
    
    # return nodes in a list
    nodes = nx.ancestors(riv_graph, target_id)
    nodes.add(target_id)

    return nodes

print("loading source shapefiles...")
# load the files and calculating the downstream of each segment
riv  = gpd.read_file(input_river_path/input_river_name)
cat  = gpd.read_file(input_basin_path/input_basin_name)

print("calculating upstream river segments...")
upstream_segments   = find_upstream(riv, target_segment, segid, downid)

riv_up = riv.loc[riv[segid].isin(upstream_segments)]
riv_up = riv_up.set_crs(4326)

print("getting upstream catchments...")
cat_up = cat.loc[cat[hruid].isin(riv_up[hruid])]
cat_up = cat_up.set_crs(4326)

# plot
cat_up.plot()
riv_up.plot()
# save
cat_up.to_file(str(outdir_basin)+'/'+case_name+'_cat.shp')
riv_up.to_file(str(outdir_river)+'/'+case_name+'_riv.shp')
print("Done! outputs saved to the specified folder")
# --- Code provenance
# Generates a basic log file in the domain folder and copies the control file and itself there.
 
# Set the log path and file name
logPath = outdir_basin
log_suffix = '_basin_subset.txt'
 
# Create a log folder
logFolder = '_workflow_log'
Path( logPath / logFolder ).mkdir(parents=True, exist_ok=True)
 
# Copy this script
thisFile = 'basin_subset.py'
copyfile(thisFile, logPath / logFolder / thisFile);
 
# Get current date and time
now = datetime.now()
 
# Create a log file
logFile = now.strftime('%Y%m%d') + log_suffix
with open( logPath / logFolder / logFile, 'w') as file:
     
    lines = ['Log generated by ' + thisFile + ' on ' + now.strftime('%Y/%m/%d %H:%M:%S') + '\n',
             'Generated basin subset .shp files from {} and {}.'.format(input_basin_name,input_river_name)]
    for txt in lines:
        file.write(txt)

4.2. Network Topology

Navigate to the directory 1_network_topology/ and run the following command:

python3 create_network_topology_file.py

Example network topology script:

# modules
import os
import pandas as pd
import netCDF4 as nc4
import geopandas as gpd
import numpy as np
import xarray as xs
from pathlib import Path
from shutil import copyfile
from datetime import datetime
import warnings
  
#Control file handling
# Easy access to control file folder
controlFolder = Path('../0_control_files')
  
# Store the name of the 'active' file in a variable
controlFile = 'control_active.txt'
  
#Function to extract a given setting from the control file
def read_from_control( file, setting ):
      
    # Open 'control_active.txt' and ...
    with open(file) as contents:
        for line in contents:
              
            # ... find the line with the requested setting
            if setting in line and not line.startswith('#'):
                break
      
    # Extract the setting's value
    substring = line.split('|',1)[1]      # Remove the setting's name (split into 2 based on '|', keep only 2nd part)
    substring = substring.split('#',1)[0] # Remove comments, does nothing if no '#' is found
    substring = substring.strip()         # Remove leading and trailing whitespace, tabs, newlines
         
    # Return this value  
    return substring
  
#Function to specify a default path
def make_default_path(suffix):
      
    # Get the root path
    rootPath = Path( read_from_control(controlFolder/controlFile,'root_path') )
      
    # Get the domain folder
    #domainName = read_from_control(controlFolder/controlFile,'domain_name')
    #domainFolder = 'domain_' + domainName
      
    # Specify the forcing path
    #defaultPath = rootPath / domainFolder / suffix
    defaultPath = rootPath / suffix
    return defaultPath
 
# Get the domain folder
domain_name = read_from_control(controlFolder/controlFile,'domain_name')
domainFolder = 'domain_' + domain_name
 
# Find location of river network shapefile
# River network shapefile path & name
river_network_path = read_from_control(controlFolder/controlFile,'river_network_shp_path')
river_network_name = read_from_control(controlFolder/controlFile,'river_network_shp_name')
  
# Specify default path if needed
if river_network_path == 'default':
    river_network_path = make_default_path('shape_file/river_network/') # outputs a Path()
else:
    river_network_path = Path(river_network_path) # make sure a user-specified path is a Path()
      
#Find the field names we're after
river_seg_id      = read_from_control(controlFolder/controlFile,'river_network_shp_segid')
river_down_seg_id = read_from_control(controlFolder/controlFile,'river_network_shp_downsegid')
river_slope       = read_from_control(controlFolder/controlFile,'river_network_shp_slope')
river_length      = read_from_control(controlFolder/controlFile,'river_network_shp_length')
 
# Added by MESH workflow
try:
    river_outlet_id   = float(read_from_control(controlFolder/controlFile,'river_network_shp_outlet_id'))
except ValueError:
    print ('The ID of the most downstream segment was not found')
    river_outlet_id = []
  
#Find location of river basin shapefile (routing catchments)
# River network shapefile path & name
river_basin_path = read_from_control(controlFolder/controlFile,'river_basin_shp_path')
river_basin_name = read_from_control(controlFolder/controlFile,'river_basin_shp_name')
  
#Specify default path if needed
if river_basin_path == 'default':
    river_basin_path = make_default_path('shape_file/catchment/') # outputs a Path()
else:
    river_basin_path = Path(river_basin_path) # make sure a user-specified path is a Path()
  
#Find the field names we're after
basin_hru_id     = read_from_control(controlFolder/controlFile,'river_basin_shp_rm_hruid')
basin_hru_area   = read_from_control(controlFolder/controlFile,'river_basin_shp_area')
basin_hru_to_seg = read_from_control(controlFolder/controlFile,'river_basin_shp_hru_to_seg')
  
#Find where the topology file needs to go
# Topology .nc path and name
topology_path = read_from_control(controlFolder/controlFile,'settings_routing_path')
topology_name = read_from_control(controlFolder/controlFile,'settings_routing_topology')
  
#Specify default path if needed
if topology_path == 'default':
    topology_path = make_default_path('vector_based_workflow/workflow_data/domain_'+domain_name+'/topology/') # outputs a Path()
else:
    topology_path = Path(topology_path) # make sure a user-specified path is a Path()
  
#Make the folder if it doesn't exist
topology_path.mkdir(parents=True, exist_ok=True)
  
# Make the river network topology file
# Open the shapefile
shp_river = gpd.read_file(river_network_path/river_network_name)
shp_basin = gpd.read_file(river_basin_path/river_basin_name)
  
# Added by MESH workflow
# sort basin to be consistent with river
shp_basin = shp_basin.sort_values(by=basin_hru_id)
  
# convert area to m^2
# Note: if area unit is already based on m**2, it is not requried to covert m**2
shp_basin[basin_hru_area].values[:] = shp_basin[basin_hru_area].values[:]*10**6
  
# covert river_length to m
# Note: if length unit is already based on m, it is not requried to covert m
shp_river[river_length].values[:]   = shp_river[river_length].values[:]*1000
  
# adding centroid of each subbasin.
# Note: the more accurate should be done in equal area projection
warnings.simplefilter('ignore') # silent the warning
shp_basin['lon'] = shp_basin.centroid.x
shp_basin['lat'] = shp_basin.centroid.y
warnings.simplefilter('default') # back to normal
  
# specifying other variables
# Note: the river width and manning is optional. The manning coefficient is specified in the MESH
# hydrology configuration file
shp_river['width']   = 50
shp_river['manning'] = 0.03
  
# finished edit by MESH workflow
  
# Find the number of segments and subbasins
num_seg = len(shp_river)
num_hru = len(shp_basin)
  
# Ensure that the most downstream segment in the river network has a downstream_ID of 0
# This indicates to routing that this segment has no downstream segment attached to it
# Modified by MESH workflow
if (np.size(river_outlet_id)!= 0):
    shp_river.loc[shp_river[river_seg_id] == river_outlet_id, river_down_seg_id] = 0
  
# Function to create new nc variables
def create_and_fill_nc_var(ncid, var_name, var_type, dim, fill_val, fill_data, long_name, units):
      
    # Make the variable
    ncvar = ncid.createVariable(var_name, var_type, (dim,), fill_val)
      
    # Add the data
    ncvar[:] = fill_data  
      
    # Add meta data
    ncvar.long_name = long_name
    ncvar.unit = units
      
    return
  
# Make the netcdf file
with nc4.Dataset(topology_path/topology_name, 'w', format='NETCDF4') as ncid:
      
    # Set general attributes
    now = datetime.now()
    ncid.setncattr('Author', "Created by MESH vector-based workflow scripts")
    ncid.setncattr('History','Created ' + now.strftime('%Y/%m/%d %H:%M:%S'))
    ncid.setncattr('Purpose','Create a river network .nc file for WATROUTE routing')
      
    # Define the seg and hru dimensions
    # it can be renamed to 'subbasin'
    # Modified by MESH workflow
    ncid.createDimension('n', num_seg)
    # ncid.createDimension('hru', num_hru)
    # finished edit by MESH workflow
  
    # --- Variables
    # renaming variable and adding lat, lon, manning and width
    # create_and_fill_nc_var(ncid, 'seg_id', 'int', 'n', False, \
                           # shp_river[river_seg_id].values.astype(int), \
                           # 'Unique ID of each stream segment', '-')
    # Modified by MESH workflow                      
    create_and_fill_nc_var(ncid, 'seg_id', 'int', 'n', False, \
                           np.ndarray.round(shp_river[river_seg_id].values).astype(int), \
                           'Unique ID of each stream segment', '-')                      
                            
    # create_and_fill_nc_var(ncid, 'tosegment', 'int', 'n', False, \
                           # shp_river[river_down_seg_id].values.astype(int), \
                           # 'ID of the downstream segment', '-')
     
    # Modified by MESH workflow
    create_and_fill_nc_var(ncid, 'tosegment', 'int', 'n', False, \
                           np.ndarray.round(shp_river[river_down_seg_id].values).astype(int), \
                           'ID of the downstream segment', '-')
     
    create_and_fill_nc_var(ncid, 'slope', 'f8', 'n', False, \
                           shp_river[river_slope].values.astype(float), \
                           'Segment slope', '-')  
    # added by MESH workflow
    create_and_fill_nc_var(ncid, 'lon', 'f8', 'n', False, \
                           shp_basin['lon'].values.astype(float), \
                           'longitude', '-')    
    create_and_fill_nc_var(ncid, 'lat', 'f8', 'n', False, \
                           shp_basin['lat'].values.astype(float), \
                           'latitude', '-')
    # finished edit by MESH workflow  
    create_and_fill_nc_var(ncid, 'length', 'f8', 'n', False, \
                           shp_river[river_length].values.astype(float), \
                           'Segment length', 'm')
    create_and_fill_nc_var(ncid, 'hruid', 'int', 'n', False, \
                           shp_basin[basin_hru_id].values.astype(int), \
                           'Unique hru ID', '-')
    create_and_fill_nc_var(ncid, 'seg_hr_id', 'int', 'n', False, \
                           shp_basin[basin_hru_to_seg].values.astype(int), \
                           'ID of the stream segment to which the HRU discharges', '-')
    create_and_fill_nc_var(ncid, 'basin_area', 'f8', 'n', False, \
                           shp_basin[basin_hru_area].values.astype(float), \
                           'HRU area', 'm^2')  
    # added by MESH workflow
    create_and_fill_nc_var(ncid, 'width', 'f8', 'n', False, \
                           shp_river['width'].values.astype(float), \
                           'width', 'm')                     
    create_and_fill_nc_var(ncid, 'manning', 'f8', 'n', False, \
                           shp_river['manning'].values.astype(float), \
                           'manning', '-') 
   # finished edit by MESH workflow
  
  
# --- Code provenance
# Generates a basic log file in the domain folder and copies the control file and itself there.
  
# Set the log path and file name
logPath = topology_path
log_suffix = '_make_river_network_topology.txt'
  
# Create a log folder
logFolder = '_workflow_log'
Path( logPath / logFolder ).mkdir(parents=True, exist_ok=True)
  
# Copy this script
thisFile = 'create_network_topology_file.py'
copyfile(thisFile, logPath / logFolder / thisFile);
  
# Get current date and time
now = datetime.now()
  
# Create a log file
logFile = now.strftime('%Y%m%d') + log_suffix
with open( logPath / logFolder / logFile, 'w') as file:
      
    lines = ['Log generated by ' + thisFile + ' on ' + now.strftime('%Y/%m/%d %H:%M:%S') + '\n',
             'Generated network topology .nc file.']
    for txt in lines:
        file.write(txt)

4.3. Zonal Statistics

There are currently two methods to produce the necessary zonal statistics representing land cover types spanning a domain of interest. The preferred method makes use of the GIS Tool package because it eliminates the need to process spatial data in a desktop GIS environment. However, the GIS Tool Method requires access to the Digital Alliance of Canada cluster Graham and the Landsat data stored there. Users without Graham access can retrieve the NALCMS land cover data from the CEC website, and make use of the alternative QGIS method. Regardless of the method used, two inputs are required: the land cover TIFF file and the MERIT Hydro catchment shapefile.

Though control file support for GIS Tool in the vector-based workflow is planned for a future version, users must currently enter arguments to the landcover_zonalhist.sh script to extract the requisite zonal statistics. A general description of GIS Tool is available here.

4.3.1. Running GIS Tool for the MESH vector-based workflow

GIS Tool will only work on Graham. Users without access to Graham should consult section 4.3.2 to calculate zonal statistics with QGIS

The script landcover_zonalhist.sh has been created to allow users to conveniently specify arguments to extract-gis.sh, which is responsible for running the GIS Tool processes. The table below describes the necessary arguments, their purpose, and the recommended settings.

landcover_zonalhist.sh:

#!/bin/bash
 
# Geospatial Dataset Processing Workflow
# Copyright (C) 2022, University of Saskatchewan
#
# This file is part of the Geospatial Dataset Processing Workflow
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
# This is a simple example to extract common statistics for the
# Bow At Banff Shapefiles of Landsat landcover
# from the MERIT-Hydro rasters.
 
# Always call the script in the root directory of this repository
# cd ..
# echo "The current directory is: $(pwd)"
 
# # first download a sample shapefile - could be any shapefile
# wget -m -nd -nv -q -A "cat_pfaf_71_MERIT_Hydro_v07_Basins_v01_bugfix1.*" \
     # "http://hydrology.princeton.edu/data/mpan/MERIT_Basins/MERIT_Hydro_v07_Basins_v01_bugfix1/pfaf_level_02/";
 
 
# implement subsetting and zonal statistics
./extract-gis.sh --dataset="landsat" \
  --dataset-dir="/project/rpp-kshook/Model_Output/Landsat/" \
  --variable="NA_NALCMS_2010_v2_land_cover_30m" \
  --shape-file="$(pwd)/../../shape_file/catchment/BowAtBanff_cat.shp" \
  --print-geotiff=true \
  --output-dir="$(pwd)/../workflow_data/domain_BowAtBanff/zonalhist/" \
  --prefix="landsat_BowAtBanff_" \
  --stat="majority,minority,frac" \
  --email=your.email@something.ca \
  -j;

Ensure that execute permission is enabled for the following files before executing landsat_zonalhist.sh:

landsat_zonalhist.sh
extract-gis.sh
landsat.sh

To ensure that execute permission is enabled for all files inside the gistool folder, navigate to the 2_zonal_statisitcs/ directory and use the command:

chmod -R +x 2a_gistool

Execute landcover_zonalhist.sh to begin extracting the zonal statistics.

./landcover_zonalhist.sh

By default, the job will be sent to the scheduler. To disable job submission, remove the ‘-j’ argument from landcover_zonalhist.sh

By default, the output will be placed in root_path/vector_based_workflow/workflow_data/[domain_name]/zonalhist/

4.3.2. QGIS zonal histogram

The following steps show how to calculate the zonal histogram using the Python terminal in QGIS.

4.3.2.1. Fixing invalid geometry

The polygons of the MERIT Hydro shapefiles may exhibit a ring self-intersection problem and this will produce an 'invalid geometry' error message when running the zonal histogram. Therefore it is required that this problem is fixed before running the calculation.

To solve the invalid geometry error, QGIS has a built-in algorithm to fix this error automatically. In the Processing Toolbox, search for the fix geometry, then double click on it to run.

After fixing the geometry error, the zonal histogram can be obtained.

4.3.2.2. Obtaining the zonal histogram using the Python terminal in QGIS

The following Python code can be run inside QGIS to obtain the zonal histogram.

lc_extract.py:

lc_extract.py

"""
@ Author : MESHworkflow
 
"""
 
# %% import module
import processing
 
# %% run the algorithm
input_raster = '../../landcover.tif'
input_shape  = '../../MERITHydro_cat_fixgeo.shp'
output_shpae = '../../domain_name_intersect_NGRU.shp'
band           =  1
 
processing.run("native:zonalhistogram" , {'COLUMN_PREFIX' : 'NALCMS_',
            'INPUT_RASTER' : input_raster,
            'INPUT_VECTOR' : input_shape,
            'OUTPUT'       : output_shpae,
            'RASTER_BAND'  : band })

4.4. Calculate Landcover Fraction and Generate the Drainage Database (Basin Information File)

Navigate to the directory 3_drainage database/ and run the following command:

python3 create_MESH_drainage_database.py

Example drainage database script:

create_MESH_drainage_database.py

# -*- coding: utf-8 -*-
"""
Name
    create_MESH_drainage_database. previously called lc_vectorbased
Purpose
    The purpose of this script is to calculate land cover fractions for each
    subbasin of interest. Then the landcover is converted the (subbasin*lc_types)
    it is adhered to the driange database. 
Programmer(s)
    Ala Bahrami
    Cooper Albano
Revision History
    20210422 -- Initial version created 
    20210504 -- 1) changed dimension name from 'lc_type' to 'ngru' and variable
    name from 'lc_frac' to 'GRU'. 2) Added LandUse variable
    20210505 -- append the LandUse information to drainage_ddb
    20210506 -- changed 'ngru' dimension to 'gru' to be consistent with MESH code
    20210604 -- modified I/O and variables for Fraser application
    20210704 -- modified to the new_rank_extract function
    20220623 -- 1) Modified based on the new_rank_modi2 which is adaptable for 
                multi-outlet
                2) Consider the entire 19 land cover classes instead of regrouping
                3) save subbasin reordered metadata
                4) visualize and save subbasin selection for any outlet (optional)
    20220626 -- 1) modified the way of reindexing the zonal histogram
    20220825 -- 1) modify code to adapt NEXT variable having multiple outlets
    20221101 -- 1) added a line to accept GIS tool .csv zonal hist as pandas dataframe
             -- 2) added a line to rename prefix_0 to prefix_NOD for the GIS tool .csv zonal hist
             -- 3) removed landcover fraction calculation for GIS tool .csv file
    20221116 -- 1) removed .csv column reordering. No longer necessary due to gistool bug fix
    20221123 -- 1) added functionality to read I/O from control file.
See also
     
Reference
 
Todo:
    1) The lc_types is based on NALCMS 2010. The name list is hard-coded  
"""
 
# %% importing modules
import geopandas as gpd
import numpy as np
import xarray as xs
import pandas as pd
from   datetime import date
from pathlib import Path
from shutil import copyfile
import time
 
#Control file handling
# Easy access to control file folder
controlFolder = Path('../0_control_files')
 
# Store the name of the 'active' file in a variable
controlFile = 'control_active.txt'
 
#Function to extract a given setting from the control file
def read_from_control( file, setting ):
      
    # Open 'control_active.txt' and ...
    with open(file) as contents:
        for line in contents:
              
            # ... find the line with the requested setting
            if setting in line and not line.startswith('#'):
                break
      
    # Extract the setting's value
    substring = line.split('|',1)[1]      # Remove the setting's name (split into 2 based on '|', keep only 2nd part)
    substring = substring.split('#',1)[0] # Remove comments, does nothing if no '#' is found
    substring = substring.strip()         # Remove leading and trailing whitespace, tabs, newlines
         
    # Return this value  
    return substring
 
#Function to specify a default path
def make_default_path(suffix):
      
    # Get the root path
    rootPath = Path( read_from_control(controlFolder/controlFile,'root_path') )
      
    # Get the domain folder
    #domainName = read_from_control(controlFolder/controlFile,'domain_name')
    #domainFolder = 'domain_' + domainName
      
    # Specify the forcing path
    defaultPath = rootPath / suffix
      
    return defaultPath
 
# Get the domain folder
domain_name = read_from_control(controlFolder/controlFile,'domain_name')
domainFolder = 'domain_' + domain_name
 
## Find location of zonal statistics file
# Zonal statistics file path & name
lc_zh_path = read_from_control(controlFolder/controlFile,'input_lc_zh_path')
lc_zh_name = read_from_control(controlFolder/controlFile,'input_lc_zh_name')
 
# Specify default path if needed
if lc_zh_path == 'default':
    lc_zh_path = make_default_path('vector_based_workflow/workflow_data/domain_'+domain_name+'/zonalhist/') # outputs a Path()
else:
    lc_zh_path = Path(lc_zh_path) # make sure a user-specified path is a Path()
 
## Find location of network topology file
# Network topology file path & name
topo_path = read_from_control(controlFolder/controlFile,'input_topo_path')
topo_name = read_from_control(controlFolder/controlFile,'input_topo_name')
 
# Specify default path if needed
if topo_path == 'default':
    topo_path = make_default_path('vector_based_workflow/workflow_data/domain_'+domain_name+'/topology/') # outputs a Path()
else:
    topo_path = Path(topo_path) # make sure a user-specified path is a Path()
 
## Find location of MERIT Basin shapefile
# MERIT Bassin shapefile path & name
merit_path = read_from_control(controlFolder/controlFile,'merit_basin_path')
merit_name = read_from_control(controlFolder/controlFile,'merit_basin_name')
 
# Specify default path if needed
if merit_path == 'default':
    merit_path = make_default_path('shape_file/catchment/') # outputs a Path()
else:
    merit_path = Path(merit_path) # make sure a user-specified path is a Path()
 
## Find location of output directory
# output directory path
outdir = read_from_control(controlFolder/controlFile,'DDB_output_dir')
 
# Specify default path if needed
if outdir == 'default':
    outdir = make_default_path('vector_based_workflow/workflow_data/domain_'+domain_name+'/drainagedatabase/') # outputs a Path()
else:
    outdir = Path(outdir) # make sure a user-specified path is a Path()
 
# %% directory of input files
# Enter path to a zonal histogram file in either .csv format from GIS tool or in .shp format from QGIS
start_time = time.time()
input_lc_zh              = lc_zh_path/lc_zh_name    
input_topology           = topo_path/topo_name
Merit_catchment_shape    = merit_path/merit_name
domain_name              = read_from_control(controlFolder/controlFile,'domain_name')
lc_type_prefix           = read_from_control(controlFolder/controlFile,'lc_type_prefix')
 
#%% Function reindex to extract drainage database variables
def new_rank_extract(input_topology):
        #% Reading topology file and finding outlets
        drainage_db = xs.open_dataset(input_topology)
        drainage_db.close()
 
        segid = drainage_db['seg_id'].values
        tosegment = drainage_db['tosegment'].values
 
        # Count the number of outlets
        outlets = np.where(tosegment == 0)[0]
 
        #% Search over to extract the subbasins drain into each outlet
        rank_id_domain = np.array([]).astype(int)  
        outlet_number = np.array([]).astype(int)
        for k in range(len(outlets)):
            # initial step
            #segid_target = drainage_db['seg_id'].values[outlets[k]]
            segid_target = segid[outlets[k]]
            # set the rank of the outlet
            rank_id = outlets[k]
             
            # find upstream segids drains into downstream
            while(np.size(segid_target) >= 1):
                if (np.size(segid_target) == 1):
                    r = np.where(tosegment == segid_target)[0]
                else:
                    r = np.where(tosegment == segid_target[0])[0]   
                # updated the target segid
                segid_target = np.append(segid_target, segid[r])
                # remove the first searched target
                segid_target = np.delete(segid_target,0,0)
                if (len(segid_target) == 0):
                    break
                # update the rank_id
                rank_id = np.append(rank_id,r)
            rank_id = np.flip(rank_id)
            if (np.size(rank_id) > 1):
                outlet_number = np.append(outlet_number, (k)*np.ones((len(rank_id),1)).astype(int))
            else:
                outlet_number = np.append(outlet_number, (k))
            rank_id_domain = np.append(rank_id_domain, rank_id)
            rank_id = []
        #% reorder segid and tosegment
        segid = segid[rank_id_domain]
        tosegment = tosegment[rank_id_domain]        
               
        # rearrange outlets to be consistent with MESH outlet structure
        # NB: In MESH outlets should be placed at the end of NEXT variable
        NA = len(rank_id_domain)
        fid1 = np.where(tosegment != 0)[0]
        fid2 = np.where(tosegment == 0)[0]
        fid =  np.append(fid1,fid2)
         
        rank_id_domain = rank_id_domain[fid]
        segid =segid[fid]
        tosegment = tosegment[fid]
        outlet_number = outlet_number[fid]
         
        #% construct Rank and Next variables
        Next = np.zeros(NA).astype(np.int32)
         
        for k in range(NA):
            if (tosegment[k] != 0):
                r = np.where(tosegment[k] == segid)[0] + 1
                Next[k] = r
            else:
                Next[k] = 0
                 
        # Construct Rank from 1:NA
        Rank = np.arange(1,NA+1).astype(np.int32)
         
        #% save subbasins reordered metadata
        dt = {'Merit_reorderd_ID':rank_id_domain, 'Outlet_Number':outlet_number,
              'Rank':Rank,'Next':Next,'Segid':segid,'tosegment':tosegment}
        df = pd.DataFrame(data=dt, dtype = np.int64)
        outrank = domain_name+'_Rank_ID'+'.csv'
        df.to_csv(outdir/outrank, index=False)
         
        # % reordering network topology variables based on Rank 1:NA
        for m in ['basin_area', 'length', 'slope', 'lon', 'lat', 'hruid',
                  'seg_id', 'seg_hr_id', 'tosegment', 'width', 'manning']:
            drainage_db[m].values = drainage_db[m].values[rank_id_domain]
             
        # % check if channel slope values match the minimum threshold
        min_slope = 0.000001
        drainage_db['slope'].values[drainage_db['slope'].values < min_slope] = min_slope
         
        # % Adding Rank and Next variables to the file
        drainage_db['Rank'] = (['n'], Rank)
        drainage_db['Rank'].attrs.update(standard_name = 'Rank',
                            long_name = 'Element ID', units = '1', _FillValue = -1)
         
        drainage_db['Next'] = (['n'], Next)
        drainage_db['Next'].attrs.update(standard_name = 'Next',
                           long_name = 'Receiving ID', units = '1', _FillValue = -1)
 
        # % Adding missing attributes and renaming variables
        # Add 'axis' and missing attributes for the 'lat' variable.
        drainage_db['lat'].attrs['standard_name'] = 'latitude'
        drainage_db['lat'].attrs['units'] = 'degrees_north'
        drainage_db['lat'].attrs['axis'] = 'Y'
          
        # Add 'axis' and missing attributes for the 'lon' variable.
        drainage_db['lon'].attrs['standard_name'] = 'longitude'
        drainage_db['lon'].attrs['units'] = 'degrees_east'
        drainage_db['lon'].attrs['axis'] = 'X'
          
        # Add or overwrite 'grid_mapping' for each variable (except axes).
        for v in drainage_db.variables:
            if (drainage_db[v].attrs.get('axis') is None):
                drainage_db[v].attrs['grid_mapping'] = 'crs'
          
        # Add the 'crs' itself (if none found).
        if (drainage_db.variables.get('crs') is None):
            drainage_db['crs'] = ([], np.int32(1))
            drainage_db['crs'].attrs.update(grid_mapping_name = 'latitude_longitude', longitude_of_prime_meridian = 0.0, semi_major_axis = 6378137.0, inverse_flattening = 298.257223563)
          
        # Rename variables.
        for old, new in zip(['basin_area', 'length', 'slope', 'manning'], ['GridArea', 'ChnlLength', 'ChnlSlope', 'R2N']):
            drainage_db = drainage_db.rename({old: new})
          
        # Rename the 'subbasin' dimension (from 'n').
        drainage_db = drainage_db.rename({'n': 'subbasin'})
         
        # % Specifying the NetCDF "featureType"
        # Add a 'time' axis with static values set to today (in this case, time is not actually treated as a dimension).
        drainage_db['time'] = (['subbasin'], np.zeros(len(rank_id_domain)))
        drainage_db['time'].attrs.update(standard_name = 'time', units = ('days since %s 00:00:00' % date.today().strftime('%Y-%m-%d')), axis = 'T')
          
        # Set the 'coords' of the dataset to the new axes.
        drainage_db = drainage_db.set_coords(['time', 'lon', 'lat'])
          
        # Add (or overwrite) the 'featureType' to identify the 'point' dataset.
        drainage_db.attrs['featureType'] = 'point'
         
        return rank_id_domain, drainage_db, outlet_number
     
# %% calling the new_rank_extract
rank_id_domain, drainage_db, outlet_number = new_rank_extract(input_topology)
 
# %% reading the input zonal histogram of landcover and reindex it.
if str(input_lc_zh).endswith('.shp'):
    lc_zonal_hist = gpd.read_file(input_lc_zh)                        # read QGIS .shp zonal histogram
    lc_zonal_hist = lc_zonal_hist.sort_values(by=['COMID'])           # sort by COMID for QGIS zonal histogram
elif str(input_lc_zh).endswith('.csv'):
    lc_zonal_hist = pd.read_csv(input_lc_zh)                           # read GIS tool .csv zonal histogram
    lc_zonal_hist = lc_zonal_hist.sort_values(by=['p[[1]]'])           # sort by p[[1]] for GIS tool zonal histogram
else:
    print('Zonal histogram not recognized.')
    exit()
 
# rename frac_0 to frac_NOD for compatibility with verify lc_types. Not necessary for QGIS version.
if str(input_lc_zh).endswith('.csv'):
    lc_zonal_hist = lc_zonal_hist.rename(columns={lc_type_prefix+'0':lc_type_prefix+'NOD'})
    cols = lc_zonal_hist.columns.tolist()
    for i in cols:
        if lc_type_prefix in i:
            if 'NOD' in i:
                nod=i
                cols.remove(i)
                cols.append(nod)
    lc_zonal_hist = lc_zonal_hist[cols]
 
#%% reading source MeritHydro catchment file and visualize and save subbasin selection
## NB: this section can be uncommented if a user want to do a sanity check of the subbasin selection
## list of major segid_target outlet ids per each PFAF
## {78011862 (Fraser), 78017388(columbia), 82000048(MRB),
## 83012503, 71004266 (Hudson), 72039675 (St.Laurent),
## Mississipi (74072586), 73017442, 81018374 (Yukon), 77032206,
## 75022612, 75038087 (Hondo River), 75038096 (Usumacinta)}
 
# shape_catchment = gpd.read_file(Merit_catchment_shape)
# shape_catchment = shape_catchment.sort_values(by=['COMID'])
# shape_catchment.reset_index(drop=True, inplace=True)
 
# segid = drainage_db['seg_id'].values
# segid_target = 75038096
# r = np.where(segid == segid_target)[0]
# r2 = np.where(outlet_number == outlet_number[r])[0]
# rank_id = rank_id_domain[r2]
 
# shape_catchment.loc[rank_id].plot(color='white', edgecolor='black')
# shape_catchment.loc[rank_id].to_file(outdir+'PFAF_subselect_'+'%d'%segid_target+'.shp')
 
# %% land class types
# NB: the NOD here represent the No-data. The NALCMS data has no-data category which its values is zero
lc_type = np.array(['Temperate/sub-polar needleleaf forest','Sub-polar taiga needleleaf forest','Tropical/sub-tropical broadleaf evergreen forest','Tropical/sub-tropical broadleaf deciduous forest',
           'Temperate/sub-polar broadleaf deciduous forest','Mixed Forest','Tropical/sub-tropical shrubland', 'Temperate/sub-polar shrubland',
           'Tropical/sub-tropical grassland','Temperate/sub-polar grassland','Sub-polar/polar shrubland-lichen-moss','Sub-polar/polar grassland-lichen-moss',
           'Sub-polar/polar barren-lichen-moss','Wetland','Cropland','Barren Lands',
           'Urban-Built-up','Water','Snow-Ice','No-data'])
  
# %% verify list of lc types
m = len(lc_type) + 1
st = [];
p  = [];
for i in (range(1,m)):
    if (i < m-1):
        st1 =   lc_type_prefix+ str(i)
    else:
        st1 =   lc_type_prefix+ 'NOD'
     
    st = np.append(st, st1)
    fid = np.where(lc_zonal_hist.columns == st1)[0]
    if (fid.size == 0):
        print ('land cover %s is not presented in the list of NALCMS land cover for this PFAF' % lc_type[i-1])
        p = np.int32(np.append(p, i-1))
         
# add dummy land cover type required by MESH
lc_type = np.append(lc_type, 'Dump')   
 
# remove missing land cover types from  the list
if (len(p) != 0) :
    lc_type = np.delete(lc_type, p)
     
#%% calculate land cover fraction
# extract land cover zonal hist
lc_frac = lc_zonal_hist.filter(like=lc_type_prefix, axis = 1)
 
# NB: Based on NALCMS LANDSAT data, the open water data are classified as No-DATA.
# NB: So, if the catchments have some no-data, users should verify if it falls inside the open-water
# NB: that later be added to t he 'Water' land cover class.
 
# here the it is required to add NALCMS-NOD to land class type of Water if No-data is included in lc_type
fid = np.where(lc_type == 'No-data')[0]
if (fid.size != 0):
    r1 = np.where(lc_type == 'Water')[0]
    r2 = np.where(lc_type == 'No-data')[0] 
    print(r1,r2)
    # adding the nodata values to the water land cover type and drop it and remove from lc_type
    lc_frac.values[:,r1] = lc_frac.values[:,r1] + lc_frac.values[:,r2]
    lc_frac = lc_frac.drop(lc_frac.columns[r2], axis=1)
    lc_type = np.delete(lc_type, r2)
 
# add Dump layer for MESH application
lc_frac['Dump'] = 0
 
# calculating land cover percentage. Only calculate if input zonal histogram is a shapefile (i.e. QGIS version)
if str(input_lc_zh).endswith('.shp'):
    lc_frac = lc_frac.apply(lambda x: round(x/x.sum(),2), axis=1)
 
# %% convert the lc_frac as a dataset and save it as netcdf
lon = drainage_db['lon'].values
lat = drainage_db['lat'].values
tt = drainage_db['time'].values
 
lc_ds =  xs.Dataset(
    {
        "GRU": (["subbasin", "gru"], lc_frac.values),
        "LandUse": (["gru"], lc_type),
    },
    coords={
        "lon": (["subbasin"], lon),
        "lat": (["subbasin"], lat),
        "time": tt,
    },
)
 
# meta data attributes
lc_ds.attrs['Conventions'] = 'CF-1.6'
lc_ds.attrs['License']     = 'The data were written by Ala Bahrami'
lc_ds.attrs['history']     = 'Created on April 23, 2021'
lc_ds.attrs['featureType'] = 'point'         
 
# editing lat attribute
lc_ds['lat'].attrs['standard_name'] = 'latitude'
lc_ds['lat'].attrs['units'] = 'degrees_north'
lc_ds['lat'].attrs['axis'] = 'Y'
  
# editing lon attribute
lc_ds['lon'].attrs['standard_name'] = 'longitude'
lc_ds['lon'].attrs['units'] = 'degrees_east'
lc_ds['lon'].attrs['axis'] = 'X'
 
# editing time attribute
lc_ds['time'].attrs.update(standard_name = 'time',
                                 units = ('days since %s 00:00:00' % date.today().strftime('%Y-%m-%d')),
                                 axis = 'T')
 
# coordinate system
lc_ds['crs'] = drainage_db['crs'].copy()
 
## made changes here
#lc_ds['GRU'].values = lc_ds['GRU'].values[reorder,:]
 
# %% Append land cover information to existing drainage database
drainage_db["GRU"] = (["subbasin", "gru"], lc_frac.values)
drainage_db['GRU'].attrs['standard_name'] = 'GRU'
drainage_db['GRU'].attrs['long_name'] = 'Group Response Unit'
drainage_db['GRU'].attrs['units'] = '-'
drainage_db['GRU'].attrs['_FillValue'] = -1
 
drainage_db["LandUse"] = (["gru"], lc_type)
 
# Set the 'coords' of the dataset to the new axes.
drainage_db = drainage_db.set_coords(['time', 'lon', 'lat'])
 
# saved the drainage_database
outDDB = domain_name+'_MESH_drainage_database.nc'
drainage_db.to_netcdf(outdir/outDDB)
  
# %% Save land cover fraction (this is optional)
outFRAC = domain_name+'_MESH_LC_FRAC.nc'
lc_ds.to_netcdf(outdir/outFRAC)
print('--%s seconds--' %(time.time() - start_time))

The drainage database will be output to the folder specified in the control file (/root_path/vector_based_workflow/workflow_data/[domain_name]/drainagedatabase/ by default).

5. Climate forcing

5.1. Extracting climate forcing

For users with access to the Graham HPC cluster, the datatool.sh script can be used to extract the desired RDRSv2.1 variables into daily forcing files. Later, daily files can be combined into monthly or yearly files using forcing_merge.sh.

datatool.sh

#!/bin/bash
#SBATCH --account=rpp-kshook
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=1
#SBATCH --mem=256M
#SBATCH --time=00:10:00
#SBATCH --job-name=fextract1980
#SBATCH --error=errors80
#SBATCH --mail-user=example.email@usask.ca
#SBATCH --mail-type=BEGIN
#SBATCH --mail-type=END
 
cd
rm -rf $HOME/scratch/rdrs_output/
rm -rf datatool
git clone https://github.com/kasra-keshavarz/datatool # clone the repository
cd datatool
 
./extract-dataset.sh  --dataset=RDRS \
  --dataset-dir="/project/rpp-kshook/Model_Output/RDRSv2.1" \
  --output-dir="$HOME/scratch/rdrs_output/" \
  --start-date="1980-01-01 00:00:00" \
  --end-date="1980-12-31 00:00:00" \
  --lat-lims=0,89  \
  --lon-lims=-182,0 \
  --variable="RDRS_v2.1_P_P0_SFC,RDRS_v2.1_P_HU_09944,RDRS_v2.1_P_TT_09944,RDRS_v2.1_P_UVC_09944,RDRS_v2.1_A_PR0_SFC,RDRS_v2.1_P_FB_SFC,RDRS_v2.1_P_FI_SFC" \
  --prefix="rdrsv2.1_" \
  --email="example.email@usask.ca" \
  -j;
   
  cd ../scratch

forcing_merge.sh

#!/bin/bash
#SBATCH --account=rpp-kshook
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=1
#SBATCH --mem=60G
#SBATCH --time=01:00:00
#SBATCH --job-name=RDRS_forcing_merge
#SBATCH --error=errors_RDRS_forcing_merge
#SBATCH --mail-user=example.email@usask.ca
#SBATCH --mail-type=BEGIN
#SBATCH --mail-type=END
 
# load modules and combine forcing NetCDF
module load cdo/1.9.8
 
# cd to datatool.sh output directory
cd /scratch/$User/rdrs_output/1980


# merge into monthly files - Uncomment these lines to merge daily files into monthly files. 
#for i in 01 02 03 04 05 06 07 08 09 10 11 12
#do
#cdo -z zip -b F32 mergetime rdrsv2.1_1980$i*.nc rdrsv2.1_19080_*i.nc
#done
 
# merge into yearly files
cdo -z zip -b F32 mergetime rdrsv2.1_1980*.nc rdrsv2.1_1980.nc

5.2. Remapping climate forcing variables

The purpose of this section is to read gridded climate forcing variables and remap them as subbasin-average values to match the subbasin locations derived in the previous steps. In this application, the EASYMORE package is used for remapping. In the following code example, the RDRS_V2 climate forcing is used. Two input files are required for this application, the gridded climate forcing and the MERIT Hydro catchment shapefile that was used to derive the "network_topology" file.

easymore_remapping.py

#!/cvmfs/soft.computecanada.ca/easybuild/software/2020/avx2/Core/python/3.8.10/bin/python
# which python to find the location of installed pyhthon on graham

# %% loading EASYMORE
from easymore.easymore import easymore
from pathlib import Path
from shutil import copyfile
from datetime import datetime


# %%Control file handling
controlFolder = Path('../0_control_files') # Easy access to control file folder
controlFile = 'control_active.txt' # Store the name of the 'active' file in a variable


#Function to extract a given setting from the control file
def read_from_control( file, setting ):
# Open 'control_active.txt' and ...
with open(file) as contents:
for line in contents:
# ... find the line with the requested setting
if setting in line and not line.startswith('#'):
break 
# Extract the setting's value
substring = line.split('|',1)[1] # Remove the setting's name (split into 2 based on '|', keep only 2nd part)
substring = substring.split('#',1)[0] # Remove comments, does nothing if no '#' is found
substring = substring.strip() # Remove leading and trailing whitespace, tabs, newlines
# Return this value 
return substring


#Function to specify a default path
def make_default_path(suffix):
# Get the root path
rootPath = Path( read_from_control(controlFolder/controlFile,'root_path') )
defaultPath = rootPath / suffix 
return defaultPath


# %% Get the domain folder
domain_name = read_from_control(controlFolder/controlFile,'domain_name')
domainFolder = 'domain_' + domain_name 

# %% Get the target shapefile
target_shp_path = read_from_control(controlFolder/controlFile,'river_basin_shp_path')
# Specify default path if needed
if target_shp_path == 'default':
target_shp_path = make_default_path('shape_file/catchment/') # outputs a Path()
else:
target_shp_path = Path(target_shp_path) # make sure a user-specified path is a Path()
target_shp_name = read_from_control(controlFolder/controlFile,'river_basin_shp_name')
target_shp = target_shp_path / target_shp_name
target_shp_ID = read_from_control(controlFolder/controlFile,'river_basin_shp_rm_hruid')

# %% Get the source forcing file
source_nc_path = read_from_control(controlFolder/controlFile,'source_nc_path')
# Specify default path if needed
if source_nc_path == 'default':
source_nc_path = make_default_path('forcing/NA/yearly/') # outputs a Path()
else:
source_nc_path = Path(source_nc_path) # make sure a user-specified path is a Path()
source_nc_name = read_from_control(controlFolder/controlFile,'source_nc_name')
source_nc = str(source_nc_path / source_nc_name)


outdir = read_from_control(controlFolder/controlFile,'remapping_out')
if outdir == 'default':
outdir = make_default_path('forcing/{}'.format(domain_name)) # outputs a Path()
else:
outdir = Path(outdir) # make sure a user-specified path is a Path()


# %% Get the list of variable neames
var_names = read_from_control(controlFolder/controlFile,'var_names').split(', ')


# %% initializing EASYMORE object
esmr = easymore()
# specifying EASYMORE objects
# name of the case; the temporary, remapping and remapped file names include case name
esmr.case_name = 'RDRS_{}'.format(domain_name)
# temporary path that the EASYMORE generated GIS files and remapped file will be saved
esmr.temp_dir = 'temporary{}/'.format(domain_name)
# name of target shapefile that the source netcdf files should be remapped to
esmr.target_shp = target_shp
#'../shape_file/cat_pfaf_81_MERIT_Hydro_v07_Basins_v01_bugfix1_WGS84.shp'
esmr.target_shp_ID = target_shp_ID
# name of netCDF file(s); multiple files can be specified with *
esmr.source_nc = source_nc 
#'/project/6008034/baha2501/NA/forcing/NA_RDRS_v2.1_1980010112.nc'

# name of variables from source netCDF file(s) to be remapped
esmr.var_names = var_names

# rename the variables from source netCDF file(s) in the remapped files;
# it will be the same as source if not provided
#esmr.var_names_remapped = ['RDRS']
# name of variable longitude in source netCDF files
esmr.var_lon = read_from_control(controlFolder/controlFile,'var_lon')
# name of variable latitude in source netCDF files
esmr.var_lat = read_from_control(controlFolder/controlFile,'var_lat')
# name of variable time in source netCDF file; should be always time
esmr.var_time = read_from_control(controlFolder/controlFile,'var_time')
# location where the remapped netCDF file will be saved
esmr.output_dir = outdir
# format of the variables to be saved in remapped files,
# if one format provided it will be expanded to other variables
esmr.format_list = read_from_control(controlFolder/controlFile,'format_list').split(', ')
# fill values of the variables to be saved in remapped files,
# if one value provided it will be expanded to other variables
esmr.fill_value_list = read_from_control(controlFolder/controlFile,'fill_value_list').split(', ')
# if required that the remapped values to be saved as csv as well
#esmr.save_csv = True
#esmr.complevel = 9
# if uncommented EASYMORE will use this and skip GIS tasks
#esmr.remap_csv = 'temporary78/subbasin_select/RDRS_78_remapping.csv' # RDRS_81_remapping.csv

# # create source shapefile 
import geopandas as gpd
esmr.NetCDF_SHP_lat_lon()
# create the source shapefile for case 1 and 2 if shapefile is not provided
if (esmr.case == 1 or esmr.case == 2) and (esmr.source_shp == ''):
if esmr.case == 1:
if hasattr(esmr, 'lat_expanded') and hasattr(esmr, 'lon_expanded'):
esmr.lat_lon_SHP(esmr.lat_expanded, esmr.lon_expanded,\
esmr.temp_dir+esmr.case_name+'_source_shapefile.shp')
else:
esmr.lat_lon_SHP(esmr.lat, esmr.lon,\
esmr.temp_dir+esmr.case_name+'_source_shapefile.shp')
else:
esmr.lat_lon_SHP(esmr.lat, esmr.lon,\
esmr.temp_dir+esmr.case_name+'_source_shapefile.shp')
print('EASYMORE is creating the shapefile from the netCDF file and saving it here:')
print(esmr.temp_dir+esmr.case_name+'_source_shapefile.shp')

shp = gpd.read_file(esmr.temp_dir+esmr.case_name+'_source_shapefile.shp')
print(shp)
shp = shp [shp['lon_s']>-179]
shp.to_file(esmr.temp_dir+esmr.case_name+'_source_shapefile.shp')

# add the source shapefile 
esmr.source_shp = esmr.temp_dir+esmr.case_name+'_source_shapefile.shp'
esmr.source_shp_lat = 'lat_s' # name of column latitude in the source shapefile
esmr.source_shp_lon = 'lon_s' # name of column longitude in the source shapefile
esmr.source_shp_ID = 'ID_s' # name of column ID in the source shapefile

# execute EASYMORE
esmr.nc_remapper()
# %%

# --- Code provenance
# Generates a basic log file in the domain folder and copies the control file and itself there.

# Set the log path and file name
logPath = outdir
log_suffix = '_easymore_remapping.txt'

# Create a log folder
logFolder = '_workflow_log'
Path( logPath / logFolder ).mkdir(parents=True, exist_ok=True)

# Copy this script
thisFile = '2_easymore_remapping.py'
copyfile(thisFile, logPath / logFolder / thisFile);

# Get current date and time
now = datetime.now()

# Create a log file
logFile = now.strftime('%Y%m%d') + log_suffix
with open( logPath / logFolder / logFile, 'w') as file:

lines = ['Log generated by ' + thisFile + ' on ' + now.strftime('%Y/%m/%d %H:%M:%S') + '\n',
'Generated remapped climate forcing .nc file.']
for txt in lines:
file.write(txt)

5.5 Changing the forcing units

By default, the RDRSv2.1 dataset provides forcing data in units that are not acceptable for MESH. Therefore, the following script is required to ensure the correct units are used.

Climateforcing_change_unity.sh

#!/bin/bash
#SBATCH --account=rpp-kshook
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=1
#SBATCH --mem=10G
#SBATCH --time=00:30:00
#SBATCH --job-name=BowAtBanff_change_unity
#SBATCH --error=errors_BowAtBanff
#SBATCH --mail-user=example.email@usask.ca
#SBATCH --mail-type=BEGIN
#SBATCH --mail-type=END

# laod modules 
module load cdo/1.9.8

# input folder
PFAF=78
infolder=OutputBowAtBanff/ # /subbasin_select  

#fl=RDRS_"${PFAF}"_remapped_1980-01-01-13-00-00.nc
fl=BowAtBanff_remapped_1980-01-01-13-00-00.nc
echo $fl

for var in RDRS_v2.1_P_HU_09944 RDRS_v2.1_A_PR0_SFC RDRS_v2.1_P_P0_SFC RDRS_v2.1_P_FB_SFC RDRS_v2.1_P_FI_SFC RDRS_v2.1_P_TT_09944 RDRS_v2.1_P_UVC_09944 latitude longitude 
	do
		echo $var
		cdo selname,$var $infolder/$fl $infolder/$var'_'$fl		
	done

## Adjust Units
## Pressure from "mb" to "Pa"
cdo mulc,100 $infolder/'RDRS_v2.1_P_P0_SFC_'$fl $infolder/tmp.nc
cdo setattribute,RDRS_v2.1_P_P0_SFC@units=Pa $infolder/tmp.nc $infolder/'RDRS_v2.1_P_P0_SFC_'$fl
rm $infolder/tmp.nc

## Temperature from "deg_C" to "K"
cdo addc,273.16 $infolder/'RDRS_v2.1_P_TT_09944_'$fl $infolder/tmp.nc
cdo setattribute,RDRS_v2.1_P_TT_09944@units=K $infolder/tmp.nc $infolder/'RDRS_v2.1_P_TT_09944_'$fl
rm $infolder/tmp.nc

## Wind speed from "knts" to "m/s"
cdo mulc,0.5144444444444444 $infolder/'RDRS_v2.1_P_UVC_09944_'$fl $infolder/tmp.nc
cdo setattribute,RDRS_v2.1_P_UVC_09944@units="m s-1" $infolder/tmp.nc $infolder/'RDRS_v2.1_P_UVC_09944_'$fl
rm $infolder/tmp.nc

## Precipitation from "m" over the hour to a rate "mm/s" = "kg m-2 s-1"
cdo divc,3.6 $infolder/'RDRS_v2.1_A_PR0_SFC_'$fl $infolder/tmp.nc
cdo setattribute,RDRS_v2.1_A_PR0_SFC@units="mm s-1" $infolder/tmp.nc $infolder/'RDRS_v2.1_A_PR0_SFC_'$fl
rm $infolder/tmp.nc

#merge files 
cdo -z zip -b F32 merge $infolder/'RDRS_v2.1_P_HU_09944_'$fl $infolder/'RDRS_v2.1_A_PR0_SFC_'$fl $infolder/'RDRS_v2.1_P_P0_SFC_'$fl $infolder/'RDRS_v2.1_P_FB_SFC_'$fl $infolder/'RDRS_v2.1_P_FI_SFC_'$fl $infolder/'RDRS_v2.1_P_TT_09944_'$fl $infolder/'RDRS_v2.1_P_UVC_09944_'$fl $infolder/'latitude_'$fl $infolder/'longitude_'$fl $infolder/u$fl

for var in RDRS_v2.1_P_HU_09944 RDRS_v2.1_A_PR0_SFC RDRS_v2.1_P_P0_SFC RDRS_v2.1_P_FB_SFC RDRS_v2.1_P_FI_SFC RDRS_v2.1_P_TT_09944 RDRS_v2.1_P_UVC_09944 latitude longitude
	do
		rm $infolder/$var'_'$fl		
	done 

rm 	$infolder/$fl
mv $infolder/u$fl $infolder/$fl

5.4. Reordering the remapped forcing

In the final postprocessing part of the forcing dataset, it is required to reorder the forcing variables based on the remapped "Rank" IDs from the basin information file "drainage_database". Because EASYMORE remaps the forcing variables based on the MERIT Hydro catchment IDs (COMID), the order of forcing variables may not match the order of the "Rank" variable. Therefore, the fields in the remapped forcing files must be remapped to be compatible with the "drainage_database" file used for MESH. Three input data files are required for this process, the "drainage_databse", remapped forcing files, and the MERIT Hydro catchment shapefile used in the previous steps. The following section code block executes the reordering operation.

MESH_vectorbased_forcing.py

# -*- coding: utf-8 -*-
# """
# NAME
    # MESH_vectorbased_forcing 
# PURPOSE
    # The purpose of this script is to extract vector-based forcing files 
    # remapped from easymore and reorder 
    # them based rank, then save them to a netcdf format that can be read by 
    # MESH model
# PROGRAMMER(S)
    # Ala Bahrami
# REVISION HISTORY
    # 20210513 -- Initial version created and posted online
    # 20220626 -- 1) Instead of reading the shape file, the MeritHydro subbasin metadata 
                # is read as an input for reindexing forcing. 
                # 2) changed variable names to match RDRSV2.1
    
# See also 
        # easymore_extarct, create_MESH_drainage_database
        
# REFERENCES

# """
# %% load modules 
#import os
import numpy as np
import xarray as xs
import time 
from shutil import copyfile
from datetime import datetime
from pathlib import Path
#import pandas as pd
#import geopandas as gpd 

# %% define input file 
start_time = time.time()
forcing_dir     = './OutputBowAtBanff/'
forcing_name    = 'RDRS_BowAtBanff_remapped_1980-01-01-13-00-00.nc'
domain_name     = 'BowAtBanff' 
outdir          = '../workflow_data/domain_BowAtBanff/drainagedatabase/'

# %% reading input basin 
# the segids are stored in the remapped forcing, so it is not necessary to read input shape file
#basin = gpd.read_file(inut_basin)

# %% reading input netcdf files db
db = xs.open_dataset(outdir+domain_name+'_MESH_drainage_database.nc')
db.close()
segid =  db.variables['seg_id'].values
# reading for control check 
lon = db.variables['lon'].values
lat = db.variables['lat'].values

# %% reading input forcing 
forc = xs.open_dataset(forcing_dir+forcing_name)
forc.close()
lon_ease = forc.variables['longitude'].values
lat_ease = forc.variables['latitude'].values

# %% extract indices of forcing ids based on the drainage database
# Note: the input forcing is arranaged based on COMID of input MeritHydro Catchment 
# So, a reordering of forcing indices is required to match order Rank1 to RankN in MESH

n = len(segid)
ind = []

for i in range(n):
    fid = np.where(np.int32(forc['ID'].values) == segid[i])[0]
    ind = np.append(ind, fid)

ind = np.int32(ind)  

# %% reorder input forcing 
# should check if dim[0] is time 
# Note : name of variables is hard coded and can be modified regarding the input forcing 
# PR, FI','FB','TT','UV','P0','HU' 

forc_vec = xs.Dataset(
    {
        "RDRS_v2.1_A_PR0_SFC": (["subbasin", "time"], forc['RDRS_v2.1_A_PR0_SFC'].values[:,ind].transpose()),
    },
    coords={
        "time": forc['time'].values.copy(),
        "lon": (["subbasin"], lon),
        "lat": (["subbasin"], lat),
    }
    )

forc_vec['RDRS_v2.1_A_PR0_SFC'].encoding['coordinates'] = 'time lon lat'
forc_vec['RDRS_v2.1_A_PR0_SFC'].attrs["units"]          = forc['RDRS_v2.1_A_PR0_SFC'].units
forc_vec['RDRS_v2.1_A_PR0_SFC'].attrs["grid_mapping"]   = 'crs'

for n in ['RDRS_v2.1_P_FI_SFC','RDRS_v2.1_P_FB_SFC','RDRS_v2.1_P_TT_09944',
          'RDRS_v2.1_P_UVC_09944','RDRS_v2.1_P_P0_SFC','RDRS_v2.1_P_HU_09944']:
    forc_vec[n] = (("subbasin", "time"), forc[n].values[: , ind].transpose()) 
    forc_vec[n].coords["time"]          = forc['time'].values.copy()
    forc_vec[n].coords["lon"]           = (["subbasin"], lon)
    forc_vec[n].coords["lat"]           = (["subbasin"], lat)
    forc_vec[n].attrs["units"]          = forc[n].units
    forc_vec[n].attrs["grid_mapping"]   = 'crs'
    forc_vec[n].encoding['coordinates'] = 'time lon lat'

# %% update meta data attribuetes 
forc_vec.attrs['Conventions'] = 'CF-1.6'
forc_vec.attrs['License']     = 'The data were written by Ala Bahrami'
forc_vec.attrs['history']     = 'Created on June 07, 2021'
forc_vec.attrs['featureType'] = 'timeSeries'          

# editing lat attribute
forc_vec['lat'].attrs['standard_name'] = 'latitude'
forc_vec['lat'].attrs['units']         = 'degrees_north'
forc_vec['lat'].attrs['axis']          = 'Y'
 
# editing lon attribute
forc_vec['lon'].attrs['standard_name'] = 'longitude'
forc_vec['lon'].attrs['units']         = 'degrees_east'
forc_vec['lon'].attrs['axis']          = 'X'

# editing time attribute
forc_vec['time'].attrs['standard_name'] = 'time'
forc_vec['time'].attrs['axis']          = 'T'
forc_vec['time'].encoding['calendar']   = 'gregorian'
forc_vec.encoding.update(unlimited_dims = 'time')

# coordinate system
forc_vec['crs'] = db['crs'].copy()

# Define a variable for the points and set the 'timeseries_id' (required for some viewers).
forc_vec['subbasin'] = (['subbasin'], db['seg_id'].values.astype(np.int32).astype('S20'))
forc_vec['subbasin'].attrs['cf_role'] = 'timeseries_id'

#%% save to netcdf 
comp = dict(zlib=True, complevel=6)
encoding = {var: comp for var in forc_vec.data_vars}
forc_vec.to_netcdf(forcing_dir+domain_name+'_MESH_'+forcing_name,encoding=encoding)
print('--%s seconds--' %(time.time() - start_time))


# --- Code provenance
# Generates a basic log file in the domain folder and copies the control file and itself there.
 
# Set the log path and file name
logPath = Path(forcing_dir)
log_suffix = '_MESH_vectorbased_forcing.txt'
 
# Create a log folder
logFolder = '_workflow_log'
Path( logPath / logFolder ).mkdir(parents=True, exist_ok=True)
 
# Copy this script
thisFile = '4_MESH_vectorbased_forcing.py'
copyfile(thisFile, logPath / logFolder / thisFile);
 
# Get current date and time
now = datetime.now()
 
# Create a log file
logFile = now.strftime('%Y%m%d') + log_suffix
with open( logPath / logFolder / logFile, 'w') as file:
     
    lines = ['Log generated by ' + thisFile + ' on ' + now.strftime('%Y/%m/%d %H:%M:%S') + '\n',
             'Generated remapped climate forcing .nc file, reordered to MESH requirements.']
    for txt in lines:
        file.write(txt)

6. MESH configuration

The configuration input files required to run the vector-based configuration with MESH is similar to the grid-based one. The "nc_subbasin" option is used to specify the file format for both SHDFILEFLAG and BASINFORCINGFLAG. The following is an example of the MESH_input_run_options.ini file. For more information on the input files required to run MESH, consult MESH Input Files.

MESH input run options file                             # comment line 1                                | *
##### Control Flags #####                               # comment line 2                                | *
----#                                                   # comment line 3                                | *
   N                                                    # Number of control flags                       | I5
BASINFORCINGFLAG       nc_subbasin hf=60 start_date=2000010112 time_shift=-8 fname=BowAtBanff_MESH_RDRSV2_input_1980   
BASINSHORTWAVEFLAG     name_var=RDRS_v2.1_P_FB_SFC
BASINHUMIDITYFLAG      name_var=RDRS_v2.1_P_HU_09944
BASINRAINFLAG          name_var=RDRS_v2.1_A_PR0_SFC
BASINPRESFLAG          name_var=RDRS_v2.1_P_P0_SFC
BASINLONGWAVEFLAG      name_var=RDRS_v2.1_P_FI_SFC
BASINWINDFLAG          name_var=RDRS_v2.1_P_UVC_09944
BASINTEMPERATUREFLAG   name_var=RDRS_v2.1_P_TT_09944
TIMESTEPFLAG           60
INPUTPARAMSFORMFLAG    txt
SHDFILEFLAG            nc_subbasin
""
The rest of other required flags
""
##### Output Grid selection #####                      
----#                                                  
    0   #Maximum 5 points                              
---------#---------#---------#---------#---------#     
      1936                                             
         1                                             
CLASSOUT                                                # Output directory                            | 5A10
##### Output Directory #####                           
---------#                                             
results                                                 # Output Directory for total-basin files      | A10
##### Simulation Run Times #####                       
---#---#---#---#                                       
   1980   2   0   0                                        # Start year, day, hour, minute               | 4I4
   1980   365   0   0                                        # Stop year, day, hour, minute                | 4I4

7. Sample setup

A sample setup of the BowBanff Basin with required input files in addition to the results has been provided for users on the MESH GitHub page.

1. Introduction

2. Downloading Scripts & Input Data

2.1. Downloading the MESH Vector-Based Workflow Scripts

2.2. Downloading the MERIT-Hydro Shapefiles

3. Providing Settings to the Control File

Example Control File:

3.1. Root Path and Domain Name

3.2 Subbasin Selection Settings

3.3 Network Topology Settings

3.3.1. River Network File Settings

3.3.2. Catchment File Settings

3.3.3. Routing Settings

3.4 Drainage Database Settings

3.5 Climate Forcing Settings

4. Generating the Drainage Database (Basin Information File)

4.1. Subbasin Selection

Example subbasin selection script:

4.2. Network Topology

Example network topology script:

4.3. Zonal Statistics

4.3.1. Running GIS Tool for the MESH vector-based workflow

landcover_zonalhist.sh:

4.3.2. QGIS zonal histogram

4.3.2.1. Fixing invalid geometry

4.3.2.2. Obtaining the zonal histogram using the Python terminal in QGIS

lc_extract.py:

4.4. Calculate Landcover Fraction and Generate the Drainage Database (Basin Information File)

Example drainage database script:

5. Climate forcing

5.1. Extracting climate forcing

5.2. Remapping climate forcing variables

5.5 Changing the forcing units

5.4. Reordering the remapped forcing

6. MESH configuration

7. Sample setup

8. Related pages

Browser not supported