0
# Visualization and Plotting
1

2
This module provides comprehensive plotting utilities for map views, cross-sections, and 3D visualization with matplotlib integration and export capabilities. These tools enable visualization of model grids, input arrays, results, and analysis outputs for both structured and unstructured grids.
3

4
## Map View Plotting
5

6
### PlotMapView
7

8
Main class for creating map view plots of MODFLOW models and results.
9

10
```python { .api }
11
class PlotMapView:
12
    """Map view plotting for MODFLOW models"""
13
    def __init__(
14
        self,
15
        model: object = None,
16
        ax: object = None,
17
        layer: int = 0,
18
        extent: list[float] = None,
19
        modelgrid: object = None,
20
        **kwargs
21
    ): ...
22
    
23
    def plot_grid(
24
        self,
25
        **kwargs
26
    ) -> object:
27
        """Plot model grid"""
28
        ...
29
    
30
    def plot_array(
31
        self,
32
        a: np.ndarray,
33
        masked_values: list = None,
34
        **kwargs
35
    ) -> object:
36
        """Plot 2D array data"""
37
        ...
38
    
39
    def contour_array(
40
        self,
41
        a: np.ndarray,
42
        masked_values: list = None,
43
        tri_mask: bool = False,
44
        **kwargs
45
    ) -> object:
46
        """Create contour plot of array data"""
47
        ...
48
    
49
    def plot_ibound(
50
        self,
51
        ibound: np.ndarray = None,
52
        color_noflow: str = 'black',
53
        color_ch: str = 'blue',
54
        color_vpt: str = 'red',
55
        color_active: str = None,
56
        **kwargs
57
    ) -> object:
58
        """Plot ibound array"""
59
        ...
60
    
61
    def plot_inactive(
62
        self,
63
        ibound: np.ndarray = None,
64
        color_noflow: str = 'black',
65
        **kwargs
66
    ) -> object:
67
        """Plot inactive cells"""
68
        ...
69
    
70
    def plot_bc(
71
        self,
72
        name: str = None,
73
        package: object = None,
74
        kper: int = 0,
75
        color: str = None,
76
        plotAll: bool = False,
77
        boundname: str = None,
78
        **kwargs
79
    ) -> object:
80
        """Plot boundary conditions"""
81
        ...
82
    
83
    def plot_shapefile(
84
        self,
85
        shp,
86
        **kwargs
87
    ) -> object:
88
        """Plot shapefile data"""
89
        ...
90
    
91
    def plot_shapes(
92
        self,
93
        obj,
94
        **kwargs
95
    ) -> object:
96
        """Plot geometric shapes"""
97
        ...
98
    
99
    def plot_centers(
100
        self,
101
        **kwargs
102
    ) -> object:
103
        """Plot cell centers"""
104
        ...
105
    
106
    def plot_pathline(
107
        self,
108
        pl: list = None,
109
        travel_time: np.ndarray = None,
110
        **kwargs
111
    ) -> object:
112
        """Plot MODPATH pathlines"""
113
        ...
114
    
115
    def plot_endpoint(
116
        self,
117
        ep: np.ndarray = None,
118
        direction: str = 'ending',
119
        selection: np.ndarray = None,
120
        **kwargs
121
    ) -> object:
122
        """Plot MODPATH endpoints"""
123
        ...
124
    
125
    def plot_timeseries(
126
        self,
127
        ts: np.ndarray = None,
128
        travel_time: float = None,
129
        **kwargs
130
    ) -> object:
131
        """Plot MODPATH timeseries data"""
132
        ...
133
    
134
    def plot_vector(
135
        self,
136
        vx: np.ndarray,
137
        vy: np.ndarray,
138
        istep: int = 1,
139
        jstep: int = 1,
140
        normalize: bool = False,
141
        masked_values: list = None,
142
        **kwargs
143
    ) -> object:
144
        """Plot vector field"""
145
        ...
146
    
147
    @property
148
    def extent(self) -> tuple[float, float, float, float]:
149
        """Plot extent (xmin, xmax, ymin, ymax)"""
150
        ...
151
    
152
    @property
153
    def mg(self) -> object:
154
        """Model grid object"""
155
        ...
156
```
157

158
## Cross-Section Plotting
159

160
### PlotCrossSection
161

162
Class for creating cross-sectional plots through MODFLOW models.
163

