CtrlK
BlogDocsLog inGet started
Tessl Logo

tessl/pypi-flopy

FloPy is a Python package to create, run, and post-process MODFLOW-based models

66

1.13x
Overview
Eval results
Files

discretization.mddocs/

Grid and Discretization Classes

This module provides grid classes for structured, unstructured, and vertex-based discretizations with coordinate transformations, intersection utilities, and export capabilities. These classes form the spatial foundation for all FloPy models, handling geometric calculations, coordinate systems, and spatial data management.

Base Grid Class

Grid

Abstract base class that provides common functionality for all grid types.

class Grid:
    """Base class for structured or unstructured model grids"""
    def __init__(
        self,
        grid_type: str = None,
        top: ArrayData = None,
        botm: ArrayData = None,
        idomain: ArrayData = None,
        lenuni: int = 2,
        epsg: int = None,
        proj4: str = None,
        prj: str = None,
        xoff: float = 0.0,
        yoff: float = 0.0,
        angrot: float = 0.0,
        **kwargs
    ): ...
    
    def is_valid(self) -> bool:
        """Check if grid is valid and complete"""
        ...
    
    def is_complete(self) -> bool:
        """Check if all required grid data is present"""
        ...
    
    @property
    def grid_type(self) -> str:
        """Grid type identifier"""
        ...
    
    @property
    def extent(self) -> tuple[float, float, float, float]:
        """Grid spatial extent (xmin, xmax, ymin, ymax)"""
        ...
    
    @property
    def shape(self) -> tuple[int, ...]:
        """Grid dimensions"""
        ...
    
    @property
    def size(self) -> int:
        """Total number of cells"""
        ...
    
    @property
    def xcellcenters(self) -> np.ndarray:
        """X-coordinates of cell centers"""
        ...
    
    @property
    def ycellcenters(self) -> np.ndarray:
        """Y-coordinates of cell centers"""
        ...
    
    @property
    def zcellcenters(self) -> np.ndarray:
        """Z-coordinates of cell centers"""
        ...
    
    @property
    def xvertices(self) -> np.ndarray:
        """X-coordinates of cell vertices"""
        ...
    
    @property
    def yvertices(self) -> np.ndarray:
        """Y-coordinates of cell vertices"""
        ...
    
    @property
    def zvertices(self) -> np.ndarray:
        """Z-coordinates of cell vertices"""
        ...
    
    def get_coords(self, x: float, y: float) -> tuple[float, float]:
        """Transform coordinates to model coordinate system"""
        ...
    
    def get_local_coords(self, x: float, y: float) -> tuple[float, float]:
        """Transform model coordinates to local coordinate system"""
        ...
    
    def intersect(
        self,
        x: ArrayData,
        y: ArrayData,
        z: ArrayData = None,
        local: bool = False,
        forgive: bool = False
    ) -> np.ndarray:
        """Find grid cells intersected by points"""
        ...
    
    def cross_section_vertices(
        self,
        line: dict,
        local: bool = True
    ) -> np.ndarray:
        """Extract vertices for cross-section along line"""
        ...
    
    def write_shapefile(
        self,
        filename: str,
        array_dict: dict = None,
        **kwargs
    ) -> None:
        """Export grid geometry to shapefile"""
        ...
    
    @classmethod
    def from_binary_grid_file(
        cls,
        file_path: str,
        **kwargs
    ) -> 'Grid':
        """Load grid from binary grid file"""
        ...
    
    def plot(self, **kwargs) -> object:
        """Plot the grid"""
        ...

Structured Grid

StructuredGrid

Regular rectangular grid discretization for traditional MODFLOW models.

