Skip to content

Commit

Permalink
fix several issue related to id, after essd testing
Browse files Browse the repository at this point in the history
  • Loading branch information
@changliao1025
changliao1025 committed Jul 12, 2024
1 parent 165f360 commit b2745cd
Show file tree
Hide file tree
Showing 15 changed files with 427 additions and 346 deletions.
115 changes: 58 additions & 57 deletions pyflowline/algorithms/auxiliary/calculate_area_of_difference.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,32 +3,32 @@
from osgeo import ogr, osr
import importlib.util

from pyflowline.algorithms.auxiliary.find_index_in_list import find_list_in_list
from pyflowline.algorithms.auxiliary.find_index_in_list import find_list_in_list

from pyearth.gis.geometry.calculate_polygon_area import calculate_polygon_area
from pyearth.gis.geometry.calculate_angle_between_vertex_normal import calculate_angle_between_vertex_normal

iFlag_cython = importlib.util.find_spec("cython")
iFlag_cython = importlib.util.find_spec("cython")
if iFlag_cython is not None:
from pyflowline.algorithms.cython.kernel import find_vertex_in_list
else:
from pyflowline.algorithms.auxiliary.find_vertex_in_list import find_vertex_in_list




def calculate_area_of_difference_raw(sFilename_a, sFilename_b):
#not yet supported
return

def calculate_area_of_difference_simplified(aFlowline_in, aVertex_all_in, sFilename_output_in):
iFlag_shapely = importlib.util.find_spec("shapely")
if iFlag_shapely is not None:
def calculate_area_of_difference_simplified(aFlowline_in, aVertex_all_in, sFilename_output_in):
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
if os.path.exists(sFilename_output_in):
if os.path.exists(sFilename_output_in):
os.remove(sFilename_output_in)
pass

