diff --git a/pyflowline/formats/export_flowline.py b/pyflowline/formats/export_flowline.py index b176b70..f73f562 100644 --- a/pyflowline/formats/export_flowline.py +++ b/pyflowline/formats/export_flowline.py @@ -179,17 +179,21 @@ def export_flowline_to_shapefile(iFlag_projected_in, aFlowline_in, pSpatial_refe for i in range(nFlowline): pFlowline = aFlowline_in[i] dummy =pFlowline.aVertex - aPoint=list() + #aPoint=list() + pLine = ogr.Geometry(ogr.wkbLineString) for j in dummy: if iFlag_projected_in ==1: - aPoint.append( Point( j.dx, j.dy ) ) + #aPoint.append( Point( j.dx, j.dy ) ) + pLine.AddPoint(j.dx, j.dy) pass else: - aPoint.append( Point( j.dLongitude_degree, j.dLatitude_degree ) ) + #aPoint.append( Point( j.dLongitude_degree, j.dLatitude_degree ) ) + pLine.AddPoint(j.dLongitude_degree, j.dLatitude_degree) pass - dummy1= LineString( aPoint ) - pGeometry_out = ogr.CreateGeometryFromWkb(dummy1.wkb) + #dummy1= LineString( aPoint ) + pGeometry_out = ogr.CreateGeometryFromWkb(pLine.ExportToWkb()) + #pGeometry_out = ogr.CreateGeometryFromWkb(dummy1.wkb) pFeature_out.SetGeometry(pGeometry_out) pFeature_out.SetField("id", lID) diff --git a/pyflowline/formats/export_vertex.py b/pyflowline/formats/export_vertex.py index b00aab8..4bb7099 100644 --- a/pyflowline/formats/export_vertex.py +++ b/pyflowline/formats/export_vertex.py @@ -117,14 +117,18 @@ def export_vertex_to_shapefile(aVertex_in, sFilename_shapefile_out,\ lID = 0 for i in range(nVertex): pVertex = aVertex_in[i] + pPoint = ogr.Geometry(ogr.wkbPoint) if iFlag_projected_in ==1: - dummy1= Point( pVertex.dx, pVertex.dy ) + #dummy1= Point( pVertex.dx, pVertex.dy ) + pPoint.AddPoint(pVertex.dx, pVertex.dy) pass else: - dummy1= Point( pVertex.dLongitude_degree, pVertex.dLatitude_degree ) + #dummy1= Point( pVertex.dLongitude_degree, pVertex.dLatitude_degree ) + pPoint.AddPoint(pVertex.dLongitude_degree, pVertex.dLatitude_degree) pass - pGeometry_out = ogr.CreateGeometryFromWkb(dummy1.wkb) + #pGeometry_out = ogr.CreateGeometryFromWkb(dummy1.wkb) + pGeometry_out = ogr.CreateGeometryFromWkb(pPoint.ExportToWkb()) pFeature_out.SetGeometry(pGeometry_out) pFeature_out.SetField("id", lID) if iFlag_attribute ==1: diff --git a/pyflowline/mesh/hexagon/create_hexagon_mesh.py b/pyflowline/mesh/hexagon/create_hexagon_mesh.py index 7fff0df..8f561d3 100644 --- a/pyflowline/mesh/hexagon/create_hexagon_mesh.py +++ b/pyflowline/mesh/hexagon/create_hexagon_mesh.py @@ -438,6 +438,7 @@ def add_cell_into_list2(aList, lCellID, iRow, iColumn, dLongitude_center, dLatit aHexagon_out = list() if iFlag_fill_hole ==1: + #find the virtual neighbors for pCell in aHexagon: aNeighbor_land = pCell.aNeighbor_land #including both holes and maps land cutoff by boundary @@ -449,9 +450,11 @@ def add_cell_into_list2(aList, lCellID, iRow, iColumn, dLongitude_center, dLatit if lNeighbor in aHexagon_dict: nNeighbor_land_update = nNeighbor_land_update + 1 aNeighbor_land_update.append(lNeighbor) + pass else: #a hole or boundary mpas land cell aNeighbor_land_virtual.append(lNeighbor) + pass pCell.nNeighbor= len(aNeighbor_land_update) pCell.aNeighbor = aNeighbor_land_update @@ -460,11 +463,10 @@ def add_cell_into_list2(aList, lCellID, iRow, iColumn, dLongitude_center, dLatit pCell.aNeighbor_land_virtual = aNeighbor_land_virtual pCell.nNeighbor_land_virtual = len(aNeighbor_land_virtual) - aHexagon_middle.append(pCell) + #add holes back, we also need to add into the unsorted map for pCell in aHexagon: - if pCell.nNeighbor_land_virtual ==1: #only one virtual land means it is likely next to a hole lNeighbor_hole = pCell.aNeighbor_land_virtual[0] #now find its row and column indices @@ -528,7 +530,7 @@ def add_cell_into_list2(aList, lCellID, iRow, iColumn, dLongitude_center, dLatit dLatitude_center = np.mean(aCoords[0:6,1]) if lCellID not in aHexagon_dict: - aHexagon_middle, dArea = add_cell_into_list1(aHexagon_middle, lCellID, iRow, iColumn, dLongitude_center,dLatitude_center, aCoords) + aHexagon, dArea = add_cell_into_list1(aHexagon, lCellID, iRow, iColumn, dLongitude_center,dLatitude_center, aCoords) aHexagon_dict[lCellID] = lCellIndex lCellIndex = lCellIndex + 1 pFeature.SetGeometry(pPolygon) @@ -598,7 +600,7 @@ def add_cell_into_list2(aList, lCellID, iRow, iColumn, dLongitude_center, dLatit dLatitude_center = np.mean(aCoords[0:6,1]) if lCellID not in aHexagon_dict: - aHexagon_middle, dArea = add_cell_into_list1(aHexagon_middle, lCellID, iRow, iColumn, dLongitude_center,dLatitude_center, aCoords) + aHexagon, dArea = add_cell_into_list1(aHexagon, lCellID, iRow, iColumn, dLongitude_center,dLatitude_center, aCoords) aHexagon_dict[lCellID] = lCellIndex lCellIndex = lCellIndex + 1 pFeature.SetGeometry(pPolygon) diff --git a/pyflowline/mesh/mpas/create_mpas_mesh.py b/pyflowline/mesh/mpas/create_mpas_mesh.py index a39e6ab..0269d95 100644 --- a/pyflowline/mesh/mpas/create_mpas_mesh.py +++ b/pyflowline/mesh/mpas/create_mpas_mesh.py @@ -421,7 +421,8 @@ def add_cell_into_list(aList, i, lCellID, dArea, dElevation_mean, dElevation_pro aMpas_out = list() ncell = len(aMpas) #generate the list of cell ID that are already certain - if iFlag_fill_hole == 1: + if iFlag_fill_hole == 1: + #first update neighbor information because some cell should have vitual land neighbor (not present in the mesh) #this operation does not increase the number of cells, but it update the neighbor information #specifically, it divided the land neighbor into two parts: land and virtual land @@ -485,7 +486,7 @@ def add_cell_into_list(aList, i, lCellID, dArea, dElevation_mean, dElevation_pro dArea = float(aCellArea[j]) if lCellID not in aMpas_dict: - aMpas_middle = add_cell_into_list(aMpas_middle, j, lCellID, dArea, dElevation_mean, dElevation_profile0, aCoords ) + aMpas = add_cell_into_list(aMpas, j, lCellID, dArea, dElevation_mean, dElevation_profile0, aCoords ) aMpas_dict[lCellID] = lCellIndex lCellIndex = lCellIndex + 1 #now we need to update the neightboring information as well