164
```python { .api }
165
class PlotCrossSection:
166
    """Cross-section plotting for MODFLOW models"""
167
    def __init__(
168
        self,
169
        model: object = None,
170
        ax: object = None,
171
        line: dict = None,
172
        extent: list[float] = None,
173
        modelgrid: object = None,
174
        geographic_coords: bool = False,
175
        **kwargs
176
    ): ...
177
    
178
    def plot_grid(
179
        self,
180
        **kwargs
181
    ) -> object:
182
        """Plot cross-section grid"""
183
        ...
184
    
185
    def plot_array(
186
        self,
187
        a: np.ndarray,
188
        masked_values: list = None,
189
        **kwargs
190
    ) -> object:
191
        """Plot array data in cross-section"""
192
        ...
193
    
194
    def contour_array(
195
        self,
196
        a: np.ndarray,
197
        masked_values: list = None,
198
        **kwargs
199
    ) -> object:
200
        """Create contour plot in cross-section"""
201
        ...
202
    
203
    def plot_surface(
204
        self,
205
        a: np.ndarray,
206
        masked_values: list = None,
207
        **kwargs
208
    ) -> object:
209
        """Plot surface elevation"""
210
        ...
211
    
212
    def plot_fill_between(
213
        self,
214
        a: np.ndarray,
215
        colors: list = None,
216
        masked_values: list = None,
217
        **kwargs
218
    ) -> object:
219
        """Plot filled areas between surfaces"""
220
        ...
221
    
222
    def plot_bc(
223
        self,
224
        name: str = None,
225
        package: object = None,
226
        kper: int = 0,
227
        color: str = 'red',
228
        **kwargs
229
    ) -> object:
230
        """Plot boundary conditions in cross-section"""
231
        ...
232
    
233
    def plot_ibound(
234
        self,
235
        ibound: np.ndarray = None,
236
        color_noflow: str = 'black',
237
        color_ch: str = 'blue',
238
        **kwargs
239
    ) -> object:
240
        """Plot ibound in cross-section"""
241
        ...
242
    
243
    def plot_discharge(
244
        self,
245
        frf: np.ndarray,
246
        fff: np.ndarray,
247
        flf: np.ndarray = None,
248
        head: np.ndarray = None,
249
        normalize: bool = False,
250
        **kwargs
251
    ) -> object:
252
        """Plot discharge vectors in cross-section"""
253
        ...
254
    
255
    def plot_pathline(
256
        self,
257
        pl: list = None,
258
        travel_time: np.ndarray = None,
259
        **kwargs
260
    ) -> object:
261
        """Plot pathlines in cross-section"""
262
        ...
263
    
264
    def plot_endpoint(
265
        self,
266
        ep: np.ndarray = None,
267
        direction: str = 'ending',
268
        **kwargs
269
    ) -> object:
270
        """Plot endpoints in cross-section"""
271
        ...
272
    
273
    @property
274
    def extent(self) -> tuple[float, float, float, float]:
275
        """Cross-section extent"""
276
        ...
277
    
278
    @property
279
    def direction(self) -> str:
280
        """Cross-section direction"""
281
        ...
282
```
283

284
## Utility Classes
285

286
### PlotUtilities
287

288
General plotting utilities and helper functions.
289

290
```python { .api }
291
class PlotUtilities:
292
    """General plotting utilities"""
293
    
294
    @staticmethod
295
    def plot_text(
296
        ax: object,
297
        model: object = None,
298
        text_dict: dict = None,
299
        **kwargs
300
    ) -> None:
301
        """Add text to plot"""
302
        ...
303
    
304
    @staticmethod
305
    def plot_shapefile(
306
        shp: str,
307
        ax: object = None,
308
        **kwargs
309
    ) -> object:
310
        """Plot shapefile data"""
311
        ...
312
    
313
    @staticmethod
314
    def add_colorbar(
315
        mappable: object,
316
        ax: object = None,
317
        **kwargs
318
    ) -> object:
319
        """Add colorbar to plot"""
320
        ...
321
    
322
    @staticmethod
323
    def set_aspect(
324
        ax: object,
325
        aspect: str = 'equal',
326
        **kwargs
327
    ) -> None:
328
        """Set plot aspect ratio"""
329
        ...
330
```
331

332
### SwiConcentration
333

334
Specialized plotting for seawater intrusion concentration data.
335

336
```python { .api }
337
class SwiConcentration:
338
    """Seawater intrusion concentration plotting"""
339
    def __init__(
340
        self,
341
        model: object = None,
342
        botm: np.ndarray = None,
343
        istrat: int = 1,
344
        nu: np.ndarray = None,
345
        **kwargs
346
    ): ...
347
    
348
    def plot_concentration(
349
        self,
350
        zeta: np.ndarray,
351
        layer: int = None,
352
        **kwargs
353
    ) -> object:
354
        """Plot concentration distribution"""
355
        ...
356
    
357
    def plot_zeta_surfaces(
358
        self,
359
        zeta: np.ndarray,
360
        **kwargs
361
    ) -> object:
362
        """Plot zeta surfaces"""
363
        ...
364
```
365

366
## Standalone Plotting Functions
367

