fix area function bug

changliao1025 · changliao1025 · commit 20146a18eb38 · 2022-06-01T12:09:33.000-07:00
diff --git a/pyflowline/algorithms/auxiliary/gdal_functions.py b/pyflowline/algorithms/auxiliary/gdal_functions.py
@@ -278,20 +278,24 @@ def calculate_distance_based_on_lon_lat(lon1, lat1, lon2, lat2):
     r = 6378137.0
     return c * r
 
-def calculate_polygon_area(lons, lats,  algorithm = 0, radius = 6378137.0):
+def calculate_polygon_area(aLon_in, aLat_in,  algorithm = 0):
     """
     Computes area of spherical polygon, assuming spherical Earth. 
     Returns result in ratio of the sphere's area if the radius is specified. Otherwise, in the units of provided radius.
     lats and lons are in degrees.
     """
     #TODO: take into account geodesy (i.e. convert latitude to authalic sphere, use radius of authalic sphere instead of mean radius of spherical earth)
-    lats = np.deg2rad(lats)
-    lons = np.deg2rad(lons)
-
+    #for i in range(len(aLon_in)):
+    #    aLon_in[i] = aLon_in[i] + 180.0
+    
+    lons = np.deg2rad(aLon_in)
+    lats = np.deg2rad(aLat_in)
     if algorithm==0:
         # Line integral based on Green's Theorem, assumes spherical Earth      
         #close polygon
-        if lats[0]!=lats[-1]:
+        if lats[0]==lats[-1] and lons[0]==lons[-1] :
+            pass
+        else:
             lats = np.append(lats, lats[0])
             lons = np.append(lons, lons[0])
 
@@ -319,10 +323,12 @@ def calculate_polygon_area(lons, lats,  algorithm = 0, radius = 6378137.0):
         area = abs(sum(integrands))/(4*np.pi) # Could be area of inside or outside
 
         area = min(area,1-area)
-        if radius is not None: #return in units of radius
-            return area * 4*np.pi*radius**2
-        else: #return in ratio of sphere total area
-            return area
+        radius= 6378137.0
+       
+        dArea  = area * 4*np.pi*(radius**2)
+        return dArea
+
+
     elif algorithm==2:
         #L'Huilier Theorem, assumes spherical earth
         #see:
