Development

README.md

@@ -5,9 +5,13 @@
 
 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 
 
-1. traditional rectangle mesh 
+1. traditional rectangle projected mesh 
 2. latitude-longitude 
 3. hexagon
 4. dggs ([dggrid](https://github.com/sahrk/DGGRID))
@@ -23,27 +27,21 @@ This package generates the mesh cell-based conceptual river networks using the f
 2. `Mesh generation`: PyFlowline generates structured meshes (e.g., rectangle, hexagon) or imports user-provided unstructured meshes into the PyFlowline-compatible GEOJSON format.
 3. `Topological relationship reconstruction`: PyFlowline reconstructs the topological relationship using the mesh and flowline intersections.
 
-
 ### Dependency
 
 PyFlowline depends on the following packages
 
 1. `numpy`
 2. `gdal`
 3. `netCDF4`
-4. `shapely`
 
 PyFlowline also has three optional dependency packages
 
 1. `cython` for performance 
-2. `matplotlin` for visualization
+2. `matplotlib` for visualization
 3. `cartopy` for visulization
 4. `simplekml` for Google Earth KML support
 
-### 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.
-
 ### Installation
 
 Please refer to the [official documentation](https://pyflowline.readthedocs.io/) for details on how to install the PyFlowline package.

conda-recipe/conda_build_config.yaml

@@ -23,7 +23,7 @@ metadata:
   # Package name
   name: pyflowline
   # Package version
-  version: "0.2.5"
+  version: "0.2.6"
   # Package summary
   summary: A mesh-independent river network generator for hydrologic models.
   # Package homepage
@@ -39,7 +39,6 @@ metadata:
     - python >=3.8
     - numpy
     - gdal
-    - shapely
     - netCDF4
     - matplotlib-base
   # Package build dependencies

conda-recipe/meta.yaml

@@ -1,6 +1,6 @@
 {% set name = "hexwatershed" %}
 {% set git_rev = "main" %}
-{% set version = "0.2.5" %}
+{% set version = "0.2.6" %}
 
 package:
   name: {{ name|lower }}

data/susquehanna/readme.md

@@ -0,0 +1,5 @@
+In this directory, you will find the input and output of an example PyFlowline simulation. 
+
+- input The input files for a MPAS mesh-based simulation.
+
+- output The output file from a MPAS mesh-based simulation. 

docs/source/data/data.rst

@@ -48,10 +48,10 @@ File I/O
 Inputs
 ==============================
 
-
-PyFlowline uses two configuration files to manage all the input information. Within this configuration file, it stores major model input parameters and paths to input files. 
+PyFlowline uses two configuration files to manage all the input information, within which major model input parameters and paths are specified. 
 
 These two configuration files have a parent-child relationship:
+
 1. The parent configuration file stores parameters for the whole domain, and
 2. The child configuration file stores parameters for every single watershed.
 

docs/source/glossary.rst

@@ -3,6 +3,30 @@ Glossary
 ########
 
 
+*****************
+Structured mesh
+*****************
+
+In PyFlowline, structured mesh refers to meshes that have a repeating pattern or structure.
+
+The following meshes are considered as structured:
+
+1. Projected raster meshes (e.g. 100m by 100m)
+2. GCS-based rectangle meshes  (e.g. 0.5 degree by 0.5 degree)
+3. Hexagon meshes (e.g. 100m by edge)
+4. DGGS meshes (e.g., DGGrid meshes)
+
+*****************
+Unstructured mesh
+*****************
+
+In PyFlowline, unstructured mesh refers to meshes that don't have a repeating pattern or structure and the cell size varies from cell to cell.
+
+The following meshes are considered as unstructured:
+
+1. Model for Prediction Across Scales (MPAS) meshes
+2. Triangulated irregular network (TIN) meshes
+
 ************
 Great circle
 ************

docs/source/installation/installation.rst

@@ -31,6 +31,16 @@ Option A
 In this option, you will use conda to install the released PyFlowline package, but not necessarily the latest version.
 Conda will automatically install all the dependency packages.
 
+Before you install the package, it is highly recommended that you start from a new conda environment using the following command:
+
+    conda create -n pyflowline_test
+
+After activating the environment with:
+
+    conda activate pyflowline_test
+
+You can then install it with:
+
     conda install -c conda-forge pyflowline
 
 
@@ -51,8 +61,6 @@ The following dependency packages will be installed during the process.
 * `numpy`
 * `gdal`
 * `netCDF4`
-* `shapely`
-
 
 =============
 Visualization

docs/source/quickstart.rst

@@ -2,6 +2,10 @@
 Quickstart
 #####################
 
+Installing and running PyFlowline requires some basic knowledge of the Python ecosystem.
+
+Besides, configuring a PyFlowline simulation requires some knowledge of Geographic Information System (GIS) and computational hydrology.
+
 Users can run a PyFlowline simulation in the following steps:
 
 1. Create a new Python environment using Conda, and activate the new environment.

docs/source/readme.rst

@@ -8,11 +8,11 @@ Overview
 
 PyFlowline is a mesh-independent river network generator for hydrologic models.
 
-Traditionally, the river network is represented by the vector-based river network. These vector-based river networks usually cannot be directly used by hydrologic models, which use the mesh as the spatial discretization. 
+River networks are landscape features typically represented using vector layers. However, most hydrologic models rely on regular grids to discretize the spatial domain and cannot directly ingest vector features into the model. As a result, hydrologic models usually implement a so-called stream-burning process to convert the vector-based river network into a mesh-based river network. 
 
-In most cases, the vector-based river network datasets are used in a so-called stream-burning process to convert the vector-based river network to the mesh-based river network. However, all the existing stream-burning methods only support the structured meshes and there are also some other limitations.
+However, all the existing stream-burning methods only support the structured meshes and there are also some other limitations. For example, existing stream-burning methods always treat the vector river networks as a binary mask and cannot describe the topology near river confluences and meanders. 
 
-To close this gap, PyFlowline was developed using a mesh-independent approach. At its core, PyFlowline uses the intersection between the vector river network and mesh to reconstruct the conceptual river network. It also addresses several limitations of the existing stream-burning methods.
+PyFlowline solves this issue by using a mesh-independent approach that intersects the vector river network and mesh to reconstruct the conceptual river network.
 
 ***********
 Development
@@ -31,6 +31,13 @@ As a result, if a spatially-distributed hydrologic model uses the unstructured m
 To close this gap, PyFlowline was developed using a mesh-independent approach. At its core, PyFlowline uses the intersection between the vector river network and mesh to reconstruct the conceptual river network.
 
 
+*****************
+Target audience
+*****************
+
+PyFlowline is an advanced modeling tool for hydrologists and hydrologic modelers. 
+Users of PyFlowline should be familiar with basic concepts in Geographic Information System (GIS), including vector and raster data, coordinate systems, and projections.
+
 *****************
 Important notice
 *****************

examples/susquehanna/run_simulation_mpas.py

@@ -27,6 +27,7 @@
 #setup case information
 #===================================
 iCase_index = 1
+iFlag_simulation = 0
 iFlag_visualization = 1
 sMesh = 'mpas'
 sDate='20230701'
@@ -38,28 +39,32 @@
 oPyflowline.aBasin[0].dLatitude_outlet_degree=39.462000
 oPyflowline.aBasin[0].dLongitude_outlet_degree=-76.009300
 
-oPyflowline.setup()
+#oPyflowline.setup()
 if iFlag_visualization ==1:
     #oPyflowline.plot(sVariable_in = 'flowline_filter', sFilename_output_in = 'filter_flowline.png'  )
     pass
 
-oPyflowline.flowline_simplification()
+if iFlag_simulation == 1:
+    #oPyflowline.flowline_simplification()
+    pass
 
 if iFlag_visualization == 1:
 
     aExtent_meander = [-76.5,-76.2, 41.6,41.9] 
-    oPyflowline.plot( sVariable_in='flowline_simplified' , sFilename_output_in = 'flowline_simplified.png' ) 
-
-    oPyflowline.plot( sVariable_in='flowline_simplified' , sFilename_output_in = 'flowline_simplified_zoom.png', aExtent_in =aExtent_meander ) 
+    #oPyflowline.plot( sVariable_in='flowline_simplified' , sFilename_output_in = 'flowline_simplified.png' ) 
+    #oPyflowline.plot( sVariable_in='flowline_simplified' , sFilename_output_in = 'flowline_simplified_zoom.png', aExtent_in =aExtent_meander ) 
 
     pass
-aCell = oPyflowline.mesh_generation()
+
+if iFlag_simulation == 1:
+    aCell = oPyflowline.mesh_generation()
 
 if iFlag_visualization == 1:
     oPyflowline.plot( sVariable_in='mesh', sFilename_output_in = 'mesh.png' ) 
     pass
 
-oPyflowline.reconstruct_topological_relationship(aCell)
+if iFlag_simulation == 1:
+    oPyflowline.reconstruct_topological_relationship(aCell)
 
 if iFlag_visualization == 1:
     oPyflowline.plot(  sVariable_in='overlap', sFilename_output_in = 'mesh_w_flowline.png',)

paper.bib

@@ -1,5 +1,5 @@
 @article{Engwirda:2021,
-  title={‘UNIFIED’LAGUERRE-POWER MESHES FOR COUPLED EARTH SYSTEM MODELLING},
+  title={'Unified' laguerre-Power Meshes For Coupled Earth System Modelling},
   author={Engwirda, Darren and Liao, Chang},
   year = {2021},
   doi = {10.5281/zenodo.5558988},
@@ -23,11 +23,11 @@ @article{Feng:2022
 }
 
 
-@article{Liao:2022,
+@article{Liao:2023,
 author = {Liao, Chang and Zhou, Tian and Xu, Donghui and Cooper, Matthew G and Engwirda, Darren and Li, Hong-Yi and Leung, L. Ruby},
 title = {Topological relationship-based flow direction modeling: Mesh-independent river networks representation},
 journal = {Journal of Advances in Modeling Earth Systems},
-year = {2022},
+year = {2023},
 volume = {n/a},
 number = {n/a},
 pages = {e2022MS003089},

paper.md

@@ -34,7 +34,7 @@ River networks are crucial in hydrologic and Earth system models. Accurately rep
 For hydrologic modelers, river networks are a key input for hydrologic models. 
 While some hydrologic models accept vector-based river networks [@Schwenk:2021], others only accept mesh cell-based, which requires a generation method from the vector-based river network. 
 Currently, generating a mesh cell-based river network from a given vector-based river network and arbitrary computational mesh is a major challenge.
-Existing methods are typically limited to structured rectangular meshes, such as 30m x 30m cartesian grids for high-resolution watershed-scale modeling or 0.5 degree x 0.5 degree geographic grids for global climate modeling. 
+Existing methods are typically limited to structured rectangular meshes, such as 30m x 30m cartesian grids for high-resolution watershed-scale modeling or 0.5 degree x 0.5 degree geographic grids for global climate modeling. In PyFlowline, we define structured meshes (e.g., lat-lon, raster files with projections, and hexagon) as those with fixed cell sizes and shapes and unstructured meshes as those with variable cell sizes and shapes.
 
 Structured mesh-based methods use fixed cartesian or geographic cell sizes, which have several limitations: (1) they perform poorly at coarse resolution (>1km), and (2) they cannot be seamlessly coupled with other unstructured mesh-based hydrologic models such as oceanic models [@Engwirda:2021]. In contrast, unstructured meshes offer a flexible structure with variable grid-cell sizes and shapes. This flexibility makes them ideal for adapting to complex geometry such as river channels and coastlines. Besides, unstructured meshes provide the flexibility to couple different hydrologic models under a unified framework.
 Thus, unstructured meshes are increasingly being adopted in hydrologic modeling.
@@ -45,14 +45,14 @@ A mesh-independent river network representation method that preserves topologica
 PyFlowline is a Python package that provides a framework for generating river networks for hydrologic models, meeting the identified need. Using an object-oriented programming approach, PyFlowline represents river network elements and mesh cell relationships. It relies on open-source Python libraries like GDAL and Cython for data input/output and spatial data operations.
 
 The computational geometry algorithms used in PyFlowline are designed and implemented using a unified spherical framework, making it suitable for regional and global-scale simulations. PyFlowline uses topological relationships to capture the river networks so they are preserved even at coarse spatial resolutions.
-Moreover, PyFlowline is mesh-independent, supporting both structured and unstructured meshes. It can quickly adopt other mesh types, such as triangulated irregular networks (TIN) or discrete global grid systems (DGGs) [@Sahr:2011]. PyFlowline is a core component of the HexWatershed model, a mesh-independent flow direction model. Several scientific studies focused on coupled Earth system models [@Feng:2022; @Liao:2022] have utilized PyFlowline. A workshop tutorial has also been provided online and in person to support its implementation.
+Moreover, PyFlowline is mesh-independent, supporting both structured and unstructured meshes. It can quickly adopt other mesh types, such as triangulated irregular networks (TIN) or discrete global grid systems (DGGs) [@Sahr:2011]. PyFlowline is a core component of the HexWatershed model, a mesh-independent flow direction model. Several scientific studies focused on coupled Earth system models [@Feng:2022; @Liao:2023] have utilized PyFlowline. A workshop tutorial has also been provided online and in person to support its implementation.
 
 
 # Model features
 
 PyFlowline uses Python's OOP architecture to describe river networks using three essential elements: segments, reaches, and confluences. When applicable, river networks are processed as objects throughout the package.
 
-![The data model. \label{fig:oop}](https://github.com/changliao1025/pyflowline/blob/main/docs/figures/basic_element.png?raw=true)
+![The data model. A vertex class object represents a point on the Earth surface. It have three coordiantes. An edge class object represents a directed line between two points. Besides, it has a length attribute. A flowline class object represents a list of connected lines. \label{fig:oop}](https://github.com/changliao1025/pyflowline/blob/main/docs/figures/basic_element.png?raw=true)
 
 PyFlowline provides several key features, including
 

pyflowline/algorithms/auxiliary/calculate_area_of_difference.py

@@ -3,10 +3,7 @@
 from osgeo import ogr, osr
 import importlib
 
-from shapely.ops import polygonize
-
 from pyflowline.algorithms.auxiliary.find_index_in_list import find_list_in_list 
-
 from pyflowline.external.pyearth.gis.gdal.gdal_functions import calculate_angle_betwen_vertex_normal
 from pyflowline.external.pyearth.gis.gdal.gdal_functions import calculate_polygon_area
 
@@ -16,6 +13,13 @@
 else:
     from pyflowline.algorithms.auxiliary.find_vertex_in_list import find_vertex_in_list
 
+iFlag_shapely = importlib.util.find_spec("shapely") 
+if iFlag_shapely is not None:    
+    from shapely.ops import polygonize
+else:
+    print('shapely is required for this function')
+    pass
+
 
 def calculate_area_of_difference_raw(sFilename_a, sFilename_b):
     #not yet supported

pyflowline/algorithms/intersect/intersect_flowline_with_flowline.py

@@ -1,10 +1,9 @@
 import os
 import numpy as np
 from osgeo import ogr, osr
-from shapely.wkt import loads
+#from shapely.wkt import loads
 from pyflowline.classes.vertex import pyvertex
 
-
 import importlib
 iFlag_cython = importlib.util.find_spec("cython") 
 if iFlag_cython is not None:
@@ -22,7 +21,7 @@ def intersect_flowline_with_flowline( sFilename_flowline_a_in, sFilename_flowlin
     if os.path.exists(sFilename_output_in): 
         os.remove(sFilename_output_in)
 
-    pDriver_geojson = ogr.GetDriverByName( "GeoJSON")
+    pDriver_geojson = ogr.GetDriverByName("GeoJSON")
     pDataset_flowline_a = pDriver_geojson.Open(sFilename_flowline_a_in, 0)
     pDataset_flowline_b = pDriver_geojson.Open(sFilename_flowline_b_in, 0)   
     pLayer_flowline_a = pDataset_flowline_a.GetLayer(0)
@@ -50,9 +49,9 @@ def intersect_flowline_with_flowline( sFilename_flowline_a_in, sFilename_flowlin
 
     pDataset_out = pDriver_geojson.CreateDataSource(sFilename_output_in)
 
-    pLayerOut = pDataset_out.CreateLayer('flowline', pSpatial_reference_b, ogr.wkbMultiPoint)
+    pLayerOut = pDataset_out.CreateLayer('intersect', pSpatial_reference_b, ogr.wkbMultiPoint)
     # Add one attribute
-    pLayerOut.CreateField(ogr.FieldDefn('id', ogr.OFTInteger64)) #long type for high resolution
+    pLayerOut.CreateField(ogr.FieldDefn('pointid', ogr.OFTInteger64)) #long type for high resolution
 
     pLayerDefn = pLayerOut.GetLayerDefn()
     pFeatureOut = ogr.Feature(pLayerDefn)    
@@ -126,7 +125,7 @@ def intersect_flowline_with_flowline( sFilename_flowline_a_in, sFilename_flowlin
                             else:                       
                                 aVertex_intersect.append(pVertex)
                                 pFeatureOut.SetGeometry(point)
-                                pFeatureOut.SetField("id", lVertexID)         
+                                pFeatureOut.SetField("pointid", lVertexID)         
                                 pLayerOut.CreateFeature(pFeatureOut)    
                                 lVertexID = lVertexID + 1
 
@@ -143,7 +142,7 @@ def intersect_flowline_with_flowline( sFilename_flowline_a_in, sFilename_flowlin
                         else:                
                             aVertex_intersect.append(pVertex)
                             pFeatureOut.SetGeometry(pGeometry_intersect)
-                            pFeatureOut.SetField("id", lVertexID)         
+                            pFeatureOut.SetField("pointid", lVertexID)         
                             pLayerOut.CreateFeature(pFeatureOut)    
                             lVertexID = lVertexID + 1