368
```python { .api }
369
def plot_shapefile(
370
    shp: str,
371
    ax: object = None,
372
    layer: str = None,
373
    **kwargs
374
) -> object:
375
    """Plot shapefile data on map"""
376
    ...
377

378
def shapefile_extents(shp: str) -> tuple[float, float, float, float]:
379
    """Get shapefile spatial extents"""
380
    ...
381

382
def plot_cvfd(
383
    verts: np.ndarray,
384
    iverts: list,
385
    xc: np.ndarray = None,
386
    yc: np.ndarray = None,
387
    head: np.ndarray = None,
388
    **kwargs
389
) -> object:
390
    """Plot control volume finite difference grid"""
391
    ...
392
```
393

394
## Styling and Formatting
395

396
### Plotting Styles
397

398
Predefined plotting styles and color schemes.
399

400
```python { .api }
401
# Style constants and color maps
402
styles = {
403
    'default': {
404
        'grid_color': 'black',
405
        'grid_linewidth': 0.5,
406
        'inactive_color': 'gray',
407
        'bc_colors': {
408
            'WEL': 'red',
409
            'DRN': 'yellow', 
410
            'RIV': 'blue',
411
            'GHB': 'cyan',
412
            'CHD': 'navy'
413
        }
414
    },
415
    'publication': {
416
        'grid_color': 'black',
417
        'grid_linewidth': 0.3,
418
        'font_size': 12,
419
        'title_size': 14
420
    },
421
    'presentation': {
422
        'grid_color': 'gray',
423
        'grid_linewidth': 1.0,
424
        'font_size': 14,
425
        'title_size': 18
426
    }
427
}
428

429
def set_style(style_name: str) -> None:
430
    """Set plotting style"""
431
    ...
432

433
def get_color_palette(name: str) -> list[str]:
434
    """Get predefined color palette"""
435
    ...
436
```
437

438
## Usage Examples
439

440
### Basic Model Visualization
441

442
```python
443
import flopy
444
import flopy.utils as fpu
445
import matplotlib.pyplot as plt
446
import numpy as np
447

448
# Load or create MODFLOW model
449
mf = flopy.modflow.Modflow.load('model.nam')
450

451
# Create map view plot
452
fig, ax = plt.subplots(figsize=(12, 10))
453
mapview = flopy.plot.PlotMapView(model=mf, ax=ax, layer=0)
454

455
# Plot model grid
456
mapview.plot_grid(lw=0.5, color='gray')
457

458
# Plot boundary conditions
459
mapview.plot_bc('WEL', color='red', s=50)  # Wells
460
mapview.plot_bc('RIV', color='blue', lw=3)  # Rivers
461
mapview.plot_bc('GHB', color='cyan', lw=2)  # General head boundaries
462

463
# Plot head results if available
464
try:
465
    hds = fpu.HeadFile('model.hds')
466
    head = hds.get_data(kstpkper=(0, 0))
467
    
468
    # Contour heads
469
    contours = mapview.contour_array(
470
        head,
471
        levels=np.linspace(head.min(), head.max(), 20),
472
        colors='black',
473
        linewidths=1
474
    )
475
    ax.clabel(contours, inline=True, fontsize=8, fmt='%.1f')
476
    
477
    # Fill contours with colors
478
    filled = mapview.plot_array(
479
        head,
480
        cmap='viridis',
481
        alpha=0.7
482
    )
483
    cb = plt.colorbar(filled, ax=ax, shrink=0.8)
484
    cb.set_label('Head Elevation (m)', fontsize=12)
485
    
486
    hds.close()
487
    
488
except FileNotFoundError:
489
    print("Head file not found, skipping head contours")
490

491
# Plot inactive cells
492
if hasattr(mf, 'bas6'):
493
    mapview.plot_inactive(ibound=mf.bas6.ibound.array, color_noflow='gray')
494

495
# Add title and labels
496
ax.set_title('MODFLOW Model - Layer 1 Heads', fontsize=14, fontweight='bold')
497
ax.set_xlabel('X Coordinate (m)', fontsize=12)
498
ax.set_ylabel('Y Coordinate (m)', fontsize=12)
499

500
# Set equal aspect ratio
501
ax.set_aspect('equal')
502

503
plt.tight_layout()
504
plt.show()
505
```
506

507
### Flow Vector Visualization
508