class StructuredGrid(Grid):
    """Regular rectangular grid discretization"""
    def __init__(
        self,
        delc: ArrayData = None,
        delr: ArrayData = None,
        top: ArrayData = None,
        botm: ArrayData = None,
        idomain: ArrayData = None,
        lenuni: int = 2,
        epsg: int = None,
        proj4: str = None,
        prj: str = None,
        xoff: float = 0.0,
        yoff: float = 0.0,
        angrot: float = 0.0,
        nlay: int = None,
        nrow: int = None,
        ncol: int = None,
        **kwargs
    ): ...
    
    @property
    def nlay(self) -> int:
        """Number of layers"""
        ...
    
    @property
    def nrow(self) -> int:
        """Number of rows"""
        ...
    
    @property
    def ncol(self) -> int:
        """Number of columns"""
        ...
    
    @property
    def delr(self) -> np.ndarray:
        """Column spacing (along rows)"""
        ...
    
    @property
    def delc(self) -> np.ndarray:
        """Row spacing (along columns)"""
        ...
    
    @property
    def delz(self) -> np.ndarray:
        """Layer thickness"""
        ...
    
    @property
    def cell_thickness(self) -> np.ndarray:
        """Cell thickness array"""
        ...
    
    @property
    def saturated_thickness(self, array: ArrayData) -> np.ndarray:
        """Calculate saturated thickness given head array"""
        ...
    
    def is_regular_x(self, rtol: float = 1e-5) -> bool:
        """Check if column spacing is regular"""
        ...
    
    def is_regular_y(self, rtol: float = 1e-5) -> bool:
        """Check if row spacing is regular"""
        ...
    
    def is_regular_z(self, rtol: float = 1e-5) -> bool:
        """Check if layer spacing is regular"""
        ...
    
    def get_lrc(self, nodes: list[int]) -> list[tuple[int, int, int]]:
        """Convert node numbers to layer-row-column indices"""
        ...
    
    def get_node(self, lrc: list[tuple[int, int, int]]) -> list[int]:
        """Convert layer-row-column indices to node numbers"""
        ...
    
    def get_vertices(self, i: int, j: int) -> list[tuple[float, float]]:
        """Get vertices for cell at row i, column j"""
        ...
    
    def get_cellcenters(self) -> tuple[np.ndarray, np.ndarray]:
        """Get cell center coordinates"""
        ...

Parameters:

  • delr (ArrayData): Column spacing array (length ncol)
  • delc (ArrayData): Row spacing array (length nrow)
  • top (ArrayData): Top elevation of layer 1
  • botm (ArrayData): Bottom elevations of each layer
  • nlay (int): Number of layers
  • nrow (int): Number of rows
  • ncol (int): Number of columns

Unstructured Grid

UnstructuredGrid

Flexible unstructured mesh discretization for complex geometries.

class UnstructuredGrid(Grid):
    """Flexible unstructured mesh discretization"""
    def __init__(
        self,
        vertices: ArrayData = None,
        iverts: ArrayData = None,
        xcenters: ArrayData = None,
        ycenters: ArrayData = None,
        top: ArrayData = None,
        botm: ArrayData = None,
        idomain: ArrayData = None,
        lenuni: int = 2,
        epsg: int = None,
        proj4: str = None,
        prj: str = None,
        xoff: float = 0.0,
        yoff: float = 0.0,
        angrot: float = 0.0,
        ncpl: ArrayData = None,
        **kwargs
    ): ...
    
    @property
    def nnodes(self) -> int:
        """Number of nodes"""
        ...
    
    @property
    def nvertices(self) -> int:
        """Number of vertices"""
        ...
    
    @property
    def ncpl(self) -> np.ndarray:
        """Number of cells per layer"""
        ...
    
    @property
    def vertices(self) -> np.ndarray:
        """Vertex coordinates array"""
        ...
    
    @property
    def iverts(self) -> list[list[int]]:
        """Vertex indices for each cell"""
        ...
    
    def get_vertices(self, cellid: int) -> np.ndarray:
        """Get vertices for specified cell"""
        ...
    
    def get_cell_vertices(self, cellid: int) -> list[tuple[float, float]]:
        """Get vertex coordinates for cell"""
        ...
    
    def get_cellcenters(self) -> tuple[np.ndarray, np.ndarray]:
        """Get cell center coordinates"""
        ...

Parameters:

  • vertices (ArrayData): Vertex coordinate array
  • iverts (ArrayData): Cell-to-vertex connectivity
  • xcenters (ArrayData): Cell center x-coordinates
  • ycenters (ArrayData): Cell center y-coordinates
  • ncpl (ArrayData): Number of cells per layer

Vertex Grid

VertexGrid

DISV-type vertex-based grid discretization combining structured layers with unstructured horizontal discretization.

