Assistance with exporting 3D layer data to shapefiles in MODFLOW 6

[YozeGholami]

Hi
I hope this message finds you well. I am currently working on exporting the bottom elevations of each layer to a shapefile in MODFLOW 6 using the following FloPy documentation code: (ml.dis.botm.export(pth, "botm.shp"). However, I have encountered an issue where the bottom elevation values exported to the shapefile are identical for every layer.

Could you assist me in determining the correct approach to export 3D results to a shapefile accurately?

Best regards,

Amin

[sebmeng]

Hello Amin,
Could you maybe share your model? Maybe there is a problem, and it indeed has all of the layers with the same elevation? You could also export it as a VTK (3D mesh) that you can load and see after in ParaView or using the PyVista library.
You can also plot 2D sections of your model to check it, I reccomend this example notebook

Cheers!
Sebastian

[YozeGholami]

Hello Sebastian,

Thank you for your suggestions. I have re-checked the exported data using the MODFLOW 6 examples available at modflow6-examples.readthedocs.io, and it appears this issue occurs consistently. If you check any model with more than one layer and import the resulting shapefile into ArcGIS Pro, you will notice that the bottom elevation of all layers is the same.

To clarify, my primary objective is not merely visualization but to extract shapefiles for use in ArcGIS Pro. Therefore, while exporting the model as a VTK for viewing in ParaView or PyVista might be helpful for other purposes, it does not meet my current needs. I require the shapefiles for conducting further spatial analysis within GIS software.

Could we explore any alternative methods or adjustments to the FloPy export functions that might ensure the accurate extraction of layer-specific data into shapefiles?

Thank you for your assistance!

Best regards,
Amin

[sebmeng]

if you export your whole simulation object with sim.export("sim.shp"), load it to GIS, and then do you still see that all of the bottoms of the layers are the same?

[YozeGholami]

Thank you for your comment. I attempted to use the following code:

"sim_sub.export(os.path.join(ws_sub, "model.shp")) ".

However, I encountered an AttributeError: 'export'.
Best regards,
Amin

[martclanor]

Hi Amin,

I can confirm this happening also for vertex grids. Iiuc, this has something to do with how the creation of flopy.mf6.data.mfdatastorage.LayerStorage object is triggered within MultiList.build_list method wherein layer 0 is constantly being passed as an argument (instead of layer 0, layer 1, etc.). I'll try to submit a PR on this.

https://github.com/martclanor/flopy/blob/a0e9219407b7be208d38f4f902cd2b5b96da1351/flopy/utils/datautil.py#L587-L610

Regards,
Mart

[martclanor]

It seems that this issue can be avoided if you would export instead at the model level or package level (similar to what @sebmeng suggested above), for example,

ml.dis.export(path / "dis.shp")
gwf.export(path / "model.shp")

So maybe, for your use case that would be the faster solution.

There is a difference in behavior because exporting the 3d array, ml.dis.botm.export makes use of flopy/export/utils.array3d_export while exporting at model or package level makes use of flopy/export/shapefile_utils/model_attributes_to_shapefile. The two handles the indexing of the MFArray differently.

[YozeGholami]

Thank you for your recommendation. Initially, I tried using "gwf.export(path / 'model.shp')" but it only gave me some polygons that I couldn't use. Then I switched to "ml.dis.export(path / 'dis.shp')" which was very helpful, allowing me to extract data from the bottom of different layers. I appreciate your helpful comment.

Best regards,
Amin

[wpbonelli]

A fix was merged with #2310, if you'd like to confirm the problem is solved on your end @YozeGholami.

[YozeGholami]

Thank you for the update! I'll pull the latest changes and test it on my end to confirm that the issue is resolved. I'll let you know if everything works as expected.

Best regards,
Amin