509
```python
510
import flopy
511
import flopy.utils as fpu
512
import matplotlib.pyplot as plt
513
import numpy as np
514

515
# Load model and results
516
mf = flopy.modflow.Modflow.load('model.nam')
517

518
# Read budget file for flow data
519
cbb = fpu.CellBudgetFile('model.cbb')
520
frf = cbb.get_data(text='FLOW RIGHT FACE')[0]  # Flow right face
521
fff = cbb.get_data(text='FLOW FRONT FACE')[0]  # Flow front face
522
flf = cbb.get_data(text='FLOW LOWER FACE')[0]  # Flow lower face (optional)
523

524
# Read head data
525
hds = fpu.HeadFile('model.hds')
526
head = hds.get_data(kstpkper=(0, 0))
527

528
# Create map view
529
fig, ax = plt.subplots(figsize=(14, 10))
530
mapview = flopy.plot.PlotMapView(model=mf, ax=ax, layer=0)
531

532
# Plot head contours as background
533
head_contours = mapview.contour_array(
534
    head,
535
    levels=15,
536
    colors='gray',
537
    alpha=0.5,
538
    linewidths=0.8
539
)
540

541
# Plot discharge vectors
542
discharge_plot = mapview.plot_discharge(
543
    frf, fff, flf,
544
    head=head,
545
    normalize=True,  # Normalize vectors for better visualization
546
    color='blue',
547
    scale=50,  # Scale factor for arrow size
548
    headwidth=3,
549
    headlength=2,
550
    width=0.002
551
)
552

553
# Plot model grid (light)
554
mapview.plot_grid(lw=0.3, color='lightgray', alpha=0.5)
555

556
# Plot boundary conditions
557
mapview.plot_bc('WEL', color='red', marker='o', s=100, label='Wells')
558
mapview.plot_bc('RIV', color='cyan', lw=4, alpha=0.7, label='Rivers')
559

560
# Add legend
561
ax.legend(loc='upper right')
562

563
ax.set_title('Groundwater Flow Vectors', fontsize=16, fontweight='bold')
564
ax.set_xlabel('X Coordinate (m)', fontsize=12)
565
ax.set_ylabel('Y Coordinate (m)', fontsize=12)
566

567
plt.tight_layout()
568
plt.show()
569

570
# Close files
571
cbb.close()
572
hds.close()
573
```
574

575
### Cross-Section Visualization
576

577
```python
578
import flopy
579
import flopy.utils as fpu
580
import matplotlib.pyplot as plt
581
import numpy as np
582

583
# Load model
584
mf = flopy.modflow.Modflow.load('model.nam')
585

586
# Define cross-section line
587
line = {'line': [(0, 2500), (5000, 2500)]}  # West-east cross-section
588

589
# Create cross-section plot
590
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(14, 12))
591

592
# First subplot - Head distribution
593
xs1 = flopy.plot.PlotCrossSection(model=mf, ax=ax1, line=line)
594

595
# Read and plot heads
596
try:
597
    hds = fpu.HeadFile('model.hds')
598
    head = hds.get_data(kstpkper=(0, 0))
599
    
600
    # Plot head array
601
    head_plot = xs1.plot_array(
602
        head,
603
        cmap='RdYlBu',
604
        alpha=0.8
605
    )
606
    
607
    # Add colorbar
608
    cb1 = plt.colorbar(head_plot, ax=ax1, shrink=0.8)
609
    cb1.set_label('Head Elevation (m)', fontsize=12)
610
    
611
    # Contour heads
612
    head_contours = xs1.contour_array(
613
        head,
614
        levels=10,
615
        colors='black',
616
        linewidths=1
617
    )
618
    ax1.clabel(head_contours, inline=True, fontsize=8)
619
    
620
    hds.close()
621

622
except FileNotFoundError:
623
    print("Head file not found")
624

625
# Plot grid and model layers
626
xs1.plot_grid(lw=1, color='black')
627

628
# Plot ibound
629
if hasattr(mf, 'bas6'):
630
    xs1.plot_ibound(
631
        ibound=mf.bas6.ibound.array,
632
        color_noflow='gray',
633
        alpha=0.5
634
    )
635

636
ax1.set_title('Cross-Section: Head Distribution', fontsize=14, fontweight='bold')
637
ax1.set_ylabel('Elevation (m)', fontsize=12)
638

639
# Second subplot - Hydraulic conductivity
640
xs2 = flopy.plot.PlotCrossSection(model=mf, ax=ax2, line=line)
641

642
# Plot hydraulic conductivity
643
if hasattr(mf, 'lpf'):
644
    k_array = mf.lpf.hk.array
645
    k_plot = xs2.plot_array(
646
        k_array,
647
        cmap='plasma',
648
        norm=plt.Normalize(vmin=k_array.min(), vmax=k_array.max()),
649
        alpha=0.8
650
    )
651
    
652
    cb2 = plt.colorbar(k_plot, ax=ax2, shrink=0.8)
653
    cb2.set_label('Hydraulic Conductivity (m/d)', fontsize=12)
654

655
# Plot grid
656
xs2.plot_grid(lw=1, color='black')
657

658
# Plot boundary conditions in cross-section
659
xs2.plot_bc('WEL', color='red', s=100)
660

661
ax2.set_title('Cross-Section: Hydraulic Conductivity', fontsize=14, fontweight='bold')
662
ax2.set_xlabel('Distance along section (m)', fontsize=12)
663
ax2.set_ylabel('Elevation (m)', fontsize=12)
664

665
plt.tight_layout()
666
plt.show()
667
```
668