class VertexGrid(Grid):
    """DISV-type vertex-based grid discretization"""
    def __init__(
        self,
        vertices: ArrayData = None,
        cell2d: ArrayData = None,
        top: ArrayData = None,
        botm: ArrayData = None,
        idomain: ArrayData = None,
        lenuni: int = 2,
        epsg: int = None,
        proj4: str = None,
        prj: str = None,
        xoff: float = 0.0,
        yoff: float = 0.0,
        angrot: float = 0.0,
        nlay: int = None,
        ncpl: int = None,
        nvert: int = None,
        **kwargs
    ): ...
    
    @property
    def nlay(self) -> int:
        """Number of layers"""
        ...
    
    @property
    def ncpl(self) -> int:
        """Number of cells per layer"""
        ...
    
    @property
    def nvert(self) -> int:
        """Number of vertices"""
        ...
    
    @property
    def vertices(self) -> np.ndarray:
        """Vertex coordinates"""
        ...
    
    @property
    def cell2d(self) -> np.ndarray:
        """2D cell definition array"""
        ...
    
    def get_vertices(self, cellid: int) -> np.ndarray:
        """Get vertices for cell"""
        ...
    
    def get_cell_vertices(self, cellid: int) -> list[tuple[float, float]]:
        """Get vertex coordinates for cell"""
        ...

Parameters:

  • vertices (ArrayData): Vertex coordinate pairs
  • cell2d (ArrayData): Cell definition with vertex indices
  • nlay (int): Number of layers
  • ncpl (int): Number of cells per layer
  • nvert (int): Number of vertices

Temporal Discretization

ModelTime

Temporal discretization management for transient simulations.

class ModelTime:
    """Temporal discretization management"""
    def __init__(
        self,
        perioddata: list[tuple] = None,
        time_units: str = 'days',
        start_datetime: str = None,
        **kwargs
    ): ...
    
    @property
    def nper(self) -> int:
        """Number of stress periods"""
        ...
    
    @property
    def perlen(self) -> list[float]:
        """Length of each stress period"""
        ...
    
    @property
    def nstp(self) -> list[int]:
        """Number of time steps in each stress period"""
        ...
    
    @property
    def tsmult(self) -> list[float]:
        """Time step multipliers"""
        ...
    
    @property
    def steady(self) -> list[bool]:
        """Steady state flags for each stress period"""
        ...
    
    @property
    def totim(self) -> list[float]:
        """Total elapsed times"""
        ...
    
    def get_totim(self, kstp: int, kper: int) -> float:
        """Get total time for time step and stress period"""
        ...
    
    def get_kstpkper(self, totim: float) -> tuple[int, int]:
        """Get time step and stress period for total time"""
        ...

Parameters:

  • perioddata (list[tuple]): List of (perlen, nstp, tsmult) tuples
  • time_units (str): Time units ('days', 'hours', 'minutes', 'seconds', 'years')
  • start_datetime (str): Starting date/time in ISO format

Caching and Performance

CachedData

Caching mechanism for expensive grid calculations.

class CachedData:
    """Caching mechanism for grid data"""
    def __init__(
        self,
        data_func: callable,
        cache_size: int = 128,
        **kwargs
    ): ...
    
    def get_data(self, *args, **kwargs) -> object:
        """Get cached or computed data"""
        ...
    
    def clear_cache(self) -> None:
        """Clear the cache"""
        ...
    
    @property
    def cache_info(self) -> dict:
        """Cache statistics"""
        ...

Grid Utility Functions

def array_at_verts_basic2d(
    array: ArrayData,
    grid: StructuredGrid
) -> np.ndarray:
    """Convert cell-centered arrays to vertex-based for 2D structured grids"""
    ...

def array_at_faces_1d(
    array: ArrayData,
    grid: StructuredGrid,
    direction: str = 'x'
) -> np.ndarray:
    """Convert cell-centered arrays to face-centered values"""
    ...

def create_structured_grid(
    delr: ArrayData,
    delc: ArrayData,
    top: ArrayData = None,
    botm: ArrayData = None,
    **kwargs
) -> StructuredGrid:
    """Factory function for creating structured grids"""
    ...

def create_unstructured_grid(
    vertices: ArrayData,
    iverts: ArrayData,
    **kwargs
) -> UnstructuredGrid:
    """Factory function for creating unstructured grids"""
    ...

def create_vertex_grid(
    vertices: ArrayData,
    cell2d: ArrayData,
    **kwargs
) -> VertexGrid:
    """Factory function for creating vertex grids"""
    ...

Usage Examples

Basic Structured Grid

import flopy
import numpy as np

# Create simple structured grid
nlay, nrow, ncol = 3, 50, 100
delr = np.ones(ncol) * 100.0  # 100m columns
delc = np.ones(nrow) * 100.0  # 100m rows
top = 50.0
botm = [30.0, 10.0, -10.0]

