We did this tutorial since the documentation or examples on the topic were not available. Flopy can export the model grid and model attributes to shapefile with a coordinate system of reference for the three types of discretizations of MODFLOW 6. We have done an applied example to export the grid of a discretized by vertices model (DISV) to the ESRI Shapefile model. The tutorial also show options for the final grid representation in geopandas.

The reference model used for this tutorial can be found in this link:

MODFLOW 6 and MODPATH7 used to simulate regional groundwater flow and advective transport of per- and polyfluoroalkyl substances, Joint Base McGuire-Dix-Lakehurst and vicinity, New Jersey, 2018

Tutorial

Code

#!pip install ipympl
import os, re, time
import flopy
import numpy as np
import geopandas as gpd

Working with the system of reference and model extension

modelRef = gpd.read_file('../Shp/ofr2022_1112.shp')
modelRef.crs

<Projected CRS: PROJCS["NAD83 / BLM 18N (ftUS)",GEOGCS["NAD83",DAT ...>
Name: NAD83 / BLM 18N (ftUS)
Axis Info [cartesian]:
- [east]: Easting (US survey foot)
- [north]: Northing (US survey foot)
Area of Use:
- undefined
Coordinate Operation:
- name: unnamed
- method: Transverse Mercator
Datum: North American Datum 1983
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

#filter active zone
active = modelRef[modelRef.AREA == 'Active']
active.plot()

<Axes: >

Open MODFLOW 6 simulation and model

#define some model paths
simName = "mfsim.nam"
simWs = "../Model/JBMDLPFAS-MF/"
exeName = "../Bin/mf6.exe"

#open simulation
sim = flopy.mf6.MFSimulation.load(simName,
                                  exe_name=exeName, 
                                  sim_ws=simWs)

loading simulation...
  loading simulation name file...
  loading tdis package...
  loading model gwf6...
    loading package disv...
WARNING: Unable to resolve dimension of ('gwf6', 'disv', 'cell2d', 'cell2d', 'icvert') based on shape "ncvert".
    loading package ic...
WARNING: Block "options" is not a valid block name for file type ic.
    loading package npf...
    loading package obs...
    loading package gnc...
    loading package oc...
    loading package wel...
    loading package drn...
    loading package maw...
    loading package rch...
  loading solution package modflow...

sim.model_names

['modflow']

gwf = sim.get_model('modflow')

Apply / verify geospatial information to model grid

gwf.modelgrid.extent

(1678000.0, 1896000.0, 14470000.0, 14604000.0)

gwf.modelgrid.crs = 4438
gwf.modelgrid.crs

<Projected CRS: EPSG:4438>
Name: NAD83 / BLM 18N (ftUS)
Axis Info [cartesian]:
- X[east]: Easting (US survey foot)
- Y[north]: Northing (US survey foot)
Area of Use:
- name: United States (USA) - between 78°W and 72°W onshore and offshore - Connecticut; Delaware; Maryland; Massachusetts; New Hampshire; New Jersey; New York; North Carolina; Pennsylvania; Virginia; Vermont.
- bounds: (-78.0, 28.28, -72.0, 45.03)
Coordinate Operation:
- name: BLM zone 18N (US survey feet)
- method: Transverse Mercator
Datum: North American Datum 1983
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

Export model and model grid as shapefiles

gridLines = '../Shp/Model/gridLines.shp'
flopy.export.shapefile_utils.write_gridlines_shapefile(gridLines,
                                                  gwf.modelgrid)

gridLinesDf = gpd.read_file(gridLines)
gridLinesDf.head()

	number	geometry
0	0	LINESTRING (1678000.000 14603000.000, 1678000....
1	1	LINESTRING (1678000.000 14604000.000, 1679000....
2	2	LINESTRING (1679000.000 14604000.000, 1679000....
3	3	LINESTRING (1679000.000 14603000.000, 1678000....
4	4	LINESTRING (1679000.000 14603000.000, 1679000....

#get lines on the active part
activeGridLines = gridLinesDf.intersection(active.iloc[0].geometry)

#plot a part of the grid
activeGridLines.iloc[:100000].plot(figsize=(12,8))

<Axes: >

#plot the whole grid
activeGridLines.plot(figsize=(12,8),lw=.1,alpha=0.5)

<Axes: >

#another plot of the meshgrid
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(12,8))
activeGridLines.plot(ax=ax)
ax.set_ylim(14540000,14560000)
ax.set_xlim(1720000,1740000)