669
### Particle Tracking Visualization
670

671
```python
672
import flopy
673
import flopy.utils as fpu
674
import matplotlib.pyplot as plt
675
import numpy as np
676

677
# Load model
678
mf = flopy.modflow.Modflow.load('model.nam')
679

680
# Create map view
681
fig, ax = plt.subplots(figsize=(14, 10))
682
mapview = flopy.plot.PlotMapView(model=mf, ax=ax, layer=0)
683

684
# Plot base head contours
685
try:
686
    hds = fpu.HeadFile('model.hds')
687
    head = hds.get_data(kstpkper=(0, 0))
688
    
689
    contours = mapview.contour_array(
690
        head,
691
        levels=15,
692
        colors='lightgray',
693
        linewidths=1,
694
        alpha=0.7
695
    )
696
    hds.close()
697
except:
698
    pass
699

700
# Plot model grid (light)
701
mapview.plot_grid(lw=0.3, color='lightgray', alpha=0.5)
702

703
# Plot MODPATH results
704
try:
705
    # Read pathline data
706
    pth_file = fpu.PathlineFile('model.mppth')
707
    pathlines = pth_file.get_alldata()
708
    
709
    # Plot pathlines colored by travel time
710
    pathline_plot = mapview.plot_pathline(
711
        pathlines,
712
        travel_time=None,  # Will use time from pathline data
713
        lw=2,
714
        alpha=0.8,
715
        cmap='viridis'
716
    )
717
    
718
    print(f"Plotted {len(pathlines)} pathlines")
719
    
720
    # Read endpoint data
721
    ep_file = fpu.EndpointFile('model.mpend')
722
    endpoints = ep_file.get_alldata()
723
    
724
    # Plot starting points (green) and ending points (red)
725
    starting_points = mapview.plot_endpoint(
726
        endpoints,
727
        direction='starting',
728
        marker='o',
729
        color='green',
730
        s=50,
731
        label='Starting points',
732
        alpha=0.8
733
    )
734
    
735
    ending_points = mapview.plot_endpoint(
736
        endpoints,
737
        direction='ending',
738
        marker='s',
739
        color='red',
740
        s=30,
741
        label='Ending points',
742
        alpha=0.8
743
    )
744
    
745
    print(f"Plotted {len(endpoints)} particle endpoints")
746
    
747
except FileNotFoundError as e:
748
    print(f"MODPATH files not found: {e}")
749

750
# Plot boundary conditions
751
mapview.plot_bc('WEL', color='blue', marker='o', s=150, 
752
               edgecolor='darkblue', linewidth=2, label='Wells')
753
mapview.plot_bc('RIV', color='cyan', lw=5, alpha=0.7, label='Rivers')
754

755
# Add colorbar for pathlines if they were plotted
756
try:
757
    if 'pathline_plot' in locals():
758
        cb = plt.colorbar(pathline_plot, ax=ax, shrink=0.8)
759
        cb.set_label('Travel Time (days)', fontsize=12)
760
except:
761
    pass
762

763
# Legend
764
ax.legend(loc='upper right', fontsize=10)
765

766
ax.set_title('Particle Tracking Results', fontsize=16, fontweight='bold')
767
ax.set_xlabel('X Coordinate (m)', fontsize=12)
768
ax.set_ylabel('Y Coordinate (m)', fontsize=12)
769
ax.set_aspect('equal')
770

771
plt.tight_layout()
772
plt.show()
773
```
774

775
### Multi-Layer Visualization
776

