Is your feature request related to a problem? Please describe.
There is no direct way to identify hfb cellid1 and cellid2 from polylines and the model grid
Describe the solution you'd like
Extension of intersect capabilities to allow
Describe alternatives you've considered
Custom one-off coding.
Additional context
The code below was developed by @askar-INTERA to determine hfb cellid1 and cellid2 using polylines and the model grid. The code relies on pandas, geopandas, and shapely for the heavy lifting.
# hfb for faults
from shapely.ops import linemerge, unary_union
from shapely import offset_curve
from flopy.utils import GridIntersect
inset_grid=gwf.modelgrid
inset_x=inset_grid.xcellcenters
inset_y=inset_grid.ycellcenters
inset_gdf=inset_grid.geo_dataframe
inset_verts=inset_grid.verts
ix = GridIntersect(inset_grid)
hfb=[]
all_gemos=[]
for idx, shp in enumerate (['regional_faults_usgs_2024.shp', 'Roca_Faults_Update_2010_utm.shp']):
ploygon_ws = os.path.join('Shapefiles',shp)
faults=gpd.read_file(os.path.join(ploygon_ws))
if idx==0:
names=[x.split(" fault")[0].replace(" ","_").lower() for x in faults.Label.values]
else:
names=["rh_"+str(x+1) for x in range(len(faults))]
faults.csr = None
x, y = 28822.363114339616, 4002417.863173714839
faults = faults.translate(xoff=-x, yoff=-y)
faults = faults.scale(xfact=1 / 0.3048, yfact=1 / 0.3048, zfact=1.0, origin=(0, 0))
faults = faults.rotate(43.75, origin=(0, 0))
for poly,nam in zip(faults.geometry.values,names):
result = ix.intersect(poly)
df=pd.DataFrame(columns=["cell1","cell2"])
count=0
for cell in result:
cell_coord=(inset_x[cell[0]],inset_y[cell[0]])
adj_cells=inset_grid.neighbors(cell[0])
for adj_cell in adj_cells:
test_line=LineString([(inset_x[adj_cell],inset_y[adj_cell]),cell_coord])
if test_line.intersects(poly):
df.loc[count]=[cell[0],adj_cell]
count+=1
temp_df = pd.DataFrame(np.sort(df[['cell1', 'cell2']].values, axis=1), columns=['temp1', 'temp2'])
temp_df['original_index'] = df.index
unique_pairs = temp_df.drop_duplicates(subset=['temp1', 'temp2'], keep='first')
df = df.loc[unique_pairs['original_index']]
test_lines=[]
for cell1,cell2 in zip(df.cell1,df.cell2):
cell1_poly=inset_gdf.loc[cell1,"geometry"]
cell2_poly=inset_gdf.loc[cell2,"geometry"]
test_lines.append(cell1_poly.boundary.intersection(cell2_poly.boundary))
all_lines = unary_union(test_lines)
merged_lines = linemerge(all_lines)
assert merged_lines.geom_type=='LineString'
df.loc[:,"name"]=nam
hfb.append(df)
all_gemos.append(merged_lines)
hfb = pd.concat(hfb, ignore_index=True)
temp_df = pd.DataFrame(np.sort(hfb[['cell1', 'cell2']].values, axis=1), columns=['temp1', 'temp2'])
temp_df['original_index'] = hfb.index
unique_pairs = temp_df.drop_duplicates(subset=['temp1', 'temp2'], keep='first')
hfb = hfb.loc[unique_pairs['original_index']]
hfb_lay=[]
for lay in range(inset_nlay):
df=hfb.copy()
df.loc[:,"lay1"]=lay
df.loc[:,"lay2"]=lay
df.loc[:,"hydchr"]=-0.1
df.loc[:,"idom1"]=[new_idom[l, cell] for l,cell in zip(df.lay1.values,df.cell1.values)]
df.loc[:,"idom2"]=[new_idom[l, cell] for l,cell in zip(df.lay2.values,df.cell2.values)]
df = df[df['idom1'] > 0]
df = df[df['idom2'] > 0]
hfb_lay.append(df[["lay1","cell1","lay2","cell2","hydchr"]])
hfb_lay = pd.concat(hfb_lay, ignore_index=True)
hfb = flopy.mf6.ModflowGwfhfb(gwf,
print_input=False,
pname="hfb",
filename="flow.hfb",
stress_period_data =hfb_lay.to_records(index=False),)
gwf.hfb.write()
Is your feature request related to a problem? Please describe.
There is no direct way to identify hfb cellid1 and cellid2 from polylines and the model grid
Describe the solution you'd like
Extension of intersect capabilities to allow
Describe alternatives you've considered
Custom one-off coding.
Additional context
The code below was developed by @askar-INTERA to determine hfb cellid1 and cellid2 using polylines and the model grid. The code relies on pandas, geopandas, and shapely for the heavy lifting.