CH-Earth · CyrilThebault · Aug 3, 2024 · Aug 6, 2024 · Aug 13, 2024
diff --git a/.gitignore b/.gitignore
@@ -11,4 +11,10 @@ core.*
 # rtd
 rtd/build
 
-rtd/source/*.md
+rtd/source/*.md
+
+# Ignore macOS system files
+.DS_Store
+
+# Ignore environment configuration files
+.env
diff --git a/0_control_files/control_Africa.txt b/0_control_files/control_Africa.txt
@@ -142,8 +142,8 @@ settings_mizu_control_file  | mizuroute.control                           # Name
 settings_mizu_routing_var   | averageRoutedRunoff                         # Name of SUMMA output variable to use for routing.
 settings_mizu_routing_units | m/s                                         # Units of the variable to be routed.
 settings_mizu_routing_dt    | 3600                                        # Size of the routing time step [s].
-settings_mizu_output_freq   | annual                                      # Frequency with which mizuRoute generates new output files. Must be one of 'single', 'day', 'month', 'annual'.
-settings_mizu_output_vars   | 0                                           # Routing output. '0' for both KWT and IRF; '1' IRF only; '2' KWT only.
+settings_mizu_output_freq   | yearly                                      # Frequency with which mizuRoute generates new output files. Must be one of 'single', 'daily', 'monthly', 'yearly'.
+settings_mizu_output_vars   | 2                                           # Option for routing schemes. 0: Sum; 1: IRF; 2: KWT; 3: KW: 4: MC; 5: DW. Saves no data if not specified 
 settings_mizu_within_basin  | 0                                           # '0' (no) or '1' (IRF routing). Flag to enable within-basin routing by mizuRoute. Should be set to 0 if SUMMA is run with "subRouting" decision "timeDlay".
 settings_mizu_make_outlet   | n/a                                         # Segment ID or IDs that should be set as network outlet. Specify multiple IDs separated by commas: X,Y,Z. Specify no IDs as: n/a. Note that this can also be done in the network shapefile.
 

diff --git a/0_control_files/control_Bow_at_Banff.txt b/0_control_files/control_Bow_at_Banff.txt
@@ -142,8 +142,8 @@ settings_mizu_control_file  | mizuroute.control                           # Name
 settings_mizu_routing_var   | averageRoutedRunoff                         # Name of SUMMA output variable to use for routing.
 settings_mizu_routing_units | m/s                                         # Units of the variable to be routed.
 settings_mizu_routing_dt    | 3600                                        # Size of the routing time step [s].
-settings_mizu_output_freq   | annual                                      # Frequency with which mizuRoute generates new output files. Must be one of 'single', 'day', 'month', 'annual'.
-settings_mizu_output_vars   | 0                                           # Routing output. '0' for both KWT and IRF; '1' IRF only; '2' KWT only.
+settings_mizu_output_freq   | yearly                                      # Frequency with which mizuRoute generates new output files. Must be one of 'single', 'daily', 'monthly', 'yearly'.
+settings_mizu_output_vars   | 2                                           # Option for routing schemes. 0: Sum; 1: IRF; 2: KWT; 3: KW: 4: MC; 5: DW. Saves no data if not specified 
 settings_mizu_within_basin  | 0                                           # '0' (no) or '1' (IRF routing). Flag to enable within-basin routing by mizuRoute. Should be set to 0 if SUMMA is run with "subRouting" decision "timeDlay".
 settings_mizu_make_outlet   | 71028585                                    # Segment ID or IDs that should be set as network outlet. Specify multiple IDs separated by commas: X,Y,Z. Specify no IDs as: n/a. Note that this can also be done in the network shapefile.
 

diff --git a/0_control_files/control_Europe.txt b/0_control_files/control_Europe.txt
@@ -142,8 +142,8 @@ settings_mizu_control_file  | mizuroute.control                           # Name
 settings_mizu_routing_var   | averageRoutedRunoff                         # Name of SUMMA output variable to use for routing.
 settings_mizu_routing_units | m/s                                         # Units of the variable to be routed.
 settings_mizu_routing_dt    | 3600                                        # Size of the routing time step [s].
-settings_mizu_output_freq   | annual                                      # Frequency with which mizuRoute generates new output files. Must be one of 'single', 'day', 'month', 'annual'.
-settings_mizu_output_vars   | 0                                           # Routing output. '0' for both KWT and IRF; '1' IRF only; '2' KWT only.
+settings_mizu_output_freq   | yearly                                      # Frequency with which mizuRoute generates new output files. Must be one of 'single', 'daily', 'monthly', 'yearly'.
+settings_mizu_output_vars   | 2                                           # Option for routing schemes. 0: Sum; 1: IRF; 2: KWT; 3: KW: 4: MC; 5: DW. Saves no data if not specified 
 settings_mizu_within_basin  | 0                                           # '0' (no) or '1' (IRF routing). Flag to enable within-basin routing by mizuRoute. Should be set to 0 if SUMMA is run with "subRouting" decision "timeDlay".
 settings_mizu_make_outlet   | n/a                                         # Segment ID or IDs that should be set as network outlet. Specify multiple IDs separated by commas: X,Y,Z. Specify no IDs as: n/a. Note that this can also be done in the network shapefile.
 

diff --git a/0_control_files/control_NorthAmerica.txt b/0_control_files/control_NorthAmerica.txt
@@ -142,8 +142,8 @@ settings_mizu_control_file  | mizuroute.control                           # Name
 settings_mizu_routing_var   | averageRoutedRunoff                         # Name of SUMMA output variable to use for routing.
 settings_mizu_routing_units | m/s                                         # Units of the variable to be routed.
 settings_mizu_routing_dt    | 3600                                        # Size of the routing time step [s].
-settings_mizu_output_freq   | annual                                      # Frequency with which mizuRoute generates new output files. Must be one of 'single', 'day', 'month', 'annual'.
-settings_mizu_output_vars   | 0                                           # Routing output. '0' for both KWT and IRF; '1' IRF only; '2' KWT only.
+settings_mizu_output_freq   | yearly                                      # Frequency with which mizuRoute generates new output files. Must be one of 'single', 'daily', 'monthly', 'yearly'.
+settings_mizu_output_vars   | 2                                           # Option for routing schemes. 0: Sum; 1: IRF; 2: KWT; 3: KW: 4: MC; 5: DW. Saves no data if not specified 
 settings_mizu_within_basin  | 0                                           # '0' (no) or '1' (IRF routing). Flag to enable within-basin routing by mizuRoute. Should be set to 0 if SUMMA is run with "subRouting" decision "timeDlay".
 settings_mizu_make_outlet   | n/a                                         # Segment ID or IDs that should be set as network outlet. Specify multiple IDs separated by commas: X,Y,Z. Specify no IDs as: n/a. Note that this can also be done in the network shapefile.
 

diff --git a/0_control_files/control_NorthAsia.txt b/0_control_files/control_NorthAsia.txt
@@ -142,8 +142,8 @@ settings_mizu_control_file  | mizuroute.control                           # Name
 settings_mizu_routing_var   | averageRoutedRunoff                         # Name of SUMMA output variable to use for routing.
 settings_mizu_routing_units | m/s                                         # Units of the variable to be routed.
 settings_mizu_routing_dt    | 3600                                        # Size of the routing time step [s].
-settings_mizu_output_freq   | annual                                      # Frequency with which mizuRoute generates new output files. Must be one of 'single', 'day', 'month', 'annual'.
-settings_mizu_output_vars   | 0                                           # Routing output. '0' for both KWT and IRF; '1' IRF only; '2' KWT only.
+settings_mizu_output_freq   | yearly                                      # Frequency with which mizuRoute generates new output files. Must be one of 'single', 'daily', 'monthly', 'yearly'.
+settings_mizu_output_vars   | 2                                           # Option for routing schemes. 0: Sum; 1: IRF; 2: KWT; 3: KW: 4: MC; 5: DW. Saves no data if not specified 
 settings_mizu_within_basin  | 0                                           # '0' (no) or '1' (IRF routing). Flag to enable within-basin routing by mizuRoute. Should be set to 0 if SUMMA is run with "subRouting" decision "timeDlay".
 settings_mizu_make_outlet   | n/a                                         # Segment ID or IDs that should be set as network outlet. Specify multiple IDs separated by commas: X,Y,Z. Specify no IDs as: n/a. Note that this can also be done in the network shapefile.
 

diff --git a/0_control_files/control_Oceania.txt b/0_control_files/control_Oceania.txt
@@ -142,8 +142,8 @@ settings_mizu_control_file  | mizuroute.control                           # Name
 settings_mizu_routing_var   | averageRoutedRunoff                         # Name of SUMMA output variable to use for routing.
 settings_mizu_routing_units | m/s                                         # Units of the variable to be routed.
 settings_mizu_routing_dt    | 3600                                        # Size of the routing time step [s].
-settings_mizu_output_freq   | annual                                      # Frequency with which mizuRoute generates new output files. Must be one of 'single', 'day', 'month', 'annual'.
-settings_mizu_output_vars   | 0                                           # Routing output. '0' for both KWT and IRF; '1' IRF only; '2' KWT only.
+settings_mizu_output_freq   | yearly                                      # Frequency with which mizuRoute generates new output files. Must be one of 'single', 'daily', 'monthly', 'yearly'.
+settings_mizu_output_vars   | 2                                           # Option for routing schemes. 0: Sum; 1: IRF; 2: KWT; 3: KW: 4: MC; 5: DW. Saves no data if not specified 
 settings_mizu_within_basin  | 0                                           # '0' (no) or '1' (IRF routing). Flag to enable within-basin routing by mizuRoute. Should be set to 0 if SUMMA is run with "subRouting" decision "timeDlay".
 settings_mizu_make_outlet   | n/a                                         # Segment ID or IDs that should be set as network outlet. Specify multiple IDs separated by commas: X,Y,Z. Specify no IDs as: n/a. Note that this can also be done in the network shapefile.
 

diff --git a/0_control_files/control_SouthAmerica.txt b/0_control_files/control_SouthAmerica.txt
@@ -142,8 +142,8 @@ settings_mizu_control_file  | mizuroute.control                           # Name
 settings_mizu_routing_var   | averageRoutedRunoff                         # Name of SUMMA output variable to use for routing.
 settings_mizu_routing_units | m/s                                         # Units of the variable to be routed.
 settings_mizu_routing_dt    | 3600                                        # Size of the routing time step [s].
-settings_mizu_output_freq   | annual                                      # Frequency with which mizuRoute generates new output files. Must be one of 'single', 'day', 'month', 'annual'.
-settings_mizu_output_vars   | 0                                           # Routing output. '0' for both KWT and IRF; '1' IRF only; '2' KWT only.
+settings_mizu_output_freq   | yearly                                      # Frequency with which mizuRoute generates new output files. Must be one of 'single', 'daily', 'monthly', 'yearly'.
+settings_mizu_output_vars   | 2                                           # Option for routing schemes. 0: Sum; 1: IRF; 2: KWT; 3: KW: 4: MC; 5: DW. Saves no data if not specified 
 settings_mizu_within_basin  | 0                                           # '0' (no) or '1' (IRF routing). Flag to enable within-basin routing by mizuRoute. Should be set to 0 if SUMMA is run with "subRouting" decision "timeDlay".
 settings_mizu_make_outlet   | n/a                                         # Segment ID or IDs that should be set as network outlet. Specify multiple IDs separated by commas: X,Y,Z. Specify no IDs as: n/a. Note that this can also be done in the network shapefile.
 

diff --git a/0_control_files/control_SouthAsia.txt b/0_control_files/control_SouthAsia.txt
@@ -142,8 +142,8 @@ settings_mizu_control_file  | mizuroute.control                           # Name
 settings_mizu_routing_var   | averageRoutedRunoff                         # Name of SUMMA output variable to use for routing.
 settings_mizu_routing_units | m/s                                         # Units of the variable to be routed.
 settings_mizu_routing_dt    | 3600                                        # Size of the routing time step [s].
-settings_mizu_output_freq   | annual                                      # Frequency with which mizuRoute generates new output files. Must be one of 'single', 'day', 'month', 'annual'.
-settings_mizu_output_vars   | 0                                           # Routing output. '0' for both KWT and IRF; '1' IRF only; '2' KWT only.
+settings_mizu_output_freq   | yearly                                      # Frequency with which mizuRoute generates new output files. Must be one of 'single', 'daily', 'monthly', 'yearly'.
+settings_mizu_output_vars   | 2                                           # Option for routing schemes. 0: Sum; 1: IRF; 2: KWT; 3: KW: 4: MC; 5: DW. Saves no data if not specified 
 settings_mizu_within_basin  | 0                                           # '0' (no) or '1' (IRF routing). Flag to enable within-basin routing by mizuRoute. Should be set to 0 if SUMMA is run with "subRouting" decision "timeDlay".
 settings_mizu_make_outlet   | n/a                                         # Segment ID or IDs that should be set as network outlet. Specify multiple IDs separated by commas: X,Y,Z. Specify no IDs as: n/a. Note that this can also be done in the network shapefile.
 

diff --git a/0_example/shapefiles/catchment/bow_distributed_elevation_zone.dbf b/0_example/shapefiles/catchment/bow_distributed_elevation_zone.dbf
diff --git a/0_example/shapefiles/catchment/bow_distributed_elevation_zone.shp b/0_example/shapefiles/catchment/bow_distributed_elevation_zone.shp
diff --git a/0_example/shapefiles/catchment/bow_distributed_elevation_zone.shx b/0_example/shapefiles/catchment/bow_distributed_elevation_zone.shx
diff --git a/4b_remapping/1_topo/1_find_HRU_elevation.ipynb b/4b_remapping/1_topo/1_find_HRU_elevation.ipynb
@@ -5,10 +5,10 @@
    "metadata": {},
    "source": [
     "# Intersect catchment with MERIT DEM\n",
-    "Finds the mean elevation of each HRU in the model setup with pyQGIS.\n",
+    "Finds the mean elevation of each HRU in the model setup with rasterstats.\n",
     "\n",
     "### Note\n",
-    "The pyQGIS function `QgsZonalStatistics` automatically adds the calculated value to the shapefile used as input to the function. The workflow is thus:\n",
+    "The workflow is thus:\n",
     "1. Find the source catchment shapefile;\n",
     "2. Copy the source catchment shapefile to the destintion location;\n",
     "3. Run the zonal statistics algorithm on the copy."
@@ -22,12 +22,11 @@
    "source": [
     "# modules\n",
     "import os\n",
+    "import geopandas as gpd\n",
+    "from rasterstats import zonal_stats\n",
     "from pathlib import Path\n",
     "from shutil import copyfile\n",
-    "from datetime import datetime\n",
-    "from qgis.core import QgsVectorLayer\n",
-    "from qgis.core import QgsRasterLayer\n",
-    "from qgis.analysis import QgsZonalStatistics"
+    "from datetime import datetime"
    ]
   },
   {
@@ -241,43 +240,40 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "#### QGIS analysis"
+    "#### Rasterstats analysis"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 15,
+   "execution_count": 16,
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Convert Path() to string for QGIS\n",
-    "catchment_file = str(intersect_path/intersect_name) # needs to be the coped file because output is automatically added to this\n",
-    "dem_file = str(dem_path/dem_name)"
+    "# Load the shapefile\n",
+    "gdf = gpd.read_file(str(intersect_path/intersect_name))"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": 17,
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Load the shape and raster\n",
-    "layer_polygon = QgsVectorLayer(catchment_file,'merit_hydro_basin','ogr') # this works\n",
-    "layer_raster  = QgsRasterLayer(dem_file,'merit_hydro_dem') # this works"
+    "# Calculate zonal statistics\n",
+    "stats = zonal_stats(gdf, \n",
+    "                    str(dem_path/dem_name), \n",
+    "                    stats=['mean'], \n",
+    "                    all_touched=True)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 17,
+   "execution_count": 18,
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Check we loaded the layers correctly\n",
-    "if not layer_raster.isValid():\n",
-    "    print('Raster layer failed to load')\n",
-    "    \n",
-    "if not layer_polygon.isValid():\n",
-    "    print('Polygon layer failed to load')"
+    "# Add the mean elevation to the GeoDataFrame\n",
+    "gdf['elev_mean'] = [stat['mean'] for stat in stats]"
    ]
   },
   {
@@ -286,13 +282,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Create a zonal statistics object, automatically saved to file\n",
-    "band = 1 # raster band with the data we are after\n",
-    "zonalstats = QgsZonalStatistics(layer_polygon,                 # shapefile\n",
-    "                                layer_raster,                  # .tif\n",
-    "                                'elev_',                       # prefix for the new column added to the shapefile  \n",
-    "                                band,                          # raster band we're interested in\n",
-    "                                stats=QgsZonalStatistics.Mean).calculateStatistics(None)"
+    "# Save the updated GeoDataFrame\n",
+    "gdf.to_file(str(intersect_path/intersect_name))"
    ]
   },
   {

diff --git a/4b_remapping/1_topo/1_find_HRU_elevation.py b/4b_remapping/1_topo/1_find_HRU_elevation.py
@@ -1,22 +1,18 @@
 # Intersect catchment with MERIT DEM
-# Finds the mean elevation of each HRU in the model setup with pyQGIS.
+# Finds the mean elevation of each HRU in the model setup with rasterstats.
 #
 # Note:
-# The pyQGIS function `QgsZonalStatistics` automatically adds the calculated value to the shapefile used as input to the function. The workflow is thus:
 # 1. Find the source catchment shapefile;
 # 2. Copy the source catchment shapefile to the destintion location;
 # 3. Run the zonal statistics algorithm on the copy.
 
 # modules
 import os
+import geopandas as gpd
+from rasterstats import zonal_stats
 from pathlib import Path
-from shutil import which
 from shutil import copyfile
 from datetime import datetime
-from qgis.core import QgsApplication
-from qgis.core import QgsVectorLayer
-from qgis.core import QgsRasterLayer
-from qgis.analysis import QgsZonalStatistics
 
 
 # --- Control file handling
@@ -115,39 +111,22 @@ def make_default_path(suffix):
         copyfile(catchment_path/file, intersect_path/newfile);
 
 
-# --- QGIS analysis
-# Initialize QGIS
-os.environ["QT_QPA_PLATFORM"] = "offscreen" # disable QT trying to connect to display; needed on HPC infrastructure
-qgis_path = which('qgis') # find the QGIS install location
-QgsApplication.setPrefixPath(qgis_path, True) # supply path to qgis install location
-qgs = QgsApplication([], False) # create a reference to the QgsApplication, GUI = False
-qgs.initQgis() # load providers
-
-# Convert Path() to string for QGIS
-catchment_file = str(intersect_path/intersect_name) # needs to be the copied file because output is automatically added to this
-dem_file = str(dem_path/dem_name)
-
-# Load the shape and raster
-layer_polygon = QgsVectorLayer(catchment_file,'merit_hydro_basin','ogr')
-layer_raster  = QgsRasterLayer(dem_file,'merit_hydro_dem')
-
-# Check we loaded the layers correctly
-if not layer_raster.isValid():
-    print('Raster layer failed to load')
-
-if not layer_polygon.isValid():
-    print('Polygon layer failed to load')
-
-# Create a zonal statistics object, automatically saved to file
-band = 1 # raster band with the data we are after
-zonalstats = QgsZonalStatistics(layer_polygon,                 # shapefile
-                                layer_raster,                  # .tif
-                                'elev_',                       # prefix for the new column added to the shapefile  
-                                band,                          # raster band we're interested in
-                                stats=QgsZonalStatistics.Mean).calculateStatistics(None)
+# --- Rasterstats analysis
+# Load the shapefile
+gdf = gpd.read_file(str(intersect_path/intersect_name))
+
+# Calculate zonal statistics
+stats = zonal_stats(gdf, 
+                    str(dem_path/dem_name), 
+                    stats=['mean'], 
+                    all_touched=True)
+
+# Add the mean elevation to the GeoDataFrame
+gdf['elev_mean'] = [stat['mean'] for stat in stats]
+
+# Save the updated GeoDataFrame
+gdf.to_file(str(intersect_path/intersect_name))
 
-# Clean memory
-qgs.exitQgis()
 
 
 # --- Code provenance