777
```python
778
import flopy
779
import flopy.utils as fpu
780
import matplotlib.pyplot as plt
781
import numpy as np
782

783
# Load model
784
mf = flopy.modflow.Modflow.load('model.nam')
785
nlay = mf.dis.nlay
786

787
# Create multi-panel figure
788
fig, axes = plt.subplots(2, 2, figsize=(16, 12))
789
axes = axes.flatten()
790

791
# Read head data
792
try:
793
    hds = fpu.HeadFile('model.hds')
794
    head = hds.get_data(kstpkper=(0, 0))
795
    hds.close()
796
    
797
    # Global head range for consistent coloring
798
    vmin, vmax = head.min(), head.max()
799
    levels = np.linspace(vmin, vmax, 20)
800
    
801
except FileNotFoundError:
802
    head = None
803
    vmin = vmax = None
804
    levels = None
805

806
# Plot each layer
807
for i in range(min(nlay, 4)):  # Up to 4 layers
808
    ax = axes[i]
809
    mapview = flopy.plot.PlotMapView(model=mf, ax=ax, layer=i)
810
    
811
    # Plot grid
812
    mapview.plot_grid(lw=0.3, color='black', alpha=0.5)
813
    
814
    # Plot head data if available
815
    if head is not None:
816
        head_plot = mapview.plot_array(
817
            head,
818
            cmap='RdYlBu',
819
            vmin=vmin,
820
            vmax=vmax,
821
            alpha=0.8
822
        )
823
        
824
        # Contour lines
825
        contours = mapview.contour_array(
826
            head,
827
            levels=levels[::2],  # Every other contour
828
            colors='black',
829
            linewidths=0.8,
830
            alpha=0.7
831
        )
832
    
833
    # Plot boundary conditions
834
    mapview.plot_bc('WEL', color='red', s=30, alpha=0.8)
835
    mapview.plot_bc('RIV', color='blue', lw=2, alpha=0.7)
836
    
837
    # Plot inactive cells
838
    if hasattr(mf, 'bas6'):
839
        mapview.plot_inactive(
840
            ibound=mf.bas6.ibound.array,
841
            color_noflow='lightgray'
842
        )
843
    
844
    ax.set_title(f'Layer {i+1}', fontsize=12, fontweight='bold')
845
    ax.set_xlabel('X (m)', fontsize=10)
846
    ax.set_ylabel('Y (m)', fontsize=10)
847
    ax.set_aspect('equal')
848

849
# Add shared colorbar
850
if head is not None and 'head_plot' in locals():
851
    # Create colorbar axis
852
    cbar_ax = fig.add_axes([0.92, 0.15, 0.02, 0.7])
853
    cb = plt.colorbar(head_plot, cax=cbar_ax)
854
    cb.set_label('Head Elevation (m)', fontsize=12, rotation=270, labelpad=20)
855

856
fig.suptitle('Multi-Layer Head Distribution', fontsize=16, fontweight='bold')
857
plt.tight_layout()
858
plt.subplots_adjust(right=0.9)
859
plt.show()
860
```
861

862
### Advanced Concentration Plume Visualization
863

864
```python
865
import flopy
866
import flopy.utils as fpu
867
import matplotlib.pyplot as plt
868
import numpy as np
869
from matplotlib.colors import LogNorm
870

871
# Load MT3DMS results
872
try:
873
    ucn = fpu.UcnFile('MT3D001.UCN')
874
    times = ucn.get_times()
875
    
876
    # Create time series plot
877
    fig = plt.figure(figsize=(16, 12))
878
    
879
    # Select representative times
880
    time_indices = [0, len(times)//4, len(times)//2, -1]
881
    
882
    for idx, time_idx in enumerate(time_indices):
883
        ax = plt.subplot(2, 2, idx + 1)
884
        
885
        # Get concentration data
886
        conc = ucn.get_data(idx=time_idx)
887
        time_val = times[time_idx]
888
        
889
        # Create map view (assuming structured grid)
890
        # Note: Would need model object for proper grid, using generic approach
891
        
892
        # Plot concentration with log scale for better visualization
893
        conc_masked = np.ma.masked_where(conc <= 0, conc)
894
        
895
        im = ax.imshow(
896
            conc_masked[0, :, :],  # Top layer
897
            cmap='plasma',
898
            norm=LogNorm(vmin=0.01, vmax=conc.max()),
899
            origin='lower',
900
            alpha=0.8
901
        )
902
        
903
        # Add contour lines
904
        contour_levels = np.logspace(-2, np.log10(conc.max()), 8)
905
        cs = ax.contour(
906
            conc[0, :, :],
907
            levels=contour_levels,
908
            colors='black',
909
            linewidths=1,
910
            alpha=0.6
911
        )
912
        ax.clabel(cs, inline=True, fontsize=8, fmt='%.2f')
913
        
914
        ax.set_title(f'Time = {time_val:.1f} days', fontsize=12, fontweight='bold')
915
        ax.set_xlabel('Column', fontsize=10)
916
        ax.set_ylabel('Row', fontsize=10)
917
        ax.set_aspect('equal')
918
    
919
    # Add shared colorbar
920
    cbar_ax = fig.add_axes([0.92, 0.15, 0.02, 0.7])
921
    cb = plt.colorbar(im, cax=cbar_ax)
922
    cb.set_label('Concentration (mg/L)', fontsize=12, rotation=270, labelpad=20)
923
    
924
    fig.suptitle('Contaminant Plume Evolution', fontsize=16, fontweight='bold')
925
    plt.tight_layout()
926
    plt.subplots_adjust(right=0.9)
927
    plt.show()
928
    
929
    ucn.close()
930

931
except FileNotFoundError:
932
    print("MT3DMS concentration file not found")
933

934
# Concentration breakthrough curves
935
try:
936
    ucn = fpu.UcnFile('MT3D001.UCN')
937
    
938
    # Define monitoring points
939
    monitor_points = [(0, 10, 25), (0, 20, 50), (0, 30, 75)]
940
    
941
    # Get time series
942
    conc_ts = ucn.get_ts(idx=monitor_points)
943
    times = ucn.get_times()
944
    
945
    # Plot breakthrough curves
946
    plt.figure(figsize=(12, 8))
947
    
948
    for i, (lay, row, col) in enumerate(monitor_points):
949
        plt.semilogy(
950
            times,
951
            conc_ts[:, i],
952
            'o-',
953
            linewidth=2,
954
            markersize=4,
955
            label=f'MW-{i+1} (R{row}, C{col})'
956
        )
957
    
958
    plt.xlabel('Time (days)', fontsize=12)
959
    plt.ylabel('Concentration (mg/L)', fontsize=12)
960
    plt.title('Concentration Breakthrough Curves', fontsize=14, fontweight='bold')
961
    plt.legend(fontsize=10)
962
    plt.grid(True, alpha=0.3)
963
    plt.tight_layout()
964
    plt.show()
965
    
966
    ucn.close()
967

968
except FileNotFoundError:
969
    print("MT3DMS concentration file not found")
970
```
971

