Skip to content

Commit

Permalink
upgrade the rtree method
Browse files Browse the repository at this point in the history
  • Loading branch information
@changliao1025
changliao1025 committed Jul 21, 2024
1 parent b2745cd commit d6f582d
Show file tree
Hide file tree
Showing 6 changed files with 78 additions and 81 deletions.
4 changes: 0 additions & 4 deletions pyflowline/algorithms/auxiliary/find_index_in_list.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,10 +56,6 @@ def find_vertex_on_edge(aVertex_in, pEdge_in):
else:
if diff < 1.0:
iFlag_overlap = pEdge_in.check_vertex_on_edge(pVertex)

pass


else:
for i in np.arange( nVertex):
pVertex = aVertex_in[i]
Expand Down
19 changes: 9 additions & 10 deletions pyflowline/algorithms/auxiliary/find_vertex_in_list.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,20 +40,20 @@ def find_vertex_in_list(aVertex_in, pVertex_in, dThreshold_in=1.0E-6):
pass

#now the new vertex
x=pVertex_in.dLongitude_degree
y=pVertex_in.dLatitude_degree
left = x - 1E-5
right = x + 1E-5
bottom = y-1E-5
top = y+1E-5
pBound= (left, bottom, right, top)
aIntersect = list(index_vertex.search(pBound))
#x=pVertex_in.dLongitude_degree
#y=pVertex_in.dLatitude_degree
#left = x - 1E-5
#right = x + 1E-5
#bottom = y-1E-5
#top = y+1E-5
#pBound= (left, bottom, right, top)
#aIntersect = list(index_vertex.search(pBound))
aIntersect = list(index_vertex.search_surrounding([pVertex_in.dLongitude_degree, pVertex_in.dLatitude_degree]))

for k in aIntersect:
pVertex = aVertex_in[k]
#dDistance = pVertex.calculate_distance(pVertex_in)
if pVertex == pVertex_in: #if dDistance < dThreshold_in:

iFlag_exist = 1
lIndex = k
break
Expand All @@ -67,7 +67,6 @@ def find_vertex_in_list(aVertex_in, pVertex_in, dThreshold_in=1.0E-6):
pVertex = aVertex_in[i]
#dDistance = pVertex.calculate_distance(pVertex_in)
if pVertex == pVertex_in: #if dDistance < dThreshold_in:

iFlag_exist = 1
lIndex = i
break
Expand Down
17 changes: 9 additions & 8 deletions pyflowline/algorithms/cython/kernel.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,14 +89,15 @@ cpdef find_vertex_in_list(list aVertex_in, pVertex_in, double dThreshold_in = 1
index_vertex.insert(i, pBound)

#now the new vertex
x=pVertex_in.dLongitude_degree
y=pVertex_in.dLatitude_degree
left = x - 1E-5
right = x + 1E-5
bottom = y - 1E-5
top = y + 1E-5
pBound= (left, bottom, right, top)
aIntersect = list(index_vertex.search(pBound))
#x=pVertex_in.dLongitude_degree
#y=pVertex_in.dLatitude_degree
#left = x - 1E-5
#right = x + 1E-5
#bottom = y - 1E-5
#top = y + 1E-5
#pBound= (left, bottom, right, top)
#aIntersect = list(index_vertex.search(pBound))
aIntersect = list(index_vertex.search_surrounding([pVertex_in.dLongitude_degree, pVertex_in.dLatitude_degree]))

for k in aIntersect:
pVertex = aVertex_in[k]
Expand Down
5 changes: 3 additions & 2 deletions pyflowline/algorithms/index/define_stream_order.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,8 +115,9 @@ def define_stream_order(aFlowline_in, aConfluence_in):

#update confluence
x, y = pFlowline_downstream.pVertex_end.dLongitude_degree, pFlowline_downstream.pVertex_end.dLatitude_degree
pBound = (x - 1E-5, y - 1E-5, x + 1E-5, y + 1E-5)
aIntersect = list(index_confluence.search(pBound))
#pBound = (x - 1E-5, y - 1E-5, x + 1E-5, y + 1E-5)
#aIntersect = list(index_confluence.search(pBound))
aIntersect = list(index_confluence.search_surrounding([x, y]))

for k in aIntersect:
pConfluence2 = aConfluence_in[k]
Expand Down
110 changes: 55 additions & 55 deletions pyflowline/algorithms/simplification/remove_returning_flowline.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,44 +2,44 @@
import numpy as np
from pyflowline.formats.convert_coordinates import convert_gcs_coordinates_to_flowline

def remove_returning_flowline(iMesh_type_in, aCell_intersect_in, pVertex_outlet_in):
aFlowline_out=list()
def remove_returning_flowline(iMesh_type_in, aCell_intersect_in, pVertex_outlet_in):
aFlowline_out=list()
nCell = len(aCell_intersect_in)
#lCellID_to_cell = {cell.lCellID: cell for cell in aCell_intersect_in}
def simplify_list(aCell_flowline_in):
aCell_flowline_out = copy.deepcopy(aCell_flowline_in)
nCell2 = len(aCell_flowline_out)
nCell2 = len(aCell_flowline_out)
#iFlag_unique = check_if_duplicates(aCell_flowline_out)
if int(len(aCell_flowline_out) == len(set(aCell_flowline_out))) :
return aCell_flowline_out
else:
for i in range(nCell2):
elem=aCell_flowline_out[i]
if aCell_flowline_out.count(elem) > 1:
dummy = [i for i, x in enumerate(aCell_flowline_out) if x == elem]
dummy = [i for i, x in enumerate(aCell_flowline_out) if x == elem]
start = dummy[0]
end = dummy[-1]
del aCell_flowline_out[start: end]
break
del aCell_flowline_out[start: end]
break

aCell_flowline_out = simplify_list(aCell_flowline_out)
return aCell_flowline_out

def retrieve_flowline_intersect_index(iSegment_in, iStream_order_in, pVertex_end_in):
lCellID =-1
lCellID =-1
iFlag_found = 1
pVertex_end_current = pVertex_end_in
aCell_flowline=list()
aCell_flowline=list()
aSegment_upstream = list()
aVertex_end_upstream = list()
aStream_order = list()
iFlag_skip = 0
iFlag_previous_overlap=0
while iFlag_found == 1:
iFlag_found = 0
iFlag_found = 0
for j in range(nCell):
pCell = aCell_intersect_in[j]
lCellID = pCell.lCellID
lCellID = pCell.lCellID
aFlowline= pCell.aFlowline
nFlowline = len(aFlowline)
for i in range(nFlowline):
Expand All @@ -49,49 +49,49 @@ def retrieve_flowline_intersect_index(iSegment_in, iStream_order_in, pVertex_end
iStream_segment = pFlowline.iStream_segment
iStream_order = pFlowline.iStream_order
if pVertex_end == pVertex_end_current: #be careful because this function is strict
if iStream_segment == iSegment_in:
if iStream_segment == iSegment_in:
dLength = pFlowline.dLength
iFlag_found2, dummy2 = pCell.which_edge_cross_this_vertex(pVertex_start)
iFlag_found3, dummy3 = pCell.which_edge_cross_this_vertex(pVertex_end)
if iFlag_found2 == 1 and iFlag_found3 == 1: #on the edge
#same edge
if dummy2.is_overlap(dummy3) ==1:
pVertex_end_current = pVertex_start
aCell_flowline.append(lCellID)
iFlag_found = 1
if dummy2.is_overlap(dummy3) ==1:
pVertex_end_current = pVertex_start
aCell_flowline.append(lCellID)
iFlag_found = 1
iFlag_previous_overlap =1

pass
else:
if dLength < (0.1*pCell.dLength) : #because it taks a short cut
pVertex_end_current = pVertex_start
iFlag_found = 1
if iFlag_previous_overlap ==1:
aCell_flowline.append(lCellID)
if dLength < (0.1*pCell.dLength) : #because it taks a short cut
pVertex_end_current = pVertex_start
iFlag_found = 1
if iFlag_previous_overlap ==1:
aCell_flowline.append(lCellID)
iFlag_previous_overlap=0
pass
else:
iFlag_previous_overlap=1
pass
pass
pass
else:
pVertex_end_current = pVertex_start
pVertex_end_current = pVertex_start
aCell_flowline.append(lCellID)
iFlag_found = 1
iFlag_found = 1


else:
pVertex_end_current = pVertex_start
aCell_flowline.append(lCellID)
pVertex_end_current = pVertex_start
aCell_flowline.append(lCellID)
iFlag_found = 1
iFlag_previous_overlap=0
iFlag_previous_overlap=0

break

else:
iFlag_found2, dummy2 = pCell.which_edge_cross_this_vertex(pVertex_start)
iFlag_found3, dummy3 = pCell.which_edge_cross_this_vertex(pVertex_end)
if pVertex_end == pVertex_end_in:
if pVertex_end == pVertex_end_in:
iFlag_found = 1
pass
else:
Expand All @@ -100,64 +100,64 @@ def retrieve_flowline_intersect_index(iSegment_in, iStream_order_in, pVertex_end
#save the upstream information for new step
aSegment_upstream.append(iStream_segment)
aVertex_end_upstream.append(pVertex_end)
aStream_order.append(iStream_order)
aStream_order.append(iStream_order)

pass
pass

if iFlag_found == 1:
break

#should have finished search
if iFlag_found == 0:
#reverse
aCell_flowline = aCell_flowline[::-1]
#simplify list
#reverse
aCell_flowline = aCell_flowline[::-1]
#simplify list
aCell_simple = simplify_list(aCell_flowline)
#save the output
nCell3 = len(aCell_simple)
aCoordinates = list()
aCoordinates = list()
if nCell3 >1:
for i in range(nCell3):
lCellID = aCell_simple[i]
for j in range( nCell):
if aCell_intersect_in[j].lCellID == lCellID:
if aCell_intersect_in[j].lCellID == lCellID:
x = aCell_intersect_in[j].dLongitude_center_degree
y = aCell_intersect_in[j].dLatitude_center_degree
y = aCell_intersect_in[j].dLatitude_center_degree
aCoordinates.append([x,y])

pFlowline = convert_gcs_coordinates_to_flowline(aCoordinates)
pFlowline.iStream_segment = iSegment_in
pFlowline.iStream_order=iStream_order_in
aFlowline_out.append(pFlowline)

sort_index = np.argsort(aStream_order)
sort_index = np.argsort(aStream_order)
#sort_index = sort_index[::-1] #can we process low order first?
#nUpstream = len(aSegment_upstream)
for i in sort_index:
retrieve_flowline_intersect_index(aSegment_upstream[i], aStream_order[i], aVertex_end_upstream[i])
pass
pass

return
dDiatance_min = float('inf')
lCellID_outlet = lCellID_outlet2 = lCellID_outlet3 = None
return

dDistance_min = float('inf')
lCellID_outlet = lCellIndex_outlet = lFlowlineIndex_outlet = None
for j, pCell in enumerate(aCell_intersect_in):
for i, pFlowline in enumerate(pCell.aFlowline):
dDiatance = pFlowline.pVertex_end.calculate_distance(pVertex_outlet_in)
if dDiatance < dDiatance_min:
dDiatance_min = dDiatance
if dDiatance < dDistance_min:
dDistance_min = dDiatance
lCellID_outlet = pCell.lCellID
lCellID_outlet2 = j
lCellID_outlet3 = i
lCellIndex_outlet = j
lFlowlineIndex_outlet = i

#loop through
pFlowline = aCell_intersect_in[lCellID_outlet2].aFlowline[lCellID_outlet3]
#loop through
pFlowline = aCell_intersect_in[lCellIndex_outlet].aFlowline[lFlowlineIndex_outlet]
iStream_segment =pFlowline.iStream_segment
iStream_order=pFlowline.iStream_order
pVertex_outlet=pFlowline.pVertex_end

retrieve_flowline_intersect_index(iStream_segment, iStream_order, pVertex_outlet)

#update outlet location since the conceptual flowline use center
Expand Down
4 changes: 2 additions & 2 deletions pyflowline/classes/pycase.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1168,8 +1168,8 @@ def pyflowline_run(self):
#may be read mesh
iMesh_type = self.iMesh_type
#there must be some auxiliary file associated with the mesh file
self.aCell = read_mesh_json_w_topology(iMesh_type, self.sFilename_mesh)
aCell_out = self.aCell
#self.aCell = read_mesh_json_w_topology(iMesh_type, self.sFilename_mesh)
#aCell_out = self.aCell
pass

if self.iFlag_intersect == 1:
Expand Down

0 comments on commit d6f582d

Please sign in to comment.