Merge pull request #180 from changliao1025/development

Update before opt
changliao1025 · Sep 25, 2023 · a0b58a3 · a0b58a3
2 parents d46616d + 16a3fff
commit a0b58a3
Show file tree

Hide file tree

Showing 51 changed files with 38,877 additions and 8,745 deletions.
diff --git a/.gitignore b/.gitignore
@@ -135,4 +135,7 @@ retired
 core.*
 *.nc
 *.tif
-*.tiff
+*.tiff
+
+# mac file system
+.DS_Store
diff --git a/README.md b/README.md
@@ -3,20 +3,20 @@
 [![DOI](https://zenodo.org/badge/368338554.svg)](https://zenodo.org/badge/latestdoi/368338554)
 [![Downloads](https://static.pepy.tech/badge/pyflowline)](https://pepy.tech/project/pyflowline)
 
-PyFlowline: a mesh-independent river network generator for hydrologic models. 
+PyFlowline: a mesh-independent river network generator for hydrologic models.  
 
 ### Quickstart
 
 Please refer to the [quickstart documentation](https://pyflowline.readthedocs.io/en/latest/quickstart.html) for details on how to get started using the PyFlowline package.
 
-PyFlowline is mesh independent, meaning you can apply it to both structured 
+PyFlowline is mesh independent, meaning you can apply it to both structured  
 
-1. traditional rectangle projected mesh 
-2. latitude-longitude 
+1. traditional rectangle projected mesh  
+2. latitude-longitude  
 3. hexagon
 4. dggs ([dggrid](https://github.com/sahrk/DGGRID))
 
-and unstructured mesh systems 
+and unstructured mesh systems  
 
 1. Model for Prediction Across Scales mesh ([MPAS](https://github.com/MPAS-Dev))
 2. Triangulated Irregular Network (TIN) mesh
@@ -56,7 +56,7 @@ This work was supported by the Earth System Model Development program areas of t
 
 This research was supported as part of the Next Generation Ecosystem Experiments-Tropics, funded by the U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research at Pacific Northwest National Laboratory. The study was also partly supported by U.S. Department of Energy Office of Science Biological and Environmental Research through the Earth and Environmental System Modeling program as part of the Energy Exascale Earth System Model (E3SM) project.
 
-This research was supported by the Next Generation Ecosystem Experiments-Tropics project, funded by the U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research at Pacific Northwest National Laboratory. 
+This research was supported by the Next Generation Ecosystem Experiments-Tropics project, funded by the U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research at Pacific Northwest National Laboratory.  
 
 ### License
 

diff --git a/docs/source/data/data.rst b/docs/source/data/data.rst
diff --git a/docs/source/quickstart.rst b/docs/source/quickstart.rst
@@ -12,7 +12,11 @@ Users can run a PyFlowline simulation in the following steps:
 2. Install the package using `conda install -c conda-forge pyflowline`. Conda will automatically install all the required dependencies.
 3. Clone the latest PyFlowline repository from https://github.com/changliao1025/pyflowline. 
 4. Download the additional large MPAS mesh file `lnd_cull_mesh.nc` from https://github.com/changliao1025/pyflowline/releases/tag/0.2.0 and move it under the `data/susquehanna/input` folder.
-5. Open the `examples/susquehanna/pyflowline_susquehanna_mpas.json` file and change `sWorkspace_output` to the full path to the directory where you want to save the output (e.g. `/full/path/to/pyflowline/data/susquehanna/output`), change `"sFilename_mesh_netcdf"` to the full path to `lnd_cull_mesh.nc`, `"sFilename_mesh_boundary"` to the full path to `data/susquehanna/input/boundary_wgs.geojson`, and `"sFilename_basins"` to the full path to `examples/susquehanna/pyflowline_susquehanna_basins.json`.
+5. Open the `examples/susquehanna/pyflowline_susquehanna_mpas.json` file and make the following changes:
+   - Change `sWorkspace_output` to the full path to the directory where you want to save the output (e.g. `/full/path/to/pyflowline/data/susquehanna/output`).
+   - Change `"sFilename_mesh_netcdf"` to the full path to `lnd_cull_mesh.nc`.
+   - Change `"sFilename_mesh_boundary"` to the full path to `data/susquehanna/input/boundary_wgs.geojson`.
+   - Change `"sFilename_basins"` to the full path to `examples/susquehanna/pyflowline_susquehanna_basins.json`.
 6. Open the `examples/susquehanna/pyflowline_susquehanna_basins.json` file and change `"sFilename_flowline_filter"` to the full path to `data/susquehanna/input/flowline.geojson`. Ignore the other settings in these json files for now.
 7. Open the preferred Python IDE (Visual Studio Code recommended) and run the  `examples/susquehanna/run_simulation_mpas.py` Python script. Optionally, you can also run the `notebooks/mpas_example.ipynb` notebook. The visualization of the model outputs is only experimental, and you can use other tools to visualize the model outputs.
 8. You should produce a list of model outputs in the `data/susquehanna/output` folder or the user-specified output folder.

diff --git a/examples/README.md b/examples/README.md
@@ -1,5 +1,7 @@
-In this directory, you will find several different example files for different purposes.
+# Examples
 
-- create_model_condfiguration.py: this example, create a configuration file, you can edit this configuration and then use it for the actual simulation such as in the run_simulation.py
+This directory contains several different example files that demonstrate how to use the PyFlowline software.
 
-- run_simulation.py: this example uses a configuration file to run a simulation
+- create_model_configuration.py: This example creates a configuration file. You can edit this configuration and then use it for the actual simulation such as in the run_simulation.py scripts.
+
+- run_simulation_<meshtype>.py: These examples demonstrate how to read a configuration file and use the generated PyFlowline object to run a simulation for four types of meshes.
diff --git a/examples/create_model_condfiguration.py → examples/create_model_configuration.py b/examples/create_model_condfiguration.py → examples/create_model_configuration.py
diff --git a/notebooks/mpas_example.ipynb b/notebooks/mpas_example.ipynb
@@ -59,17 +59,17 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "#check pyflowline\n",
+    "#check pyflowline installation\n",
     "iFlag_pyflowline = importlib.util.find_spec(\"pyflowline\") \n",
     "if iFlag_pyflowline is not None:\n",
     "    pass\n",
     "else:\n",
     "    print('The pyflowline package is not installed. Please install it following the quickstart document.')\n",
     "\n",
     "#check optional packages    \n",
-    "iFlag_cython = importlib.util.find_spec(\"cython\") \n",
-    "iFlag_cartopy = importlib.util.find_spec(\"cartopy\") \n",
-    "iFlag_geopandas = importlib.util.find_spec(\"geopandas\") \n",
+    "iFlag_cython = importlib.util.find_spec(\"cython\")  #check whether cython is installed\n",
+    "iFlag_cartopy = importlib.util.find_spec(\"cartopy\") #check whether cartopy is installed\n",
+    "iFlag_geopandas = importlib.util.find_spec(\"geopandas\") #check whether geopandas is installed\n",
     "if iFlag_cartopy is not None:\n",
     "    iFlag_cartopy = 1\n",
     "    pass\n",
@@ -78,6 +78,7 @@
     "    if iFlag_geopandas is not None:\n",
     "        pass\n",
     "    else:\n",
+    "        #if both cartopy and geopandas are not available, we will install the geopandas package under the current environment\n",
     "        print('We will install the geopandas package for visualization.')\n",
     "        !conda install --yes --prefix {sys.prefix} gepandas\n",
     "        iFlag_geopandas = 1\n",
@@ -166,6 +167,7 @@
    "outputs": [],
    "source": [
     "\n",
+    "#download the MPAS mesh from the github release\n",
     "sFilename_mpas = 'https://github.com/changliao1025/pyflowline/releases/download/0.2.0/lnd_cull_mesh.nc'\n",
     "\n",
     "#combind folder with filename to get the full path\n",
@@ -318,7 +320,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "#step 5\n",
+    "#step 5, set up some parameters\n",
     "sMesh_type = 'mpas'\n",
     "iCase_index = 1\n",
     "dResolution_meter=5000\n",
@@ -349,6 +351,11 @@
     }
    ],
    "source": [
+    "#the read function accepts several keyword arguments that can be used to change the default parameters.\n",
+    "#the normal keyword arguments are:\n",
+    "#iCase_index_in: this is an ID to identify the simulation case\n",
+    "#sMesh_type_in: this specifies the mesh type ('mpas' in this example)\n",
+    "#sDate_in: this specifies the date of the simulation, the final output folder will have a pattern such as 'pyflowline20230901001', where pyflowline is model, 20230901 is the date, and 001 is the case index.\n",
     "oPyflowline = pyflowline_read_model_configuration_file(sFilename_configuration_in, iCase_index_in=iCase_index, \n",
     "                                                    sMesh_type_in= sMesh_type, sDate_in=sDate)"
    ]
@@ -360,7 +367,32 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "oPyflowline.change_model_parameter('sFilename_mesh_netcdf', sFilename_download)"
+    "#other than setting the paraemeter using the read model configuration fucntion, user can also change model parameters after creating the model object\n",
+    "#rememere that, it is recommended to set output folder using the read model configuration function since the change_model_parameter function will not update output folder\n",
+    "#only a list of parameters can be changed, for full list, please check the documentation\n",
+    "#in this example, we will change the mesh file name\n",
+    "#the function will check the data type, if incorrect data type is provided, it will raise an error\n",
+    "oPyflowline.change_model_parameter('sFilename_mesh_netcdf', sFilename_download) #because the mpas mesh already contains elevation, we do not need to set the elevation file name\n",
+    "#we will set the boundary file name\n",
+    "#this file should be located in the input folder\n",
+    "sFolder_data = os.path.join(sPath_parent, 'data')\n",
+    "sFolder_data_susquehanna =  os.path.join(sFolder_data, 'susquehanna')\n",
+    "sFolder_input = os.path.join(sFolder_data_susquehanna, 'input')\n",
+    "sFilename_mesh_boundary = os.path.join(sPath_parent, 'boundary_wgs.geojson')\n",
+    "oPyflowline.change_model_parameter('sFilename_boundary_netcdf', sFilename_mesh_boundary)\n",
+    "\n",
+    "#we will set the original flowline file name\n",
+    "sFilename_flowline = os.path.join(sFolder_input, 'flowline.geojson')\n",
+    "oPyflowline.change_model_parameter('sFilename_flowline_filtered', sFilename_flowline)\n",
+    "\n",
+    "#we can also set for individual basin in the domain, in this example, we only has one basin.\n",
+    "#remember that, each basin can have different parameters, so if you want to set them different (for example, basin 1 has no dam, but basin 2 has dam), you should edit the basin json instead using this function.\n",
+    "#because change_model_parameter will set all the basin using the same parameter in current version\n",
+    "#must set iFlag_basin_in = 1 for basin parameter\n",
+    "oPyflowline.change_model_parameter('iFlag_dam', 0, iFlag_basin_in= 1)\n",
+    "\n",
+    "\n",
+    "\n"
    ]
   },
   {
@@ -443,8 +475,30 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "#another important setting for basin is the approximate outlet location\n",
     "oPyflowline.aBasin[0].dLatitude_outlet_degree=39.4620\n",
-    "oPyflowline.aBasin[0].dLongitude_outlet_degree=-76.0093"
+    "oPyflowline.aBasin[0].dLongitude_outlet_degree=-76.0093\n",
+    "#you can also set it using the change_model_parameter function\n",
+    "oPyflowline.change_model_parameter('dLongitude_outlet_degree', -76.0093, iFlag_basin_in = 1)\n",
+    "oPyflowline.change_model_parameter('dLatitude_outlet_degree', 39.4620, iFlag_basin_in = 1)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "61f5b8ed",
+   "metadata": {},
+   "source": [
+    "You can check the setting for the single basin as well"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "edd17e3e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(oPyflowline.aBasin[0].tojson())"
    ]
   },
   {
@@ -515,7 +569,7 @@
     "    oPyflowline.plot(sVariable_in = 'flowline_filter', sFilename_in = 'filter_flowline.png' )\n",
     "    pass\n",
     "else: #use the default visualization method, only experimental\n",
-    "    if iFlag_geopandas = 1:\n",
+    "    if iFlag_geopandas == 1:\n",
     "        import geopandas as gpd\n",
     "        import matplotlib.pyplot as plt\n",
     "        #use the geopanda package\n",

diff --git a/pyflowline/algorithms/auxiliary/check_head_water.py b/pyflowline/algorithms/auxiliary/check_head_water.py
@@ -1,5 +1,6 @@
 def check_head_water(aFlowline_in, pVertex_start_in):
     """[Check whether a vertex assoacited with a flowline is a headwater or not]
+    if the stream order info is avaialble, use the stream order info to check is easier and faster
 
     Args:
         aFlowline_in ([pyflowline]): [all the flowline]

diff --git a/pyflowline/algorithms/auxiliary/find_index_in_list.py b/pyflowline/algorithms/auxiliary/find_index_in_list.py
@@ -6,8 +6,10 @@
 iFlag_cython = importlib.util.find_spec("cython") 
 if iFlag_cython is not None:
     from pyflowline.algorithms.cython.kernel import find_vertex_in_list
+    from pyflowline.external.tinyr.tinyr.tinyr import RTree
+    iFlag_use_rtree = 1
 else:
-
+    iFlag_use_rtree =0
     from pyflowline.algorithms.auxiliary.find_vertex_in_list import find_vertex_in_list
 
 
@@ -19,27 +21,69 @@ def find_vertex_on_edge(aVertex_in, pEdge_in):
     nVertex= len(aVertex_in)
     npoint = 0    
     if nVertex > 0 :
-        for i in np.arange( nVertex):
-            pVertex = aVertex_in[i]
-            iFlag_overlap, dDistance, diff = pEdge_in.check_vertex_on_edge(pVertex)
-            if iFlag_overlap == 1:                
-                iFlag_exist = 1                      
-                aDistance.append(dDistance)
-                aIndex.append(i) 
-                npoint = npoint + 1          
-            else:                
-                if  diff < 1.0:
-                    iFlag_overlap = pEdge_in.check_vertex_on_edge(pVertex)
-                   
+        if iFlag_use_rtree == 1:
+            index_vertex = RTree(max_cap=5, min_cap=2)
+            for i in range(nVertex):
+                lID = i 
+                x = aVertex_in[i].dLongitude_degree
+                y = aVertex_in[i].dLatitude_degree   
+                left= x - 1E-5
+                right= x + 1E-5
+                bottom= y-1E-5
+                top=    y+1E-5
+                pBound= (left, bottom, right, top)
+                index_vertex.insert(lID, pBound)  # 
                 pass
+            #now the new vertex
+            pVertex_start = pEdge_in.pVertex_start
+            pVertex_end = pEdge_in.pVertex_end
+            x1=pVertex_start.dLongitude_degree
+            y1=pVertex_start.dLatitude_degree   
+            x2=pVertex_end.dLongitude_degree
+            y2=pVertex_end.dLatitude_degree
+            left   = np.min([x1, x2])
+            right  = np.max([x1, x2])
+            bottom = np.min([y1, y2])
+            top    = np.max([y1, y2])
+            pBound= (left, bottom, right, top)
+            aIntersect = list(index_vertex.search(pBound))
+            for k in aIntersect:
+                pVertex = aVertex_in[k]
+                iFlag_overlap, dDistance, diff = pEdge_in.check_vertex_on_edge(pVertex)
+                if iFlag_overlap == 1:                
+                    iFlag_exist = 1                      
+                    aDistance.append(dDistance)
+                    aIndex.append(k) 
+                    npoint = npoint + 1          
+                else:                
+                    if diff < 1.0:
+                        iFlag_overlap = pEdge_in.check_vertex_on_edge(pVertex)
+
+                    pass
+
 
-        #re-order 
-        if iFlag_exist == 1 :
-            x = np.array(aDistance)
-            b = np.argsort(x)
-            c = np.array(aIndex)
-            d= c[b]
-            aIndex_order = list(d)        
+        else:
+            for i in np.arange( nVertex):
+                pVertex = aVertex_in[i]
+                iFlag_overlap, dDistance, diff = pEdge_in.check_vertex_on_edge(pVertex)
+                if iFlag_overlap == 1:                
+                    iFlag_exist = 1                      
+                    aDistance.append(dDistance)
+                    aIndex.append(i) 
+                    npoint = npoint + 1          
+                else:                
+                    if  diff < 1.0:
+                        iFlag_overlap = pEdge_in.check_vertex_on_edge(pVertex)
+
+                    pass
+
+            #re-order 
+            if iFlag_exist == 1 :
+                x = np.array(aDistance)
+                b = np.argsort(x)
+                c = np.array(aIndex)
+                d= c[b]
+                aIndex_order = list(d)        
     else:
         pass
 
@@ -171,6 +215,7 @@ def add_unique_vertex(aVertex_in, pVertex_in, dThreshold_in = 1.0E-6):
         pass
 
     return aVertex_in, iFlag_exist
+
 def find_list_in_list(aList_in, pList_in):
     c = copy.deepcopy(pList_in)
     c.sort()