Does sfrmaker work with mf6 disv grids?

[porterma]

I have been working my way through this package the past week or so and things look promising. I am trying to intersect my stream shape file with an unstructured grid that is created with gridgen.exe. I am using a quadtree mesh. I am stuck at defining the grid that goes into stream_lines.to_sfr(grid, model=mdl) . It works fine with the StructuredGrid.from_sr package when I am using a structured grid. However, when I try UnstructuredGrid.from_dataframe(), I get an error that says AssertionError: structured=False can only be specified for mfusg models. Does sfrmaker work with gridgen quadtree mesh grids and if so do you have any examples? Thanks!

[aleaf]

Hi @porterma, thanks for raising this issue. In principal, SFRmaker should work with any kind of unstructured grid, but I haven't had any project work so far to motivate implementing that capability. It would be great to add a quadtree example/test case. How is your grid refined in the z-direction? Does it coarsen with depth away from the streams? The intersection and setup of the SFR reaches should be fine, but one potential hang-up I can imagine is the streams cutting into layers below layer 1, after the streambed elevations are smoothed. In this case, either the layer bottoms need to be adjusted so that the streambed bottom is above the bottom of layer 1, or the reach needs to be assigned to the appropriate layer. In the later case, if the grid coarsens downward, those reaches would be in unrefined cells, and there would be a bunch of "wasted" refined cells above them.

[porterma]

@aleaf, Thanks for the quick reply. My grid uses the same 2D quadtree mesh for each layer, so it does not coarsen with depth. Additionally, all of my streams are in the first layer. I am reading the gridgen.exe quadtree shapefile using geopandas, then using Unstructured grid to create a grid. Code is

stream_lines = lines.from_dataframe(streams_df, attr_length_units='feet', attr_height_units='feet', epsg='2243')
grid_shp = 'qtgrid.shp'
grd_df = geopandas.read_file(grid_shp)
grd = UnstructuredGrid.from_dataframe(grd_df, node_col='nodenumber', geometry_column='geometry', model_units='feet', proj_str='epsg:2243', )
sfr = stream_lines.to_sfr(grd)

sfrmaker prints out the following to the screen:

SFRmaker version 0.1.2+3.g17dfa49

Creating sfr dataset...
Model grid information
unstructured grid
nnodes: 47,871
model length units: feet
crs: epsg:2243
bounds: 2709517.94, 1150502.67, 2772838.51, 1217702.18
active area defined by: all cells

Culling hydrography to active area...
starting lines: 99
remaining lines: 99
finished in 0.03s

Intersecting 99 flowlines with 47,871 grid cells...

Building spatial index...
finished in 8.78s

Intersecting 99 features...
99
finished in 5.16s

Setting up reach data... (may take a few minutes for large grids)
finished in 16.73s

Dropping 216 reaches with length < 8.00 feet...

Repairing routing connections...
enforcing best segment numbering...

Setting up segment data...
Model grid information
unstructured grid
nnodes: 47,871
model length units: feet
crs: epsg:2243
bounds: 2709517.94, 1150502.67, 2772838.51, 1217702.18
active area defined by: all cells

Then I get an error. I followed the the execution of code into sfrdata.py to line 501. Code is
elif model is None or model.version == 'mf6':
model_ws = '.' if model is None else model.model_ws
m = fm.Modflow(model_ws=model_ws, structured=self.structured)

In my case self.structured is False and the code crashes. I don't see any other ways to specify an unstructured grid to pass into sfr = stream_lines.to_sfr(). Am I missing something?

[aleaf]

@porterma , no, doesn't look like you're missing anything. I think it's an issue with how SFRmaker is interacting with Flopy in sfrdata. It mostly uses the structure of "old-style" (pre-MODFLOW 6) SFR packages internally (i.e. segments and reaches), and converts this to MODFLOW 6 input on write. Part of that is creating a flopy.modflow.ModflowSfr2 package instance with an attached flopy.modflow.Modflow model instance. The fix might be as simple as just adding version='mfusg' to m = fm.Modflow(model_ws=model_ws, structured=self.structured). I'll look into it. As I mentioned, the handling of unstructured grids is largely untested, so there might be some other issues.

[porterma]