Expand All @@ -37,32 +37,33 @@ def calculate_area_of_difference_simplified(aFlowline_in, aVertex_all_in, sFilen
nVeretx = len(aVertex_all_in)
#rebuild index
for i in range(nFlowline):
aFlowline_in[i].lFlowlineID = i
aFlowline_in[i].lFlowlineID = i + 1
for i in range(nVeretx):
aVertex_all_in[i].lVertexID = i
aVertex_all_in[i].lVertexID = i + 1

for i in range(nVeretx):
for i in range(nVeretx):
for j in range(nFlowline):
if aVertex_all_in[i] == aFlowline_in[j].pVertex_start:
aFlowline_in[j].pVertex_start.lVertexID = i
aFlowline_in[j].pVertex_start.lVertexID = i + 1

elif aVertex_all_in[i] == aFlowline_in[j].pVertex_end:
aFlowline_in[j].pVertex_end.lVertexID = i
aFlowline_in[j].pVertex_end.lVertexID = i + 1

for pFlowline in aFlowline_in:
for pFlowline in aFlowline_in:
for j in range(len(aFlowline_in)):
pFlowline2 = aFlowline_in[j]
if pFlowline.lFlowlineID != pFlowline2.lFlowlineID:
if (pFlowline.pVertex_start == pFlowline2.pVertex_start ):
pFlowline.aFlowlineID_start_start.append(pFlowline2.lFlowlineID)
pFlowline.aFlowlineID_start_start.append(pFlowline2.lFlowlineID-1)

if (pFlowline.pVertex_start==pFlowline2.pVertex_end):
pFlowline.aFlowlineID_start_end.append(pFlowline2.lFlowlineID)
pFlowline.aFlowlineID_start_end.append(pFlowline2.lFlowlineID-1)

if (pFlowline.pVertex_end==pFlowline2.pVertex_start ):
pFlowline.aFlowlineID_end_start.append(pFlowline2.lFlowlineID)
pFlowline.aFlowlineID_end_start.append(pFlowline2.lFlowlineID-1)

if (pFlowline.pVertex_end==pFlowline2.pVertex_end):
pFlowline.aFlowlineID_end_end.append(pFlowline2.lFlowlineID)
pFlowline.aFlowlineID_end_end.append(pFlowline2.lFlowlineID-1)

def get_next_branch(iFlag_rightleft, iFlag_reverse, pFlowline_in):
"""
Expand All @@ -75,19 +76,19 @@ def get_next_branch(iFlag_rightleft, iFlag_reverse, pFlowline_in):
"""

nEdge = pFlowline_in.nEdge
iFlag_exist = 0
iFlag_exist = 0
if iFlag_reverse ==0 :
#normal direction
pVertex_stop = pFlowline_in.pVertex_end
pEdge = pFlowline_in.aEdge[nEdge-1]
x1 = pEdge.pVertex_start.dLongitude_degree
y1 = pEdge.pVertex_start.dLatitude_degree
x2 = pEdge.pVertex_end.dLongitude_degree
y2 = pEdge.pVertex_end.dLatitude_degree
y2 = pEdge.pVertex_end.dLatitude_degree
iIndex_right = -1
n = len(pFlowline_in.aFlowlineID_end_start) + len(pFlowline_in.aFlowlineID_end_end)
if n == 0:
iFlag_exist = 0
iFlag_exist = 0
iFlag_reverse_new= 0
pFlowline_out = None
pVertex_stop_out = None
Expand All @@ -105,10 +106,10 @@ def get_next_branch(iFlag_rightleft, iFlag_reverse, pFlowline_in):
pEdge2 = pFlowline_dummy.aEdge[nEdge2-1]
pVertex_dummy = pEdge2.pVertex_start

#calculate angle
#calculate angle
x3 = pVertex_dummy.dLongitude_degree
y3 = pVertex_dummy.dLatitude_degree
angle_dummy = calculate_angle_between_vertex_normal( x1, y1, x2, y2, x3, y3 )
angle_dummy = calculate_angle_between_vertex_normal( x1, y1, x2, y2, x3, y3 )
aAngle.append(angle_dummy)
for i in pFlowline_in.aFlowlineID_end_end:
pFlowline_dummy = aFlowline_in[i]
Expand All @@ -121,17 +122,17 @@ def get_next_branch(iFlag_rightleft, iFlag_reverse, pFlowline_in):
pEdge2 = pFlowline_dummy.aEdge[nEdge2-1]
pVertex_dummy = pEdge2.pVertex_start

#calculate angle
#calculate angle
x3 = pVertex_dummy.dLongitude_degree
y3 = pVertex_dummy.dLatitude_degree
angle_dummy = calculate_angle_between_vertex_normal( x1, y1, x2, y2, x3, y3 )
angle_dummy = calculate_angle_between_vertex_normal( x1, y1, x2, y2, x3, y3 )
aAngle.append(angle_dummy)

#mini
dummy = pFlowline_in.aFlowlineID_end_start + pFlowline_in.aFlowlineID_end_end
if iFlag_rightleft ==0 :
min_angle = min(aAngle)
iIndex_right = aAngle.index(min_angle)
iIndex_right = aAngle.index(min_angle)
pFlowline_out = aFlowline_in[ dummy[iIndex_right] ]
pass
else:
Expand All @@ -140,10 +141,10 @@ def get_next_branch(iFlag_rightleft, iFlag_reverse, pFlowline_in):
pFlowline_out = aFlowline_in[ dummy[iIndex_left]]
pass

if pFlowline_out.pVertex_start == pVertex_stop:
if pFlowline_out.pVertex_start == pVertex_stop:
iFlag_reverse_new=0
pVertex_stop_out = pFlowline_out.pVertex_end
else:
else:
iFlag_reverse_new=1
pVertex_stop_out=pFlowline_out.pVertex_start

Expand All @@ -153,7 +154,7 @@ def get_next_branch(iFlag_rightleft, iFlag_reverse, pFlowline_in):
else:
#rever direction
pVertex_stop = pFlowline_in.pVertex_start
pEdge = pFlowline_in.aEdge[0]
pEdge = pFlowline_in.aEdge[0]
x1 = pEdge.pVertex_end.dLongitude_degree
y1 = pEdge.pVertex_end.dLatitude_degree
x2 = pEdge.pVertex_start.dLongitude_degree
Expand All @@ -162,7 +163,7 @@ def get_next_branch(iFlag_rightleft, iFlag_reverse, pFlowline_in):
iIndex_right = -1
n = len(pFlowline_in.aFlowlineID_start_start) + len(pFlowline_in.aFlowlineID_start_end)
if n == 0:
iFlag_exist = 0
iFlag_exist = 0
iFlag_reverse_new= 0
pFlowline_out = None
pVertex_stop_out = None
Expand Down Expand Up @@ -199,7 +200,7 @@ def get_next_branch(iFlag_rightleft, iFlag_reverse, pFlowline_in):
y3 = pVertex_dummy.dLatitude_degree
angle_dummy = calculate_angle_between_vertex_normal( x1, y1, x2, y2, x3, y3 )
aAngle.append(angle_dummy)

#mini
dummy = pFlowline_in.aFlowlineID_start_start + pFlowline_in.aFlowlineID_start_end
if iFlag_rightleft ==0 :
Expand All @@ -213,16 +214,16 @@ def get_next_branch(iFlag_rightleft, iFlag_reverse, pFlowline_in):
pFlowline_out = aFlowline_in[dummy[iIndex_left]]
pass

if pFlowline_out.pVertex_start == pVertex_stop:
if pFlowline_out.pVertex_start == pVertex_stop:
iFlag_reverse_new=0
pVertex_stop_out = pFlowline_out.pVertex_end
else:
else:
iFlag_reverse_new=1
pVertex_stop_out=pFlowline_out.pVertex_start

iFlag_exist = 1
return iFlag_exist, iFlag_reverse_new, pFlowline_out, pVertex_stop_out


def walk_cycle(iFlag_rightleft, iFlag_reverse, pFlowline_in, pVertex_origin, aFlowline_list, aVertex_list ):
iFlag_loop = 1
Expand All @@ -231,14 +232,14 @@ def walk_cycle(iFlag_rightleft, iFlag_reverse, pFlowline_in, pVertex_origin, aFl
if iFlag_right == 1:
iFlag_loop = 0
return iFlag_loop, aFlowline_list
else:
iFlag_exist, iFlag_reverse_new, pFlowline_next, pVertex_stop = get_next_branch(iFlag_rightleft, 0, pFlowline_in)
else:
iFlag_exist, iFlag_reverse_new, pFlowline_next, pVertex_stop = get_next_branch(iFlag_rightleft, 0, pFlowline_in)
if iFlag_exist ==1:
aFlowline_list.append(pFlowline_next.lFlowlineID)
if (pVertex_stop == pVertex_origin):
aFlowline_list.append(pFlowline_next.lFlowlineID-1)
if (pVertex_stop == pVertex_origin):
iFlag_loop = 1
return iFlag_loop, aFlowline_list

else:
iFlag_vertex_exist, index_dummy = find_vertex_in_list(aVertex_list, pVertex_stop)
if iFlag_vertex_exist ==1:
Expand All @@ -260,14 +261,14 @@ def walk_cycle(iFlag_rightleft, iFlag_reverse, pFlowline_in, pVertex_origin, aFl
if iFlag_right == 1:
iFlag_loop = 0
return iFlag_loop, aFlowline_list
else:
iFlag_exist, iFlag_reverse_new, pFlowline_next, pVertex_stop = get_next_branch(iFlag_rightleft,1, pFlowline_in)
if iFlag_exist ==1:
aFlowline_list.append(pFlowline_next.lFlowlineID)
else:
iFlag_exist, iFlag_reverse_new, pFlowline_next, pVertex_stop = get_next_branch(iFlag_rightleft,1, pFlowline_in)
if iFlag_exist ==1:
aFlowline_list.append(pFlowline_next.lFlowlineID-1)
if (pVertex_stop == pVertex_origin):
#a loop is finished
iFlag_loop = 1
return iFlag_loop, aFlowline_list
return iFlag_loop, aFlowline_list
else:
iFlag_vertex_exist, index_dummy = find_vertex_in_list(aVertex_list, pVertex_stop)
if iFlag_vertex_exist ==1:
Expand All @@ -283,7 +284,7 @@ def walk_cycle(iFlag_rightleft, iFlag_reverse, pFlowline_in, pVertex_origin, aFl
else:
iFlag_loop = 0
return iFlag_loop, aFlowline_list

return iFlag_loop, aFlowline_list

aList_all = list()
Expand All @@ -292,23 +293,23 @@ def walk_cycle(iFlag_rightleft, iFlag_reverse, pFlowline_in, pVertex_origin, aFl
aFlowline_list=list()
aVertex_list=list()
iFlag_loop , aFlowline_list = walk_cycle(0, 0, pFlowline_in, pFlowline_in.pVertex_start, aFlowline_list, aVertex_list)
if iFlag_loop == 1:
if iFlag_loop == 1:
aFlowline_list.insert(0, i)
iFlag = find_list_in_list(aList_all,aFlowline_list )
if iFlag == 1:
pass
else:
aList_all.append(aFlowline_list)
pDataset = pDriver_geojson.CreateDataSource(sFilename_output_in)
pSpatial_reference_gcs = osr.SpatialReference()
aList_all.append(aFlowline_list)

pDataset = pDriver_geojson.CreateDataSource(sFilename_output_in)
pSpatial_reference_gcs = osr.SpatialReference()
pSpatial_reference_gcs.ImportFromEPSG(4326) # WGS84 lat/lon
pSpatial_reference_gcs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
pLayer = pDataset.CreateLayer('aod', pSpatial_reference_gcs, ogr.wkbPolygon)
# Add one attribute
pLayer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger64)) #long type for high resolution
pLayerDefn = pLayer.GetLayerDefn()
pFeature = ogr.Feature(pLayerDefn)
pFeature = ogr.Feature(pLayerDefn)
lCellID = 0
dArea = 0.0
#shapely method, but not complete
Expand All @@ -324,16 +325,16 @@ def walk_cycle(iFlag_rightleft, iFlag_reverse, pFlowline_in, pVertex_origin, aFl
for k in range(nVertex):
aCoords_gcs[k,0] = pFlowline.aVertex[k].dLongitude_degree
aCoords_gcs[k,1] = pFlowline.aVertex[k].dLatitude_degree
aCoords_gcs= tuple(j for j in aCoords_gcs)
aFlowline_polygon.append(aCoords_gcs)

aCoords_gcs= tuple(j for j in aCoords_gcs)
aFlowline_polygon.append(aCoords_gcs)

dummy = polygonize(aFlowline_polygon)
aPolygon_out = list(dummy)
aPolygon_out = list(dummy)
for po in aPolygon_out:
ring = ogr.Geometry(ogr.wkbLinearRing)
aCoords_gcs = po.exterior.coords
aCoords_gcs= np.array(aCoords_gcs)
aCoords_gcs= np.array(aCoords_gcs)
nPoint = (aCoords_gcs.shape)[0]
lons=list()
lats=list()
Expand Down
Loading

0 comments on commit b2745cd

Please sign in to comment.