# Create grid
grid = flopy.discretization.StructuredGrid(
    delr=delr,
    delc=delc,
    top=top,
    botm=botm,
    nlay=nlay,
    nrow=nrow,
    ncol=ncol
)

# Grid properties
print(f"Grid shape: {grid.shape}")
print(f"Grid extent: {grid.extent}")
print(f"Total cells: {grid.size}")

# Cell centers
xc, yc = grid.get_cellcenters()
print(f"Cell centers shape: {xc.shape}")

# Check regularity
print(f"Regular X spacing: {grid.is_regular_x()}")
print(f"Regular Y spacing: {grid.is_regular_y()}")

Advanced Structured Grid with Coordinate System

import flopy
import numpy as np

# Create grid with coordinate system and rotation
nlay, nrow, ncol = 5, 100, 200

# Variable column spacing (finer near pumping well)
delr = np.ones(ncol) * 50.0
delr[80:120] = 25.0  # Finer spacing in middle

# Variable row spacing
delc = np.ones(nrow) * 50.0
delc[40:60] = 25.0  # Finer spacing in middle

# Complex layer structure
top = 100.0
layer_thicks = [10.0, 15.0, 20.0, 25.0, 30.0]
botm = [top - sum(layer_thicks[:i+1]) for i in range(nlay)]

# Create grid with UTM coordinates and rotation
grid = flopy.discretization.StructuredGrid(
    delr=delr,
    delc=delc, 
    top=top,
    botm=botm,
    nlay=nlay,
    nrow=nrow,
    ncol=ncol,
    xoff=500000.0,  # UTM easting offset
    yoff=4000000.0,  # UTM northing offset
    angrot=15.0,  # 15 degree rotation
    epsg=32633  # UTM Zone 33N
)

# Coordinate transformations
local_x, local_y = 2500.0, 2500.0  # Local coordinates
model_x, model_y = grid.get_coords(local_x, local_y)
print(f"Model coordinates: {model_x:.1f}, {model_y:.1f}")

# Find intersecting cells
points_x = [2000.0, 3000.0, 4000.0]
points_y = [2000.0, 3000.0, 4000.0]
cells = grid.intersect(points_x, points_y)
print(f"Intersecting cells: {cells}")

# Export to shapefile
grid.write_shapefile(
    'model_grid.shp',
    array_dict={'layer_thickness': grid.delz}
)

Unstructured Grid Creation

import flopy
import numpy as np

# Create triangular mesh vertices
nvert = 100
vertices = []

# Generate random vertices in domain
np.random.seed(42)
x_coords = np.random.uniform(0, 5000, nvert)
y_coords = np.random.uniform(0, 5000, nvert)
vertices = list(zip(x_coords, y_coords))

# Define cell connectivity (simplified - would use mesh generator)
# Each cell defined by vertex indices
iverts = []
ncells = 150
for i in range(ncells):
    # Triangular cells with random vertices
    vert_indices = np.random.choice(nvert, 3, replace=False)
    iverts.append(list(vert_indices))

# Cell centers (computed from vertices)
xcenters = []
ycenters = []
for cell_verts in iverts:
    x_cell = np.mean([vertices[i][0] for i in cell_verts])
    y_cell = np.mean([vertices[i][1] for i in cell_verts])
    xcenters.append(x_cell)
    ycenters.append(y_cell)

# Layer structure
nlay = 4
top = np.ones(ncells) * 50.0
botm = np.array([30.0, 10.0, -10.0, -30.0])
botm_3d = np.repeat(botm.reshape(nlay, 1), ncells, axis=1)

# Create unstructured grid
grid = flopy.discretization.UnstructuredGrid(
    vertices=vertices,
    iverts=iverts,
    xcenters=xcenters,
    ycenters=ycenters,
    top=top,
    botm=botm_3d,
    ncpl=[ncells] * nlay
)

print(f"Unstructured grid nodes: {grid.nnodes}")
print(f"Grid extent: {grid.extent}")

# Get vertices for specific cell
cell_vertices = grid.get_cell_vertices(0)
print(f"Cell 0 vertices: {cell_vertices}")

Vertex Grid for MODFLOW 6 DISV

import flopy
import numpy as np

# Create vertex grid for DISV package
nlay = 3
ncpl = 50  # cells per layer
nvert = 75  # vertices