diff --git a/pyflowline/algorithms/simplification/remove_returning_flowline.py b/pyflowline/algorithms/simplification/remove_returning_flowline.py
@@ -62,6 +62,7 @@ def retrieve_flowline_intersect_index(iSegment_in, iStream_order_in, pVertex_end
                                 if dummy2.is_overlap(dummy3) ==1:                                    
                                     pVertex_end_current = pVertex_start    
                                     aCell_flowline.append(lCellID)
+                                    #print(lCellID)
                                     iFlag_found = 1                
                                     iFlag_previous_overlap =1
                                     
@@ -72,6 +73,7 @@ def retrieve_flowline_intersect_index(iSegment_in, iStream_order_in, pVertex_end
                                         iFlag_found = 1                                          
                                         if iFlag_previous_overlap ==1: 
                                             aCell_flowline.append(lCellID)              
+                                            #print(lCellID)
                                             iFlag_previous_overlap=0
                                             pass
                                         else:
@@ -81,6 +83,7 @@ def retrieve_flowline_intersect_index(iSegment_in, iStream_order_in, pVertex_end
                                     else:
                                         pVertex_end_current = pVertex_start    
                                         aCell_flowline.append(lCellID)
+                                        #print(lCellID)
                                         iFlag_found = 1                                      
                                         pass
 
diff --git a/pyflowline/classes/basin.py b/pyflowline/classes/basin.py
@@ -13,8 +13,8 @@
 from pyflowline.formats.read_nhdplus_flowline_shapefile import extract_nhdplus_flowline_shapefile_by_attribute
 from pyflowline.formats.read_nhdplus_flowline_shapefile import track_nhdplus_flowline
 from pyflowline.formats.convert_flowline_to_geojson import convert_flowline_to_geojson
-from pyflowline.formats.export_flowline import export_flowline_to_json
-from pyflowline.formats.export_vertex import export_vertex_to_json
+from pyflowline.formats.export_flowline import export_flowline_to_geojson
+from pyflowline.formats.export_vertex import export_vertex_to_geojson
 from pyflowline.algorithms.auxiliary.text_reader_string import text_reader_string
 from pyflowline.algorithms.split.find_flowline_vertex import find_flowline_vertex
 from pyflowline.algorithms.split.find_flowline_confluence import find_flowline_confluence
@@ -238,15 +238,15 @@ def flowline_simplification(self):
             #not used anymore                
             #aThreshold = np.full(2, 300.0, dtype=float)
             #aFlowline_basin_filtered = connect_disconnect_flowline(aFlowline_basin_filtered, aVertex, aThreshold)
-            #sFilename_out = 'flowline_connect.json'
+            #sFilename_out = 'flowline_connect.geojson'
             #sFilename_out = os.path.join(sWorkspace_output_basin, sFilename_out)    
-            #export_flowline_to_json(iFlag_projected, aFlowline_basin_filtered,pSpatial_reference_gcs, sFilename_out)
+            #export_flowline_to_geojson(iFlag_projected, aFlowline_basin_filtered,pSpatial_reference_gcs, sFilename_out)
             pass
         else:
             pass
             
         if self.iFlag_debug ==1:
-            sFilename_out = 'flowline_before_intersect.json'
+            sFilename_out = 'flowline_before_intersect.geojson'
             sFilename_out = os.path.join(sWorkspace_output_basin, sFilename_out)            
             self.export_flowline(aFlowline_basin_filtered, sFilename_out)
         #calculate length
@@ -256,14 +256,14 @@ def flowline_simplification(self):
         #simplification started
         aVertex = find_flowline_vertex(aFlowline_basin_filtered)
         if self.iFlag_debug ==1:
-            sFilename_out = 'flowline_vertex_without_confluence_before_intersect.json'
+            sFilename_out = 'flowline_vertex_without_confluence_before_intersect.geojson'
             sFilename_out = os.path.join(sWorkspace_output_basin, sFilename_out)
-            export_vertex_to_json( aVertex, sFilename_out)
+            export_vertex_to_geojson( aVertex, sFilename_out)
         aFlowline_basin_simplified = split_flowline(aFlowline_basin_filtered, aVertex)
         if self.iFlag_debug ==1:
-            sFilename_out = 'flowline_split_by_point_before_intersect.json'
+            sFilename_out = 'flowline_split_by_point_before_intersect.geojson'
             sFilename_out = os.path.join(sWorkspace_output_basin, sFilename_out)
-            export_flowline_to_json(aFlowline_basin_simplified, sFilename_out)
+            export_flowline_to_geojson(aFlowline_basin_simplified, sFilename_out)
         #ues location to find outlet
         point= dict()   
         point['dLongitude_degree'] = self.dLongitude_outlet_degree
@@ -273,41 +273,41 @@ def flowline_simplification(self):
         pVertex_outlet = aFlowline_basin_simplified[0].pVertex_end
         self.pVertex_outlet = pVertex_outlet
         if self.iFlag_debug ==1:
-            sFilename_out = 'flowline_direction_before_intersect.json'
+            sFilename_out = 'flowline_direction_before_intersect.geojson'
             sFilename_out = os.path.join(sWorkspace_output_basin, sFilename_out)
-            export_flowline_to_json( aFlowline_basin_simplified,  sFilename_out)
+            export_flowline_to_geojson( aFlowline_basin_simplified,  sFilename_out)
         #step 4: remove loops
         aFlowline_basin_simplified = remove_flowline_loop(aFlowline_basin_simplified)    
         if self.iFlag_debug ==1:
-            sFilename_out = 'flowline_loop_before_intersect.json'
+            sFilename_out = 'flowline_loop_before_intersect.geojson'
             sFilename_out = os.path.join(sWorkspace_output_basin, sFilename_out)
-            export_flowline_to_json( aFlowline_basin_simplified, sFilename_out)
+            export_flowline_to_geojson( aFlowline_basin_simplified, sFilename_out)
         #using loop to remove small river, here we use 5 steps
         for i in range(3):
             sStep = "{:02d}".format(i+1)
             aFlowline_basin_simplified = remove_small_river(aFlowline_basin_simplified, self.dThreshold_small_river)
             if self.iFlag_debug ==1:
-                sFilename_out = 'flowline_large_'+ sStep +'_before_intersect.json'
+                sFilename_out = 'flowline_large_'+ sStep +'_before_intersect.geojson'
                 sFilename_out =os.path.join(sWorkspace_output_basin, sFilename_out)
-                export_flowline_to_json( aFlowline_basin_simplified,  sFilename_out)
+                export_flowline_to_geojson( aFlowline_basin_simplified,  sFilename_out)
             aVertex, lIndex_outlet, aIndex_headwater,aIndex_middle, aIndex_confluence, aConnectivity = find_flowline_confluence(aFlowline_basin_simplified,  pVertex_outlet)
             if self.iFlag_debug ==1:
-                sFilename_out = 'flowline_vertex_with_confluence_'+ sStep +'_before_intersect.json'
+                sFilename_out = 'flowline_vertex_with_confluence_'+ sStep +'_before_intersect.geojson'
                 sFilename_out = os.path.join(sWorkspace_output_basin, sFilename_out)
-                export_vertex_to_json( aVertex,  sFilename_out, aAttribute_data=aConnectivity)
+                export_vertex_to_geojson( aVertex,  sFilename_out, aAttribute_data=aConnectivity)
             aFlowline_basin_simplified = merge_flowline( aFlowline_basin_simplified,aVertex, pVertex_outlet, aIndex_headwater,aIndex_middle, aIndex_confluence  )  
             if self.iFlag_debug ==1:
-                sFilename_out = 'flowline_merge_'+ sStep +'_before_intersect.json'
+                sFilename_out = 'flowline_merge_'+ sStep +'_before_intersect.geojson'
                 sFilename_out = os.path.join(sWorkspace_output_basin, sFilename_out)
-                export_flowline_to_json( aFlowline_basin_simplified,  sFilename_out)
+                export_flowline_to_geojson( aFlowline_basin_simplified,  sFilename_out)
             if len(aFlowline_basin_simplified) == 1:
                 break
         
         #the final vertex info
         aVertex, lIndex_outlet, aIndex_headwater,aIndex_middle, aIndex_confluence, aConnectivity = find_flowline_confluence(aFlowline_basin_simplified,  pVertex_outlet)
-        sFilename_out = 'vertex_simplified.json'
+        sFilename_out = 'vertex_simplified.geojson'
         sFilename_out = os.path.join(sWorkspace_output_basin, sFilename_out)
-        export_vertex_to_json( aVertex,  sFilename_out, aAttribute_data=aConnectivity)
+        export_vertex_to_geojson( aVertex,  sFilename_out, aAttribute_data=aConnectivity)
         
         #self.dLength_flowline_simplified = self.calculate_flowline_length(aFlowline_basin_simplified)
         aVertex = np.array(aVertex)
@@ -322,13 +322,13 @@ def flowline_simplification(self):
         if self.iFlag_debug ==1:
             sFilename_out = self.sFilename_flowline_segment_index_before_intersect
             sFilename_out = os.path.join(sWorkspace_output_basin, sFilename_out)
-            export_flowline_to_json(  aFlowline_basin_simplified, sFilename_out, \
+            export_flowline_to_geojson(  aFlowline_basin_simplified, sFilename_out, \
                 aAttribute_data=[aStream_segment], aAttribute_field=['iseg'], aAttribute_dtype=['int'])
         #build stream order 
         aFlowline_basin_simplified, aStream_order = define_stream_order(aFlowline_basin_simplified)
         sFilename_out = self.sFilename_flowline_simplified
         sFilename_out = os.path.join(sWorkspace_output_basin, sFilename_out)
-        export_flowline_to_json(  aFlowline_basin_simplified, sFilename_out, \
+        export_flowline_to_geojson(  aFlowline_basin_simplified, sFilename_out, \
                 aAttribute_data=[aStream_segment, aStream_order], aAttribute_field=['iseg','iord'], aAttribute_dtype=['int','int'])
         
 
@@ -347,6 +347,10 @@ def reconstruct_topological_relationship(self, iMesh_type, sFilename_mesh):
         aCell, aCell_intersect_basin, aFlowline_intersect_all = intersect_flowline_with_mesh(iMesh_type, sFilename_mesh, \
             sFilename_flowline_in, sFilename_flowline_intersect_out)
         sFilename_flowline_filter_geojson = self.sFilename_flowline_filter
+        if self.iFlag_debug ==1:
+            sFilename_out = 'flowline_intersect_flowline_with_mesh.geojson'
+            sFilename_out = os.path.join(sWorkspace_output_basin, sFilename_out)  
+            export_flowline_to_geojson(aFlowline_intersect_all,  sFilename_out)
         
         point= dict()
         point['dLongitude_degree'] = self.dLongitude_outlet_degree
@@ -357,22 +361,38 @@ def reconstruct_topological_relationship(self, iMesh_type, sFilename_mesh):
         #segment based
         aFlowline_basin_conceptual, lCellID_outlet, pVertex_outlet = remove_returning_flowline(iMesh_type, aCell_intersect_basin, pVertex_outlet_initial)
         if self.iFlag_debug ==1:
-            sFilename_out = 'flowline_simplified_after_intersect.json'
+            sFilename_out = 'flowline_simplified_after_intersect.geojson'
             sFilename_out = os.path.join(sWorkspace_output_basin, sFilename_out)  
-            export_flowline_to_json(aFlowline_basin_conceptual,  sFilename_out)
+            export_flowline_to_geojson(aFlowline_basin_conceptual,  sFilename_out)
 
         #edge based
         aFlowline_basin_conceptual, aEdge = split_flowline_to_edge(aFlowline_basin_conceptual)
+        if self.iFlag_debug ==1:
+            sFilename_out = 'flowline_edge_split_flowline_to_edge.geojson'
+            sFilename_out = os.path.join(sWorkspace_output_basin, sFilename_out)
+            export_flowline_to_geojson( aFlowline_basin_conceptual,  sFilename_out)
         aFlowline_basin_conceptual = remove_duplicate_flowline(aFlowline_basin_conceptual)
+        if self.iFlag_debug ==1:
+            sFilename_out = 'flowline_edge_remove_duplicate_flowline.geojson'
+            sFilename_out = os.path.join(sWorkspace_output_basin, sFilename_out)
+            export_flowline_to_geojson( aFlowline_basin_conceptual,  sFilename_out)
         aFlowline_basin_conceptual = correct_flowline_direction(aFlowline_basin_conceptual,  pVertex_outlet )
+        if self.iFlag_debug ==1:
+            sFilename_out = 'flowline_edge_correct_flowline_direction.geojson'
+            sFilename_out = os.path.join(sWorkspace_output_basin, sFilename_out)
+            export_flowline_to_geojson( aFlowline_basin_conceptual,  sFilename_out)
         aFlowline_basin_conceptual = remove_flowline_loop(  aFlowline_basin_conceptual )  
+        if self.iFlag_debug ==1:
+            sFilename_out = 'flowline_edge_remove_flowline_loop.geojson'
+            sFilename_out = os.path.join(sWorkspace_output_basin, sFilename_out)
+            export_flowline_to_geojson( aFlowline_basin_conceptual,  sFilename_out)
   
         aVertex, lIndex_outlet, aIndex_headwater,aIndex_middle, aIndex_confluence, aConnectivity\
             = find_flowline_confluence(aFlowline_basin_conceptual,  pVertex_outlet)
         if self.iFlag_debug ==1:
-            sFilename_out = 'flowline_vertex_with_confluence_after_intersect.json'
+            sFilename_out = 'flowline_vertex_with_confluence_after_intersect.geojson'
             sFilename_out = os.path.join(sWorkspace_output_basin, sFilename_out)
-            export_vertex_to_json( aVertex,  sFilename_out, aAttribute_data=aConnectivity)
+            export_vertex_to_geojson( aVertex,  sFilename_out, aAttribute_data=aConnectivity)
         
         #segment based
         aFlowline_basin_conceptual = merge_flowline( aFlowline_basin_conceptual,aVertex, pVertex_outlet, aIndex_headwater,aIndex_middle, aIndex_confluence  )                          
@@ -393,11 +413,11 @@ def reconstruct_topological_relationship(self, iMesh_type, sFilename_mesh):
         aFlowline_basin_edge, aEdge = split_flowline_to_edge(aFlowline_basin_conceptual)
         sFilename_out = self.sFilename_flowline_edge
         sFilename_out = os.path.join(sWorkspace_output_basin, sFilename_out)
-        export_flowline_to_json(  aFlowline_basin_edge, sFilename_out)
+        export_flowline_to_geojson(  aFlowline_basin_edge, sFilename_out)
 
         sFilename_out = self.sFilename_flowline_conceptual
         sFilename_out = os.path.join(sWorkspace_output_basin, sFilename_out)
-        export_flowline_to_json(  aFlowline_basin_conceptual, sFilename_out, \
+        export_flowline_to_geojson(  aFlowline_basin_conceptual, sFilename_out, \
             aAttribute_data=[aStream_segment, aStream_order], aAttribute_field=['iseg','iord'], aAttribute_dtype=['int','int'])
 
         self.aFlowline_basin_conceptual = aFlowline_basin_conceptual     
@@ -480,7 +500,7 @@ def export(self):
         return
 
     def export_flowline(self, aFlowline_in, sFilename_json_in,iFlag_projected_in = None,  pSpatial_reference_in = None):
-        export_flowline_to_json(aFlowline_in, sFilename_json_in,\
+        export_flowline_to_geojson(aFlowline_in, sFilename_json_in,\
             iFlag_projected_in= iFlag_projected_in, \
             pSpatial_reference_in = pSpatial_reference_in)
 
@@ -580,7 +600,7 @@ def export_config_to_json(self, sFilename_output_in = None):
                 cls=BasinClassEncoder)
         return
 
-    def convert_flowline_to_json(self):
+    def convert_flowline_to_geojson(self):
         sFilename_raw = self.sFilename_flowline_filter            
         sFilename_out = self.sFilename_flowline_filter_geojson
         print('This is the filtered flowline:', sFilename_raw )
diff --git a/pyflowline/classes/pycase.py b/pyflowline/classes/pycase.py
@@ -514,7 +514,7 @@ def analyze(self):
 
     def setup(self):
         if self.iFlag_flowline == 1:
-            self.convert_flowline_to_json()
+            self.convert_flowline_to_geojson()
             pass
         return
 
diff --git a/pyflowline/formats/convert_flowline_to_geojson.py b/pyflowline/formats/convert_flowline_to_geojson.py
@@ -68,7 +68,7 @@ def convert_shapefile_to_geojson(iFlag_type_in, sFilename_shapefile_in, sFilenam
             pSpatial_reference_gcs = osr.SpatialReference()
             pSpatial_reference_gcs.ImportFromEPSG(4326)
             pSpatial_reference_gcs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
-            export_flowline_to_json(aFlowline_basin, sFilename_geojson_in)
+            export_flowline_to_geojson(aFlowline_basin, sFilename_geojson_in)
         else:
             if iFlag_type_in == 2:
                 #polygon
@@ -95,7 +95,7 @@ def convert_shapefile_to_geojson_swat(iFlag_type_in, sFilename_shapefile_in, sFi
             pSpatial_reference_gcs = osr.SpatialReference()
             pSpatial_reference_gcs.ImportFromEPSG(4326)
             pSpatial_reference_gcs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
-            export_flowline_to_json(aFlowline_basin, sFilename_geojson_in)
+            export_flowline_to_geojson(aFlowline_basin, sFilename_geojson_in)
         else:
             if iFlag_type_in == 2:
                 #polygon
diff --git a/pyflowline/formats/export_vertex.py b/pyflowline/formats/export_vertex.py
@@ -2,7 +2,7 @@
 from osgeo import ogr, osr
 from shapely.geometry import Point
 
-def export_vertex_to_json(aVertex_in, \
+def export_vertex_to_geojson(aVertex_in, \
         sFilename_json_in,\
         iFlag_projected_in=None,\
         pSpatial_reference_in=None, \
diff --git a/pyflowline/formats/read_flowline.py b/pyflowline/formats/read_flowline.py