972
### Custom Styling and Publication-Quality Plots
973

974
```python
975
import flopy
976
import matplotlib.pyplot as plt
977
import matplotlib as mpl
978
import numpy as np
979

980
# Set publication style
981
plt.style.use('seaborn-v0_8-paper')  # Clean, publication-ready style
982
mpl.rcParams['font.size'] = 11
983
mpl.rcParams['axes.labelsize'] = 12
984
mpl.rcParams['axes.titlesize'] = 14
985
mpl.rcParams['xtick.labelsize'] = 10
986
mpl.rcParams['ytick.labelsize'] = 10
987
mpl.rcParams['legend.fontsize'] = 10
988
mpl.rcParams['figure.titlesize'] = 16
989

990
# Custom color schemes
991
custom_colors = {
992
    'wells': '#d62728',      # Red
993
    'rivers': '#1f77b4',     # Blue
994
    'ghb': '#ff7f0e',        # Orange
995
    'drn': '#2ca02c',        # Green
996
    'inactive': '#7f7f7f',   # Gray
997
    'grid': '#cccccc'        # Light gray
998
}
999

1000
# Load model
1001
mf = flopy.modflow.Modflow.load('model.nam')
1002

1003
# Create publication figure
1004
fig, ax = plt.subplots(figsize=(10, 8), dpi=300)  # High DPI for publication
1005

1006
# Create map view
1007
mapview = flopy.plot.PlotMapView(model=mf, ax=ax, layer=0)
1008

1009
# Plot with custom styling
1010
mapview.plot_grid(
1011
    lw=0.5,
1012
    color=custom_colors['grid'],
1013
    alpha=0.7
1014
)
1015

1016
# Plot boundary conditions with custom styling
1017
mapview.plot_bc(
1018
    'WEL',
1019
    color=custom_colors['wells'],
1020
    marker='o',
1021
    s=80,
1022
    edgecolor='darkred',
1023
    linewidth=1.5,
1024
    label='Pumping Wells',
1025
    zorder=5
1026
)
1027

1028
mapview.plot_bc(
1029
    'RIV',
1030
    color=custom_colors['rivers'],
1031
    lw=3,
1032
    alpha=0.8,
1033
    label='Rivers',
1034
    zorder=4
1035
)
1036

1037
mapview.plot_bc(
1038
    'GHB',
1039
    color=custom_colors['ghb'],
1040
    lw=2,
1041
    alpha=0.8,
1042
    label='General Head Boundary',
1043
    zorder=3
1044
)
1045

1046
# Plot head data with custom colormap
1047
try:
1048
    hds = fpu.HeadFile('model.hds')
1049
    head = hds.get_data(kstpkper=(0, 0))
1050
    
1051
    # Use a perceptually uniform colormap
1052
    head_plot = mapview.plot_array(
1053
        head,
1054
        cmap='viridis',
1055
        alpha=0.7,
1056
        zorder=1
1057
    )
1058
    
1059
    # Clean contours
1060
    contours = mapview.contour_array(
1061
        head,
1062
        levels=np.linspace(head.min(), head.max(), 15),
1063
        colors='white',
1064
        linewidths=0.8,
1065
        alpha=0.9,
1066
        zorder=2
1067
    )
1068
    
1069
    # Format colorbar
1070
    cb = plt.colorbar(
1071
        head_plot,
1072
        ax=ax,
1073
        shrink=0.8,
1074
        aspect=20,
1075
        pad=0.02
1076
    )
1077
    cb.set_label('Hydraulic Head (m above MSL)', fontsize=12, labelpad=15)
1078
    cb.ax.tick_params(labelsize=10)
1079
    
1080
    hds.close()
1081
    
1082
except FileNotFoundError:
1083
    pass
1084

1085
# Plot inactive cells
1086
if hasattr(mf, 'bas6'):
1087
    mapview.plot_inactive(
1088
        ibound=mf.bas6.ibound.array,
1089
        color_noflow=custom_colors['inactive'],
1090
        alpha=0.3
1091
    )
1092

1093
# Professional formatting
1094
ax.set_xlabel('Easting (m)', fontsize=12, fontweight='bold')
1095
ax.set_ylabel('Northing (m)', fontsize=12, fontweight='bold')
1096
ax.set_title('Groundwater Flow Model - Steady State Heads', 
1097
             fontsize=14, fontweight='bold', pad=20)
1098

1099
# Clean legend
1100
legend = ax.legend(
1101
    loc='upper left',
1102
    frameon=True,
1103
    fancybox=True,
1104
    shadow=True,
1105
    fontsize=10
1106
)
1107
legend.get_frame().set_facecolor('white')
1108
legend.get_frame().set_alpha(0.9)
1109

1110
# Set aspect ratio and tight layout
1111
ax.set_aspect('equal')
1112

1113
# Add north arrow and scale bar (custom functions would be needed)
1114
# add_north_arrow(ax, xy=(0.9, 0.9), size=0.05)
1115
# add_scale_bar(ax, length=1000, location='lower right')
1116

1117
# Grid and ticks
1118
ax.grid(True, alpha=0.3, linestyle='--', linewidth=0.5)
1119
ax.tick_params(axis='both', which='major', labelsize=10)
1120

1121
# Tight layout with padding
1122
plt.tight_layout(pad=1.5)
1123

1124
# Save in multiple formats for publication
1125
plt.savefig('model_heads.png', dpi=300, bbox_inches='tight', 
1126
            facecolor='white', edgecolor='none')
1127
plt.savefig('model_heads.pdf', bbox_inches='tight', 
1128
            facecolor='white', edgecolor='none')
1129
plt.savefig('model_heads.svg', bbox_inches='tight', 
1130
            facecolor='white', edgecolor='none')
1131

1132
plt.show()
1133
```
1134