@aleaf, that did indeed fix the problem. So it appears that in my script I need to instantiate a flopy modflow model with version='mfusg', then create a dis package and possibly a bas6 file for the ibound, then call:
` stream_lines = lines.from_dataframe(streams_df, attr_length_units='feet', attr_height_units='feet',
epsg='2243')

grd = UnstructuredGrid.from_dataframe(grd_df, node_col='nodenumber', geometry_column='geometry', model_units='feet', proj_str='epsg:2243', )

sfr = stream_lines.to_sfr(grd, model=mdl)

sfr.set_streambed_top_elevations_from_dem(dem_file, dem_z_units='feet', method='cell polygons', smooth=True)

sfr.create_ModflowSfr2(model=mfusg)
sfr.ModflowSfr2.check()

sfr.create_mf6sfr(model=mf6)
sfr.write_package(filename='testsfrmaker.sfr', version='mf6')

`

[nikobenho]

My case is very similar to @porterma
A 3-layer disv quadtree grid refined around a single stream. Refinement is exactly the same in all layers and no risk of the stream cutting through the bottom of the first layer.

I tried the workaround suggested in this issue, but it may not work for the current version of SFRmaker?

I instantiate a Lines object using sfrmaker.Lines.from_shapefile and upon running sfr = custom_lines.to_sfr(grid, model=gwf) i get AttributeError: 'VertexGrid' object has no attribute 'nrow'

If i instead instantiate a regular flopy mf model i get AttributeError: 'NoneType' object has no attribute 'nrow'

Could you give a more detailed description on how you managed the workaround, or perhaps I'm just a lost cause? :)

I think the possibility to use this package with MF6 DISV grids would be a great feature. Unfortunately I have no idea how to contribute towards this goal.

[aleaf]

Hi @nikobenho, it sounds like SFRmaker might be trying to work with the flopy modelgrid attached to your flopy model, and also assuming that the flopy grid is structured. Can you post a more complete version of your script? Did you do the first two lines that @porterma has above:

grd = UnstructuredGrid.from_dataframe(grd_df, node_col='nodenumber', geometry_column='geometry', model_units='feet', proj_str='epsg:2243', )

sfr = stream_lines.to_sfr(grd, model=mdl)

This creates an SFRmaker unstructured grid instance from a DataFrame of shapely polygons representing the cells, and then passes that to the to_sfr() method.

[nikobenho]

Hi @aleaf, thanks for replying. I really appreciate it!

So my model consist of a 3-layer quadtree disv grid, which is refined along the creek (that I want to represent using SFR) as well as around a few wells.

I eventually figured out that I had worked with the Flopy modelgrid rather than sfrmaker.grid.UnstructuredGrid.

It got me a little bit further, but then i eventually got stuck again:

1. I read a polyline shapefile, this shapefile consist of a single small creek that i manually subdivided into 13 sections:

stream_lines = sfrmaker.Lines.from_shapefile(shapefile=os.path.join(shp_pth, 'orbacken_sfr.shp'),
                                             id_column='id',
                                             routing_column='toid',
                                             width1_column='width1',
                                             width2_column='width2',
                                             up_elevation_column='elev_up',
                                             dn_elevation_column='elev_dn',
                                             attr_length_units='meters',
                                             attr_height_units='meters',
                                             )

The DF of stream_lines looks as follows:

  | id | toid | asum1 | asum2 | width1 | width2 | elevup | elevdn | name | geometry
