remove the shapely dependency and update visual

README.md

@@ -31,7 +31,6 @@ PyFlowline depends on the following packages
 1. `numpy`
 2. `gdal`
 3. `netCDF4`
-4. `shapely`
 
 PyFlowline also has three optional dependency packages
 

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 }}

docs/source/installation/installation.rst

@@ -51,8 +51,6 @@ The following dependency packages will be installed during the process.
 * `numpy`
 * `gdal`
 * `netCDF4`
-* `shapely`
-
 
 =============
 Visualization

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                                        
 

pyflowline/algorithms/intersect/intersect_flowline_with_mesh.py

@@ -2,10 +2,11 @@
 import os
 import numpy as np
 from osgeo import ogr, osr
-from shapely.wkt import loads
+#from shapely.wkt import loads
 
 from pyflowline.formats.convert_coordinates import convert_gcs_coordinates_to_cell
 from pyflowline.formats.convert_coordinates import convert_gcs_coordinates_to_flowline
+from pyflowline.external.pyearth.gis.gdal.gdal_functions import get_geometry_coords
 
 def intersect_flowline_with_mesh(iMesh_type_in, sFilename_mesh_in, sFilename_flowline_in, sFilename_output_in):
 
@@ -55,9 +56,10 @@ def intersect_flowline_with_mesh(iMesh_type_in, sFilename_mesh_in, sFilename_flo
     aFlowline_intersect_all=list()   
     for pFeature_mesh in pLayer_mesh:       
         pGeometry_mesh = pFeature_mesh.GetGeometryRef()        
-        dummy0 = loads( pGeometry_mesh.ExportToWkt() )
-        aCoords_gcs = dummy0.exterior.coords
-        aCoords_gcs= np.array(aCoords_gcs)       
+        #dummy0 = loads( pGeometry_mesh.ExportToWkt() )
+        #aCoords_gcs = dummy0.exterior.coords
+        #aCoords_gcs= np.array(aCoords_gcs)       
+        aCoords_gcs = get_geometry_coords(pGeometry_mesh)
 
         lCellID = pFeature_mesh.GetField("cellid")
         dLon = pFeature_mesh.GetField("longitude")

pyflowline/algorithms/intersect/intersect_flowline_with_vertex.py

@@ -1,9 +1,8 @@
 
 
 import os
-import numpy as np
 from osgeo import ogr, osr
-from shapely.wkt import loads
+#from shapely.wkt import loads
 
 def intersect_flowline_with_vertex( sFilename_flowline_in, sFilename_vertex_in, sFilename_output_in):
 
@@ -41,9 +40,9 @@ def intersect_flowline_with_vertex( sFilename_flowline_in, sFilename_vertex_in,
 
     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)    

pyflowline/classes/_visual_basin.py

@@ -32,8 +32,7 @@ def basin_plot(self,
     else:
         sFilename_mesh = sFilename_mesh_in
 
-    sField_thickness = ''
-    iFlag_thickness = 0
+
     if iFlag_type_in ==1:
         #point based
         pass
@@ -190,6 +189,7 @@ def _plot_polygon_variable(self,
                              dData_max_in = None,
                              sFilename_output_in=None,
                              aExtent_in = None,
+                             aLegend_in = None,
                              pProjection_map_in = None):
     """_summary_
 
@@ -275,10 +275,13 @@ def _plot_polygon_variable(self,
         pass
 
     map_vector_polygon_data(sFilename,
+                            iFlag_color_in = 1,
                              iFlag_colorbar_in = 1,
                              sFilename_output_in=sFilename_output_in,
+                             sVariable_in= sVariable,
                              sTitle_in= sTitle,
                              aExtent_in = aExtent_in,
+                             aLegend_in = aLegend_in,
                              pProjection_map_in = pProjection_map_in)
 
 

pyflowline/classes/pycase.py

@@ -540,7 +540,7 @@ def mesh_generation(self, iFlag_antarctic_in=None):
 
                 if iFlag_mesh_boundary ==1:
                     #create a polygon based on real boundary
-                    pBoundary_wkt, pBoundary_shp = read_mesh_boundary(self.sFilename_mesh_boundary)
+                    pBoundary_wkt = read_mesh_boundary(self.sFilename_mesh_boundary)
 
                     aHexagon = create_hexagon_mesh(iFlag_rotation, dX_lowerleft, dY_lowerleft, dResolution_meter, ncolumn, nrow, 
                                                    sFilename_mesh, sFilename_spatial_reference, pBoundary_wkt)
@@ -566,7 +566,7 @@ def mesh_generation(self, iFlag_antarctic_in=None):
                     nrow= int( (dY_upperleft - dY_lowerleft) / dResolution_meter )
                     if iFlag_mesh_boundary ==1:
                         #create a polygon based on real boundary
-                        pBoundary_wkt, pBoundary_shp = read_mesh_boundary(self.sFilename_mesh_boundary)
+                        pBoundary_wkt = read_mesh_boundary(self.sFilename_mesh_boundary)
 
                         aSquare = create_square_mesh(dX_lowerleft, dY_lowerleft, dResolution_meter, ncolumn, nrow, 
                                                      sFilename_mesh, sFilename_spatial_reference, pBoundary_wkt)
@@ -593,7 +593,7 @@ def mesh_generation(self, iFlag_antarctic_in=None):
 
                         if iFlag_mesh_boundary ==1:
                             #create a polygon based on real boundary
-                            pBoundary_wkt, pBoundary_shp = read_mesh_boundary(self.sFilename_mesh_boundary)
+                            pBoundary_wkt = read_mesh_boundary(self.sFilename_mesh_boundary)
                             aLatlon = create_latlon_mesh(dLongitude_left, dLatitude_bot, dResolution_degree, ncolumn, nrow, 
                                                          sFilename_mesh, pBoundary_wkt)
                             pass
@@ -634,7 +634,7 @@ def mesh_generation(self, iFlag_antarctic_in=None):
                                 if iFlag_mesh_boundary ==1:
                                     #create a polygon based on
                                     #read boundary
-                                    pBoundary_wkt, pBoundary_shp = read_mesh_boundary(self.sFilename_mesh_boundary)
+                                    pBoundary_wkt = read_mesh_boundary(self.sFilename_mesh_boundary)
 
                                     aMpas = create_mpas_mesh(iFlag_global, iFlag_use_mesh_dem, iFlag_save_mesh,
                                                              sFilename_mesh_netcdf,  sFilename_mesh, iFlag_antarctic_in=iFlag_antarctic_in, pBoundary_in = pBoundary_wkt)

pyflowline/external/pyearth/gis/gdal/gdal_functions.py

@@ -3,11 +3,9 @@
 import numpy as np
 import osgeo
 from osgeo import ogr, osr, gdal
-from shapely.wkt import loads
+#from shapely.wkt import loads
 #most of these functions are copied from the pyearth package
 
-
-
 import importlib
 iFlag_cython = importlib.util.find_spec("cython") 
 if iFlag_cython is not None:
@@ -71,8 +69,6 @@ def calculate_angle_betwen_vertex(x1, y1, x2, y2, x3, y3):
     angle3deg = angle_between_vectors_degrees(a3vec, c3vec)
     return  angle3deg
 
-
-
 def angle_between_vectors_degrees(u, v):
     """Return the angle between two vectors in any dimension space,
     in degrees.
@@ -469,8 +465,35 @@ def read_mesh_boundary(sFilename_boundary_in):
 
 
     pBoundary_wkt = pBoundary_ogr.ExportToWkt()
-    pBoundary_shp = loads( pBoundary_ogr.ExportToWkt() )
-    return pBoundary_wkt, pBoundary_shp
+
+    return pBoundary_wkt
 
 
-
+def get_geometry_coords(geometry):
+    if geometry.GetGeometryType() == ogr.wkbPolygon:
+        return get_polygon_exterior_coords(geometry)
+    elif geometry.GetGeometryType() == ogr.wkbLineString:
+        return get_linestring_coords(geometry)
+    elif geometry.GetGeometryType() == ogr.wkbPoint:
+        return get_point_coords(geometry)
+    else:
+        raise ValueError("Unsupported geometry type.")
+
+def get_polygon_exterior_coords(polygon_geometry):
+    exterior_coords = []
+    ring = polygon_geometry.GetGeometryRef(0)  # Get the exterior ring
+    for i in range(ring.GetPointCount()):
+        point = ring.GetPoint(i)
+        exterior_coords.append((point[0], point[1]))
+    return np.array(exterior_coords)
+
+def get_linestring_coords(linestring_geometry):
+    coords = []
+    for i in range(linestring_geometry.GetPointCount()):
+        point = linestring_geometry.GetPoint(i)
+        coords.append((point[0], point[1]))
+    return np.array(coords)
+
+def get_point_coords(point_geometry):
+    point = point_geometry.GetPoint()
+    return np.array([(point[0], point[1])])

pyflowline/external/pyearth/visual/animate/animate_vector_polygon_data.py

@@ -1,22 +1,16 @@
 import os
 import json
 import numpy as np
-
 import cartopy.crs as ccrs
-import cartopy.mpl.ticker as ticker
 import matplotlib as mpl
-from shapely.wkt import loads
+#from shapely.wkt import loads
 from osgeo import osr, gdal, ogr
-
-
 import matplotlib.pyplot as plt
-import matplotlib.path as mpath
-import matplotlib.ticker as mticker
 import matplotlib.patches as mpatches
 import matplotlib.cm as cm
 from matplotlib import animation
-from pyearth.toolbox.data.cgpercentiles import cgpercentiles
 from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
+from pyflowline.external.pyearth.gis.gdal.gdal_functions import get_geometry_coords
 
 pProjection = ccrs.PlateCarree()  # for latlon data only
 
@@ -83,9 +77,10 @@ def animate_vector_polygon_data(
         sGeometry_type = pGeometry_in.GetGeometryName()
 
         if sGeometry_type == "POLYGON":
-            dummy0 = loads(pGeometry_in.ExportToWkt())
-            aCoords_gcs = dummy0.exterior.coords
-            aCoords_gcs = np.array(aCoords_gcs)
+            #dummy0 = loads(pGeometry_in.ExportToWkt())
+            #aCoords_gcs = dummy0.exterior.coords
+            #aCoords_gcs = np.array(aCoords_gcs)
+            aCoords_gcs = get_geometry_coords(pGeometry_in)
 
             dLon_max = np.max([dLon_max, np.max(aCoords_gcs[:, 0])])
             dLon_min = np.min([dLon_min, np.min(aCoords_gcs[:, 0])])