flopy icon indicating copy to clipboard operation
flopy copied to clipboard

bug: cannot export package stress period data to shapefile

Open aestrad-intera opened this issue 6 months ago • 10 comments

The following lines of code:

sim = flopy.mf6.MFSimulation.load(sim_ws='run_temp_ws') gwf = sim.get_model() mg = gwf.modelgrid mg.set_coord_info(xoff=5677631.08875377, yoff=18360574.17093341, angrot=50) chd = gwf.chd chd.stress_period_data.export(os.path.join('..', '..', 'arcpro', 'shp', 'pumptest_chd.shp'), sparse=True)

Produce this error:

896 from ..export.shapefile_utils import recarray2shp
897 from ..utils.geometry import Polygon

--> 899 df = mfl.get_dataframe(squeeze=squeeze) 900 if "kper" in kwargs or df is None: 901 ra = mfl[kper]

TypeError: MFPandasTransientList.get_dataframe() got an unexpected keyword argument 'squeeze'

It appears that the stress_period_data.export utility is using a deprecated squeeze parameter. What other method is available for exporting list package data to shapefiles?

aestrad-intera avatar Jun 13 '25 19:06 aestrad-intera

Thanks for reporting @aestrad-intera there is a patch in the develop branch if you'd like to check it on your end

wpbonelli avatar Jun 13 '25 22:06 wpbonelli

@wpbonelli I included you modifications in mflist_export() utils.py but now I am getting a pandas error upon running the same code:

Cell In[2], line 2 1 chd = gwf.chd ----> 2 chd.stress_period_data.export(os.path.join('..', '..', 'arcpro', 'shp', 'pumptest_chd.shp'), sparse=True)

File dependencies\flopy\mf6\data\mfdata.py:375, in MFData.export(self, f, **kwargs) 373 return utils.transient2d_export(f, self, **kwargs) 374 elif self.data_type == DataType.transientlist: --> 375 return utils.mflist_export(f, self, **kwargs) 376 return utils.transient2d_export(f, self, **kwargs)

File dependencies\flopy\export\utils.py:922, in mflist_export(f, mfl, **kwargs) 915 row_index_name = "i" if "i" in df.columns else "cellid_row" 916 col_index_name = "j" if "j" in df.columns else "cellid_column" 917 verts = np.array( 918 [ 919 # modelgrid.get_cell_vertices(i, df.j.values[ix]) 920 # for ix, i in enumerate(df.i.values) 921 modelgrid.get_cell_vertices(i, df[col_index_name].values[ix]) --> 922 for ix, i in enumerate(df[row_index_name].values) 923 ] 924 ) 925 ra = df.to_records(index=False) 926 crs = kwargs.get("crs", None)

File ~\AppData\Local\anaconda3\envs\mesa\lib\site-packages\pandas\core\frame.py:4102, in DataFrame.getitem(self, key) 4100 if self.columns.nlevels > 1: 4101 return self._getitem_multilevel(key) -> 4102 indexer = self.columns.get_loc(key) 4103 if is_integer(indexer): 4104 indexer = [indexer]

File ~\AppData\Local\anaconda3\envs\mesa\lib\site-packages\pandas\core\indexes\base.py:3812, in Index.get_loc(self, key) 3807 if isinstance(casted_key, slice) or ( 3808 isinstance(casted_key, abc.Iterable) 3809 and any(isinstance(x, slice) for x in casted_key) 3810 ): 3811 raise InvalidIndexError(key) -> 3812 raise KeyError(key) from err 3813 except TypeError: 3814 # If we have a listlike key, _check_indexing_error will raise 3815 # InvalidIndexError. Otherwise we fall through and re-raise 3816 # the TypeError. 3817 self._check_indexing_error(key)

KeyError: 'cellid_row'

aestrad-intera avatar Jun 14 '25 14:06 aestrad-intera

@aestrad-intera to clarify how to reproduce this, could you pull the develop branch and have a look at the test I added in #2532? The only difference I see is setting coordinate info on the grid which is inconsequential here.

wpbonelli avatar Jun 18 '25 13:06 wpbonelli

@wpbonelli The model that I am using uses DISU package, could it be that this package is not covered by the patch?

aestrad-intera avatar Jun 18 '25 16:06 aestrad-intera

yeah, disv/u weren't covered. sorry @aestrad-intera we are discovering the state of all this as we go. does #2566 solve it for you?

wpbonelli avatar Aug 04 '25 13:08 wpbonelli

@aestrad-intera please reopen if this does not solve the issue for you

wpbonelli avatar Aug 07 '25 14:08 wpbonelli

It did not.

On Mon, Aug 4, 2025, 8:59 AM wpbonelli @.***> wrote:

wpbonelli left a comment (modflowpy/flopy#2531) https://github.com/modflowpy/flopy/issues/2531#issuecomment-3150837287

yeah, disv/u weren't covered. sorry @aestrad-intera https://github.com/aestrad-intera we are discovering the state of all this as we go. does #2566 https://github.com/modflowpy/flopy/pull/2566 solve it for you?

— Reply to this email directly, view it on GitHub https://github.com/modflowpy/flopy/issues/2531#issuecomment-3150837287, or unsubscribe https://github.com/notifications/unsubscribe-auth/A52MYOPVPC76MUV7OUSDDIT3L5RL5AVCNFSM6AAAAAB7I4J4V6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZTCNJQHAZTOMRYG4 . You are receiving this because you were mentioned.Message ID: @.***>

aestrad-intera avatar Aug 07 '25 14:08 aestrad-intera

Is it possible for you to share your model so we can reproduce the exact issue you're seeing?

wpbonelli avatar Aug 07 '25 15:08 wpbonelli

@aestrad-intera

I'd recommend trying to use the .geo_dataframe method on the modelgrid class and then adding the data you need to export for now. Here's an example.

sim = flopy.mf6.MFSimulation.load(sim_ws="my_disu_simulation_folder")
gwf = sim.get_model()
chd = gwf.chd

modelgrid = gwf.modelgird

# create a geopandas Geodataframe
gdf = modelgrid.geo_dataframe

data = chd.stress_period_data.data
kper = 0
ra = data[kper]

# note: if ra is a numpy recarray, convert to pandas dataframe. Otherwise skip this step
df = pd.DataFrame(ra)

# this is for DISU grids, DISV and DIS needs some additional logic to convert cellid to node
df["node"] = [i[0] + 1 for i in df.cellid.values]
gdf = pd.merge(gdf, df, how="inner", on="node")
gdf.to_file("my_shapefile.shp")

Note: this is a hack to try to help you move forward. I'm working on getting geopandas GeoDataFrame support added to data objects as a replacement for the existing shapefile export utilities, so in the near future you'll be able to do something like this gdf = chd.stress_period_data.to_geo_dataframe(kper=0, sparse=True) or gdf = chd.to_geo_dataframe(kper=0, sparse=True).

jlarsen-usgs avatar Aug 07 '25 16:08 jlarsen-usgs

@aestrad-intera any luck with @jlarsen-usgs's suggestion? and/or can you share a model we can reproduce this with?

wpbonelli avatar Dec 02 '25 18:12 wpbonelli