1135
## Common Types
1136

1137
```python { .api }
1138
# Plotting types
1139
PlotArray = np.ndarray  # 2D or 3D array for plotting
1140
ColorMap = Union[str, mpl.colors.Colormap]
1141
PlotExtent = tuple[float, float, float, float]  # (xmin, xmax, ymin, ymax)
1142

1143
# Line and geometry types
1144
CrossSectionLine = dict[str, list[tuple[float, float]]]
1145
PlotLine = dict[str, Union[list, np.ndarray]]
1146
Coordinates = tuple[float, float]
1147

1148
# Style and formatting types
1149
PlotStyle = dict[str, Union[str, int, float]]
1150
ColorSpec = Union[str, tuple[float, float, float]]
1151
MarkerStyle = str
1152
LineStyle = str
1153

1154
# Vector data types
1155
VectorField = tuple[np.ndarray, np.ndarray]  # (vx, vy)
1156
FlowData = tuple[np.ndarray, np.ndarray, np.ndarray]  # (frf, fff, flf)
1157
DischargeData = tuple[np.ndarray, np.ndarray, np.ndarray]  # (qx, qy, qz)
1158

1159
# Particle tracking types
1160
PathlineData = list[np.ndarray]
1161
EndpointData = np.ndarray
1162
TimeseriesData = np.ndarray
1163

1164
# 3D plotting types
1165
SurfaceData = np.ndarray  # 2D array for surface plots
1166
VoxelData = np.ndarray   # 3D boolean array for voxel plots
1167
ScatterData = tuple[np.ndarray, np.ndarray, np.ndarray]  # (x, y, z)
1168

1169
# File and output types
1170
ImageFormat = Literal['png', 'pdf', 'svg', 'eps', 'jpg', 'tiff']
1171
PlotSize = tuple[float, float]  # (width, height) in inches
1172
DPI = int  # Dots per inch for raster output
1173
```
1174

1175
This comprehensive documentation covers the complete plotting and visualization API for FloPy including map views, cross-sections, 3D plots, and specialized visualizations. The examples demonstrate basic to advanced plotting scenarios including publication-quality figures, multi-layer visualization, and particle tracking results.