-- | -- | -- | -- | -- | -- | -- | -- | -- | -- | --
1 | 2 | 0 | 0 | 2.0 | 2.0 | 143.62 | 143.43 | 0 | LINESTRING (427626.783017383 6655133.216913054...
2 | 3 | 0 | 0 | 2.0 | 2.0 | 143.43 | 141.07 | 0 | LINESTRING (427533.6268522348 6655090.92536188...
3 | 4 | 0 | 0 | 2.0 | 2.0 | 141.07 | 138.36 | 0 | LINESTRING (427501.7888053609 6655037.43406598...
4 | 5 | 0 | 0 | 2.0 | 2.0 | 138.36 | 138.30 | 0 | LINESTRING (427486.689886272 6654953.711706178...
5 | 6 | 0 | 0 | 2.0 | 2.0 | 138.30 | 137.10 | 0 | LINESTRING (427460.9502861826 6654902.48086205...
6 | 7 | 0 | 0 | 2.0 | 2.0 | 137.10 | 136.86 | 0 | LINESTRING (427418.219992389 6654902.216278188...
7 | 8 | 0 | 0 | 2.0 | 2.0 | 136.86 | 135.87 | 0 | LINESTRING (427393.6011787876 6654918.36930097...
8 | 9 | 0 | 0 | 2.0 | 2.0 | 135.87 | 135.27 | 0 | LINESTRING (427343.5021803295 6654961.33002525...
9 | 10 | 0 | 0 | 2.0 | 2.0 | 135.27 | 134.85 | 0 | LINESTRING (427307.1679147501 6654967.88989706...
10 | 11 | 0 | 0 | 2.0 | 2.0 | 134.85 | 134.43 | 0 | LINESTRING (427220.5262247967 6654909.09812807...
11 | 12 | 0 | 0 | 2.0 | 2.0 | 134.43 | 133.93 | 0 | LINESTRING (427201.1738741575 6654935.41509163...
12 | 13 | 0 | 0 | 2.0 | 2.0 | 133.93 | 133.28 | 0 | LINESTRING (427087.0727703692 6654957.07941936...
13 | 14 | 0 | 0 | 2.0 | 2.0 | 133.28 | 133.10 | 0 | LINESTRING (427021.0590966716 6654931.99157751...

2. I import the quadtree grid shapefile using Geopandas:

import geopandas
grid_shp = os.path.join(temp_gridgen_ws, 'qtgrid.shp')
grd_df = geopandas.read_file(grid_shp)
grd_df.head()

The DF of the grid looks like this:

nodenumber | layer | row | col | child_loca | top | bottom | delr | delc | geometry
-- | -- | -- | -- | -- | -- | -- | -- | -- | --
1 | 0.0 | 40.0 | 19.0 | None | 135.437851 | 127.354347 | 4.0 | 4.0 | POLYGON ((426921.902 6654935.578, 426924.731 6...
2 | 0.0 | 40.0 | 20.0 | None | 135.425964 | 127.389160 | 4.0 | 4.0 | POLYGON ((426924.731 6654938.406, 426927.559 6...
3 | 0.0 | 40.0 | 21.0 | None | 135.513184 | 127.372437 | 4.0 | 4.0 | POLYGON ((426927.559 6654941.235, 426930.388 6...
4 | 0.0 | 40.0 | 22.0 | None | 135.412445 | 127.410713 | 4.0 | 4.0 | POLYGON ((426930.388 6654944.063, 426933.216 6...
5 | 0.0 | 40.0 | 23.0 | None | 135.426621 | 126.775185 | 4.0 | 4.0 | POLYGON ((426933.216 6654946.891, 426936.045 6...

3. I then instantiate a SFR unstructured grid:

grd = sfrmaker.grid.UnstructuredGrid.from_dataframe(
    grd_df,
    node_col='nodenumber',
    geometry_column='geometry',
    model_units='meters',
    crs='epsg:3006'
)

Typing grd yields:

Model grid information
unstructured grid
nnodes: 70,053
model length units: meters
crs: EPSG:3006
bounds: 426921.90, 6654607.48, 427665.78, 6655354.18
active area defined by: all cells

So far so good (I think). This is where my confusion starts. When i run sfr = stream_lines.to_sfr(grid=grd, model_length_units = 'meters') I get the following AssertionError:

SFRmaker version 0.8.0

Creating sfr dataset...
Model grid information
unstructured grid
nnodes: 70,053
model length units: meters
crs: EPSG:3006
bounds: 426921.90, 6654607.48, 427665.78, 6655354.18
active area defined by: all cells

None

Culling hydrography to active area...
starting lines: 13
remaining lines: 13
finished in 0.02s


Intersecting 13 flowlines with 70,053 grid cells...

Building spatial index...
finished in 7.48s

Intersecting 13 features...
13
finished in 1.15s

Setting up reach data... (may take a few minutes for large grids)
finished in 4.99s


Dropping 120 reaches with length < 0.05 meters...

Repairing routing connections...
enforcing best segment numbering...

Setting up segment data...
Model grid information
unstructured grid
nnodes: 70,053
model length units: meters
crs: EPSG:3006
bounds: 426921.90, 6654607.48, 427665.78, 6655354.18
active area defined by: all cells

----------------------------------------------------------------
AssertionError                 Traceback (most recent call last)
<ipython-input-50-6d72d5f5967e> in <module>
----> 1 sfr = stream_lines.to_sfr(grid=grd, model_length_units = 'meters')

c:\users\nat12nho\appdata\local\programs\python\python38\lib\site-packages\sfrmaker\lines.py in to_sfr(self, grid, active_area, isfr, model, model_length_units, model_time_units, minimum_reach_length, width_from_asum_a_param, width_from_asum_b_param, minimum_reach_width, consolidate_conductance, one_reach_per_cell, add_outlets, package_name, **kwargs)
    966         # and other output
    967         rd = rd[[c for c in SFRData.rdcols if c in rd.columns]].copy()
--> 968         sfrd = SFRData(reach_data=rd, segment_data=sd, grid=grid,
    969                        model=model, model_length_units=model_length_units,
    970                        model_time_units=model_time_units,

c:\users\nat12nho\appdata\local\programs\python\python38\lib\site-packages\sfrmaker\sfrdata.py in __init__(self, reach_data, segment_data, grid, model, isfr, model_length_units, model_time_units, enforce_increasing_nsegs, package_name, **kwargs)
    176 
    177         # have to set the model last, because it also sets up a flopy sfr package instance
--> 178         self.model = model  # attached flopy model instance
    179         # self._ModflowSfr2 = None # attached instance of flopy modflow_sfr2 package object
    180 

c:\users\nat12nho\appdata\local\programs\python\python38\lib\site-packages\sfrmaker\sfrdata.py in model(self, model)
    215 
    216         # update the sfr package object as well
--> 217         self.create_modflow_sfr2(model)
    218 
    219     @property

c:\users\nat12nho\appdata\local\programs\python\python38\lib\site-packages\sfrmaker\sfrdata.py in create_modflow_sfr2(self, model, const, isfropt, unit_number, ipakcb, istcb2, **kwargs)
    618         elif model is None or model.version == 'mf6':
    619             model_ws = '.' if model is None else model.model_ws
--> 620             m = fm.Modflow(modelname=self.package_name, model_ws=model_ws,
    621                            structured=self.structured)
    622             if model is not None:

c:\users\nat12nho\appdata\local\programs\python\python38\lib\site-packages\flopy\modflow\mf.py in __init__(self, modelname, namefile_ext, version, exe_name, structured, listunit, model_ws, external_path, verbose, **kwargs)
    142         # other than mfusg is specified
    143         if not self.structured:
--> 144             assert (
    145                 "mfusg" in self.version
    146             ), "structured=False can only be specified for mfusg models"

AssertionError: structured=False can only be specified for mfusg models

I then try to pass my MF6 model instance like this: sfr = stream_lines.to_sfr(grid=grd, model=gwf, model_length_units = 'meters')
This yields AttributeError: 'VertexGrid' object has no attribute 'nrow'

I then try to do a hacky gwf._version = 'mfusg' which yields AttributeError: bas6. A bas6 package can't afaik be added to a mf6 model, so that's a dead end.

I then try to pass the base-modelgrid (an instance of flopy.modflow.Modflow() with the flopy.modflow.ModflowDis package), used by gridgen to create the disv-grid. Although this model is similar in terms of coordinates, rotation and elevations, it does not contain the same amount of cells and refinement as does the disv-grid created through gridgen. Before passing the base model i run ml.version = 'mfusg' and add a bas package using bas = flopy.modflow.ModflowBas(ml, ibound=1). This solution allow SFRmaker to run with the following output:

SFRmaker version 0.8.0

Creating sfr dataset...

c:\users\nat12nho\appdata\local\programs\python\python38\lib\site-packages\sfrmaker\gis.py:25: UserWarning: The epsg argument is deprecated, use crs instead.
  warnings.warn('The epsg argument is deprecated, use crs instead.')

Model grid information
structured grid
nnodes: 35,000
nlay: 1
nrow: 175
ncol: 200
model length units: undefined
crs: EPSG:3006
bounds: 426755.03, 6654500.00, 427815.69, 6655560.66
active area defined by: all cells

MODFLOW 3 layer(s) 175 row(s) 200 column(s) 1 stress period(s)

Culling hydrography to active area...
starting lines: 13
remaining lines: 13
finished in 0.00s


Intersecting 13 flowlines with 35,000 grid cells...

Building spatial index...
finished in 4.23s

Intersecting 13 features...
13
finished in 0.07s

Setting up reach data... (may take a few minutes for large grids)
finished in 0.08s


Dropping 10 reaches with length < 0.20 meters...

Repairing routing connections...
enforcing best segment numbering...

Setting up segment data...
Model grid information
structured grid
nnodes: 35,000
nlay: 1
nrow: 175
ncol: 200
model length units: undefined
crs: EPSG:3006
bounds: 426755.03, 6654500.00, 427815.69, 6655560.66
active area defined by: all cells


Time to create sfr dataset: 5.10s

The rows and columns stated above are from what I can tell based on the structure of the basemodel passed as the model argument in to_sfr and not based on the grd. Is this a problem?

By running sfr.reach_data.head() i can see the following output:

  | rno | node | k | i | j | iseg | ireach | rchlen | width | slope | ... | thts | thti | eps | uhc | outreach | outseg | asum | line_id | name | geometry
-- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | --
1 | 25978 | 0 | 129 | 178 | 1 | 1 | 2.475307 | 2.0 | 0.001983 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 2 | 2 | 0 | 1 | 0 | LINESTRING (427626.783017383 6655133.216913054...
2 | 25977 | 0 | 129 | 177 | 1 | 2 | 3.150006 | 2.0 | 0.001259 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 3 | 2 | 0 | 1 | 0 | LINESTRING (427624.4111155691 6655132.50894084...
3 | 25777 | 0 | 128 | 177 | 1 | 3 | 1.396430 | 2.0 | 0.003654 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 4 | 2 | 0 | 1 | 0 | LINESTRING (427621.3926997556 6655131.60799523...
4 | 25776 | 0 | 128 | 176 | 1 | 4 | 4.451344 | 2.0 | 0.001502 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 5 | 2 | 0 | 1 | 0 | LINESTRING (427620.0546052901 6655131.20859687...
5 | 25775 | 0 | 128 | 175 | 1 | 5 | 3.210688 | 2.0 | 0.001205 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 6 | 2 | 0 | 1 | 0 | LINESTRING (427615.8463303364 6655129.76001757...

sfr.segment_data.head() looks like this:

  | per | nseg | icalc | outseg | iupseg | iprior | nstrpts | flow | runoff | etsw | ... | uhc1 | hcond2 | thickm2 | elevdn | width2 | depth2 | thts2 | thti2 | eps2 | uhc2
-- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | --
0 | 1 | 1 | 2 | 0 | 0 | 0 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 143.429993 | 2.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
0 | 2 | 1 | 3 | 0 | 0 | 0 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 141.070007 | 2.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
0 | 3 | 1 | 4 | 0 | 0 | 0 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 138.360001 | 2.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
0 | 4 | 1 | 5 | 0 | 0 | 0 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 138.300003 | 2.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
0 | 5 | 1 | 6 | 0 | 0 | 0 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 137.100006 | 2.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0

len(sfr.reach_data) yields: 346

Running sfr.run_diagnostics() tells me that all checks are passed except No DIS package or SpatialReference object; cannot check reach proximities., 12 model cells with multiple non-zero SFR conductances found. and 159 reaches encountered with streambed above model top.. My guess is that this last error is the most problematic? Basically it's half of the reaches in the DF, and the differences between streambed and model top differ between 0.0015 to 1.02 meters.

I then apply sfr.set_streambed_top_elevations_from_dem(os.path.join(raster_pth, 'groundsurface.tif'), dem_z_units='meters', method='cell polygons', smooth=True) which reduces the number of reaches that contain a streambed above model top down to 2 (much easier to handle manually). However, it also introduces the following warnings which were not present before sampling the top elevations from the DEM:

Checking for streambed tops of less than -10...
isfropt setting of 1,2 or 3 requires strtop information!

Checking for streambed tops of greater than 15000...
isfropt setting of 1,2 or 3 requires strtop information!

In the docs i find the following:

isfropt0_to_1()[source]
    transfer isfropt=0 segment properties to reaches, using linear interpolation.

But I'm not really sure what to make of that. At any rate, I go ahead and run sfr.assign_layers() and sfr.reach_data.head(), which yields:

rno | node | k | i | j | iseg | ireach | rchlen | width | slope | ... | thts | thti | eps | uhc | outreach | outseg | asum | line_id | name | geometry
-- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | --
1 | 25978 | 0 | 129 | 178 | 1 | 1 | 2.475307 | 2.0 | 0.001983 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 2 | 2 | 0 | 1 | 0 | LINESTRING (427626.783017383 6655133.216913054...
2 | 25977 | 0 | 129 | 177 | 1 | 2 | 3.150006 | 2.0 | 0.001259 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 3 | 2 | 0 | 1 | 0 | LINESTRING (427624.4111155691 6655132.50894084...
3 | 25777 | 0 | 128 | 177 | 1 | 3 | 1.396430 | 2.0 | 0.003654 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 4 | 2 | 0 | 1 | 0 | LINESTRING (427621.3926997556 6655131.60799523...
4 | 25776 | 0 | 128 | 176 | 1 | 4 | 4.451344 | 2.0 | 0.001502 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 5 | 2 | 0 | 1 | 0 | LINESTRING (427620.0546052901 6655131.20859687...
5 | 25775 | 0 | 128 | 175 | 1 | 5 | 3.210688 | 2.0 | 0.001205 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 6 | 2 | 0 | 1 | 0 | LINESTRING (427615.8463303364 6655129.76001757...

The layer information (k) is set to 0 in 320 cells, and 1 in 26 cells. Shouldn't this be the same number for all cells (i.e. SFR is supposed to be defined for the top layer only)?

Given the assumption that I'm doing things right by passing the basemodel to to_sfr(), isn't it a problem that values in SFR are stored in i j columns rather than layer number + cellid (as in disv?)

Running sfr.write_package(version='mf6') yields the sfr input file modflowtest.sfr, which, including the five first reach data rows looks like this:

# MODFLOW-6 SFR input; created by SFRmaker v. 0.8.0

BEGIN Options
  save_flows
  BUDGET FILEOUT modflowtest.sfr.cbc
  STAGE FILEOUT modflowtest.sfr.stage.bin
  unit_conversion  86400.0
  auxiliary line_id
END Options

BEGIN Dimensions
  NREACHES 346
END Dimensions

BEGIN Packagedata
#rno k i j rlen rwid rgrd rtp rbth rhk man ncon ustrf ndv line_id
1 1 130 179 2.4753065 2.0 0.0019830140390995343 143.62547302246094 1 1.0 0.03700000047683716 1 1.0 0 1
2 1 130 178 3.1500058 2.0 0.0012594120204762016 143.49830627441406 1 1.0 0.03700000047683716 2 1.0 0 1
3 1 129 178 1.3964297 2.0 0.0036540934709742828 143.49830627441406 1 1.0 0.03700000047683716 2 1.0 0 1
4 1 129 177 4.4513445 2.0 0.001501969434360793 143.4317626953125 1 1.0 0.03700000047683716 2 1.0 0 1
5 1 129 176 3.210688 2.0 0.001204947673459267 143.4317626953125 1 1.0 0.03700000047683716 2 1.0 0 1

For a disv grid, the cellid is identified by layer and cell2d number, is it possible to translate from k i j to cell2d?

Sorry for the long reply, but perhaps it makes it easier to identify where I fail.

[BrockOliver]

Hi all,

Sorry to open this issue back up after several years, but I have encountered the same problem as @nikobenho: I generated an unstructured quadtree grid in gridgen as a .gsf and tried to use this grid with NHD Plus data to generate the SFR input data for my model with SFRmaker. I can provide more specifics, but the hang up that I cannot figure out is the same assertion error:

structured=False can only be specified for mfusg models

My first question is this: is it even possible to use SFRmaker with unstructured grid data? If so, Is there any progress on this front, and if so can somebody please direct me to a solution? Thank you!

-Brett

[kzeiler]

I have a suggestion for a "fix" (maybe more of a workaround). In the .to_sfr() function in lines.py between the printing of grid and model information to the screen and the check if flowlines need to be re-projected. But basically do a check if the grid is unstructured (or not structured) and there's either no model object or an mf6 model object that was supplied, and if so make a temporary/dummy mfusg model with structured=False.

        if not grid._structured and (model is None or model.version == 'mf6'):
            print('Notice - no model object or an mf6 model object was supplied for an unstructured grid, so "dummy" MFUSG model will be created for SFRData object returned by to_sfr().')
            model = flopy.mfusg.MfUsg(structured=False)

I haven't done any extensive testing, but it "just works" (so far) and would prevent users from needing to know workaround of needing an mfusg model object. The notice message could maybe use some streamlining and/or better explanation.