# Hexagonal grid vertices (simplified)
vertices = []
for i in range(nvert):
    angle = 2 * np.pi * i / nvert
    radius = 1000.0 + 500.0 * (i % 5) / 4  # Variable radius
    x = radius * np.cos(angle)
    y = radius * np.sin(angle)
    vertices.append((x, y))

# Cell definitions (cell center x, y, vertex indices)
cell2d = []
for i in range(ncpl):
    # Cell center
    angle = 2 * np.pi * i / ncpl
    x_center = 800.0 * np.cos(angle)
    y_center = 800.0 * np.sin(angle)
    
    # Connect to nearest vertices (simplified)
    vert_indices = [(i + j) % nvert for j in range(6)]
    cell_data = [i, x_center, y_center] + vert_indices
    cell2d.append(cell_data)

# Layer elevations
top = np.ones(ncpl) * 50.0
botm_layer1 = np.ones(ncpl) * 30.0
botm_layer2 = np.ones(ncpl) * 10.0  
botm_layer3 = np.ones(ncpl) * -10.0
botm = [botm_layer1, botm_layer2, botm_layer3]

# Create vertex grid
grid = flopy.discretization.VertexGrid(
    vertices=vertices,
    cell2d=cell2d,
    top=top,
    botm=botm,
    nlay=nlay,
    ncpl=ncpl,
    nvert=nvert
)

print(f"Vertex grid shape: {grid.shape}")
print(f"Cells per layer: {grid.ncpl}")
print(f"Total vertices: {grid.nvert}")

# Export for visualization
grid.write_shapefile(
    'vertex_grid.shp',
    array_dict={'top_elev': grid.top}
)

Time Discretization

import flopy
from datetime import datetime

# Create complex temporal discretization
# Stress periods: steady state, then monthly for 2 years
perioddata = []

# Initial steady state period
perioddata.append((1.0, 1, 1.0))

# Monthly periods with increasing time steps  
import calendar
for year in [2023, 2024]:
    for month in range(1, 13):
        days_in_month = calendar.monthrange(year, month)[1]
        # More time steps in first few months
        if len(perioddata) <= 6:
            nstp = 10
            tsmult = 1.2
        else:
            nstp = 5
            tsmult = 1.1
        
        perioddata.append((days_in_month, nstp, tsmult))

# Create model time object
model_time = flopy.discretization.ModelTime(
    perioddata=perioddata,
    time_units='days',
    start_datetime='2023-01-01T00:00:00'
)

print(f"Number of stress periods: {model_time.nper}")
print(f"Total simulation time: {model_time.totim[-1]} days")

# Get specific time information
totim_6months = model_time.get_totim(4, 6)  # 5th time step, 7th period
print(f"Time at 6 months: {totim_6months} days")

# Find time step for specific time
kstpkper = model_time.get_kstpkper(365.0)  # One year
print(f"Time step/period for 1 year: {kstpkper}")

# Steady state periods
steady_periods = [i for i, steady in enumerate(model_time.steady) if steady]
print(f"Steady state periods: {steady_periods}")

Common Types

# Grid and discretization types
ArrayData = Union[int, float, np.ndarray, list]
GridType = Literal['structured', 'unstructured', 'vertex']
CoordinateSystem = Union[int, str]  # EPSG code or proj4 string

# Vertex and connectivity data
VertexData = list[tuple[float, float]]
CellConnectivity = list[list[int]]
Cell2DData = list[list[Union[int, float]]]

# Temporal discretization
PeriodData = list[tuple[float, int, float]]  # (perlen, nstp, tsmult)
TimeUnits = Literal['days', 'hours', 'minutes', 'seconds', 'years']

# Spatial arrays
ElevationArray = Union[float, list[float], np.ndarray]
ThicknessArray = Union[float, list[float], np.ndarray]
SpacingArray = Union[float, list[float], np.ndarray]
DomainArray = Union[int, list[int], np.ndarray]

# Coordinate types
Coordinates = tuple[float, float]
Extent = tuple[float, float, float, float]  # xmin, xmax, ymin, ymax
CellIndex = Union[int, tuple[int, ...]]

This comprehensive documentation covers the complete discretization API including structured, unstructured, and vertex grids, temporal discretization, and utility functions. The examples demonstrate basic to advanced grid creation scenarios with coordinate systems, variable spacing, and complex geometries.

Install with Tessl CLI

npx tessl i tessl/pypi-flopy

docs

discretization.md

export.md

file-io.md

index.md

modflow6.md

modflow2005.md

particle-tracking.md

plotting.md

transport.md

utilities.md

tile.json